Line geometry support in st_extract

[kadyb]

Could you consider support for line geometry in the st_extract function? I compared it with the raster package and it is more easier and user-friendly, but maybe there is some other trick in stars? Thanks.

library("sf")
library("stars")
library("raster")

file = system.file("tif/olinda_dem_utm25s.tif", package = "stars")
ras_stars = read_stars(file)
ras_raster = raster(file)


pts_matrix = matrix(c(291511, 9113326,
                      293936, 9117740),
                    ncol = 2, byrow = TRUE)


# raster
line = spLines(pts_matrix, crs = crs(ras_raster))
elev_raster = extract(ras_raster, line)[[1]]

length(elev_raster)
#> 77
summary(elev_raster)
#> Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#> 6.00   15.00   25.00   26.88   41.00   55.00


# stars
line = st_linestring(pts_matrix)
line = st_sfc(line, crs = st_crs(ras_stars))
sample_line = st_line_sample(line, density = 1/90)
sample_line = st_cast(sample_line, "POINT")
elev_stars = st_extract(ras_stars, sample_line)
elev_stars = st_drop_geometry(elev_stars)[, 1]

length(elev_stars)
#> 56
summary(elev_stars)
#> Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#> 6.00   13.75   24.50   26.82   41.25   55.00

[edzer]

Yes, we could do that, the question is rather how:

• since we want to be able to do this for more than one line, I guess a list column with raster values (similar to what raster returns now, but than as a column in an sf object) for each line geometry should come back, possibly a stars object if the object has multiple layers and/or time instances
• a default value for the density needs to be chosen. What does raster do here? I noted that the number of points is not identical.
• similar question to st_extract with polygons

Note that in case you're going to aggregate the numbers for the line, stars already has aggregate (similar to raster::extract with fun argument specified). If you're interested in burning geometries into a raster look at st_rasterize; this gives a few options raster::extract doesn't seem to have. Also, st_extract supports bilinear interpolation, which raster doesn't seem to have.

[kadyb]

a default value for the density needs to be chosen. What does raster do here? I noted that the number of points is not identical.

raster takes all pixels that intersect with the line, hence the difference in the number of extracted values. I think, the approximation of density = 1 / pixel size will be fine in stars.

White pixels are extracted by raster, blue points by stars

Maybe the source code will be useful in some way: raster::extractLines terra::extract

BTW: There is support for bilinear interpolation in extract in raster - argument method = "bilinear".