Problem with creating elevation cross sections using R and terra

[adamkc]

The terra::extract function samples raster data from the top to the bottom and from left to right. This creates an issue if you're using a line to sample a raster with plan of plotting the data as a cross section. The returned data might be expected to be sequenced in the order the pixels are encountered by the line but that isn't the case.

Below is a reproducible example of the problem. I came up with one solution (coded below) but I suspect this is a common nuisance so maybe there is another solution to consider?

# Load packages
library(terra)

# Create a sample raster  
r <- rast(nrow=10, ncol=10)
r[] <- 1:ncell(r)

# Create a line geometry
l <- sf::st_linestring(matrix(c(-150, 150,150, -20, 20,-50), ncol=2))
l <- vect(l)

# Extract raster values along line
extract <- terra::extract(r, l,xy=TRUE)
extract$sample_order <- 1:nrow(extract)

# Plot the results
plot(r)
plot(l,add=TRUE)
text(extract$x, 

extract$y, labels=extract$sample_order, cex=0.8,add=TRUE)

My solution:

l2 <- terra::densify(l,interval =min(res(r)), flat=TRUE)
p <- as.points(l2)
extract2 <- terra::extract(r, p,xy=TRUE)

[adamkc]

A second workaround was posted here: StackExchange Post.

This appears to be a common use case of extract, and a commonly overlooked issue leading to frustrating "sometimes the plots look wrong" posts (like here: StackExchange post with incorrect answers). It would be great if terra had a native solution, such as terra::extract(r,l, sample_line_in_order=TRUE) or terra::extract_line_in_order(...).

[rhijmans]

I have added a new function extractAlong(x, y). It can probably be optimized a bit more for speed.

library(terra)
r <- rast(nrow=10, ncol=10, vals=1:100)
v <- vect(matrix(c(-150, 150,150, -20, 20,-50), ncol=2), type="lines")

extractAlong(r, v)
#   ID cell lyr.1
#1   1  338   338
#2   1  302   302
#3   1  266   266
#4   1  267   267
#5   1  231   231
#6   1  232   232
#7   1  196   196
#8   1  197   197
#9   1  161   161
#10  1  162   162
#11  1  126   126
#12  1  127   127
#13  1  163   163
#14  1  164   164
#15  1  200   200
#16  1  201   201
#17  1  237   237
#18  1  273   273
#19  1  274   274
#20  1  310   310
#21  1  310   310
#22  1  346   346
#23  1  382   382
#24  1  381   381
#25  1  417   417
#26  1  453   453
#27  1  452   452
#28  1  488   488
#29  1  488   488
#30  1  487   487
#31  1  451   451
#32  1  450   450
#33  1  414   414
#34  2  279   279
#35  2  243   243
#36  2  244   244
#37  2  208   208
#38  2  209   209
#39  2  174   174
#40  2  175   175
#41  2  139   139
#42  2  140   140
#43  2  141   141
#44  2  177   177
#45  2  213   213
#46  2  214   214
#47  2  250   250
#48  2  286   286
#49  2  322   322
#50  2  358   358
#51  2  394   394
#52  2  430   430
#53  2  429   429
#54  2  465   465
#55  2  501   501
#56  2  537   537

[Nowosad]

@rhijmans I believe the new function still can be improved: (a, first example) adding sf objects support (that is already available in extract) and (b, second example) the new function does not work for my data example.

library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.7.42

srtm = rast(system.file("raster/srtm.tif", package = "spDataLarge"))

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) |>
  st_linestring() |> 
  st_sfc(crs = crs(srtm)) |>
  st_sf(geometry = _)

zion_transect_vect = vect(zion_transect)

plot(srtm)
plot(zion_transect, add = TRUE)

exline1 = extractAlong(srtm, zion_transect)
#> Error in extractAlong(srtm, zion_transect): inherits(y, "SpatVector") is not TRUE
exline1
#> Error in eval(expr, envir, enclos): object 'exline1' not found

exline2 = extractAlong(srtm, zion_transect_vect)
#> Warning in y@pnt$rasterize(x@pnt, field, values, background, touches[1], : GDAL
#> Error 1: Failed to compute min/max, no valid pixels found in sampling.
#> Error in data.frame(i, cell, extract(x, cell)): arguments imply differing number of rows: 1, 0
exline2
#> Error in eval(expr, envir, enclos): object 'exline2' not found

[rhijmans]

Thanks for checking. I now get:

exline1 = extractAlong(srtm, zion_transect)
head(exline1)
#  ID  cell srtm
#1  1 34923 2001
#2  1 34924 1986
#3  1 35389 2033
#4  1 35390 2025
#5  1 35855 2037
#6  1 35856 2020

[Nowosad]

@rhijmans one more question: when xy=TRUE, then the X and Y coordinates of the raster cells are returned. Would not it be better to return the coordinates of the extracted points?

The current situation (when zoomed in) is:

[rhijmans]

The line touches cells, and the coordinates for these cells are returned. I think these are the coordinates for the extraction.

But it would be possible to also/alternatively compute, for each cell, the coordinates on the line that are nearest to the cell center.

[Nowosad]

Thanks @rhijmans. To explain myself: the coordinates on the line would be useful to create transects (e.g., see https://r.geocompx.org/raster-vector#fig:lineextr).

[rhijmans]

Here is an implementation. As-is, the coordinates will only be correct for one case, for cells where the line crosses multiple times. That is not expected for a transect. But it can be fixed.

library(terra)
r <- rast(nrow=10, ncol=10, vals=1:100, crs="+proj=utm +zone=1")
v <- vect(matrix(c(-150, 150,150, -20, 20,-50), ncol=2), type="lines", crs=crs(r))
e1 <- extractAlong(r, v, xy=TRUE)
e2 <- extractAlong(r, v, xy=TRUE, online=TRUE, bilinear=FALSE)

plot(as.polygons(r), border="gray", axes=FALSE); lines(v, lwd=2)
points(e1[,c("x", "y")], cex=1, col="blue"); points(e2[,c("x", "y")], col="red", pch=20) 

Note the default bilinear=TRUE when using snapped coordinates. This may give an error or warning with categorical rasters. I have not tried that yet.

e3 <- extractAlong(r, v, xy=TRUE, online=TRUE)
head(e1)
#  ID    x   y lyr.1
#1  1 -162 -27    61
#2  1 -126 -27    62
#3  1 -126  -9    52
#4  1  -90  -9    53
#5  1  -54  -9    54
#6  1  -18  -9    55

head(e2)
#  ID          x          y lyr.1
#1  1 -150.00000 -20.000000    61
#2  1 -127.33624 -16.978166    52
#3  1 -124.97817 -16.663755    52
#4  1  -89.60699 -11.947598    53
#5  1  -54.23581  -7.231441    54
#6  1  -18.86463  -2.515284    55

head(e3) # bilinear
#  ID          x          y    lyr.1
#1  1 -150.00000 -20.000000 57.44444
#2  1 -127.33624 -16.978166 56.39520
#3  1 -124.97817 -16.663755 56.28603
#4  1  -89.60699 -11.947598 54.64847
#5  1  -54.23581  -7.231441 53.01092
#6  1  -18.86463  -2.515284 51.37336

[rhijmans]

This approach places the point where the line is closest to the center of the cell. Another approach would be to place in the middle of the line segment inside a cell. That would noticeably change the location for the points in the 2nd, 6th and 9th column of the raster shown above.