6

I've combed through R's raster functions and vignettes and can't seem to get this working.

I want to specify a line/vector across a raster stack (a DEM and possibly related variables), and get a profile of values for the cells which the line intersects. I've been able to do something similar using mask with a polygon.

EDIT: Thanks to scw, I have developed the following solution. 

# I have a stack of environmental rasters in this format
new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)
res(new_r) <- 0.00225
projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"
values(new_r) <- outer(seq_len(nrow(new_r)), seq_len(ncol(new_r)), "+")
stackdata <- stack(new_r, sqrt(new_r))

# I designate two transect lines by long/lat
cds1 <- rbind(c(-156, 19), c(-155.5, 20.2))
cds2 <- rbind(c(-155, 20.2), c(-155, 19.2))
transects <- SpatialLines(list(Lines(list(Line(cds1)), ID = "one"), 
                               Lines(list(Line(cds2)), ID = "two")))

# plot the lines to confirm placement
plot(new_r)
plot(transects, add = TRUE)

# and return a list whose length is equal to the number of line segments,
# and each list element is a matrix with a column for each raster layer
e <- extract(stackdata, transects)
|improve this question
  • extract function prints an annoying long string of 1s while it runs. This can be hidden with invisible(capture.output(e <- extract(...))), but is there an easier way? – J. Win. Feb 22 '11 at 20:32
8

The extract should do the trick, but you may need to update to the version of raster on CRAN first. To use it, pull in the geometries you're interested in into SpatialLines objects like so:

require raster
require rgdal

r <- raster('dem.tif')
lines <- readOGR(dsn='lines.shp', layer='lines')

elevations <- extract(r, lines)

This works well for most analysis, but isn't fast enough if you're performing very large sets of data (I have an OGR/GDAL implementation I can post somewhere if it'd be useful).

|improve this answer
  • +1. But approximately how large is "very large"? And how fast is "fast enough"? – whuber♦ Feb 22 '11 at 15:09
  • Thanks, I have done it slightly different but you put me on the right track. This takes under a minute with my sample data, so is "fast enough." My related question is at gis.stackexchange.com/questions/6424/… – J. Win. Feb 22 '11 at 20:53
  • 1
    @whuber: I had 7.1M lines (a distance matrix between geometries), and the raster extract function without optimization was doing about 8 lines/sec or 42k points/sec or 10 days for the dataset. The OGR/GDAL version computed the entire dataset in about 4.5hrs. – scw Feb 22 '11 at 22:53
  • @whuber: raster is certainly an evolving package, a coworker found an area where a small optimization to the package generated a 100x speed increase. – scw Feb 22 '11 at 23:19
  • My faster Python based version, for anyone needing to do many extraction operations: github.com/scw/topographic-distance/blob/master/lines.py – scw Mar 5 '12 at 7:34
3

If speed is an issue, consider using RSAGA with the profiles from lines module. http://www.saga-gis.org/saga_tool_doc/7.2.0/ta_profiles_4.html

|improve this answer
  • broken link, please update. – Mouad_S Feb 25 at 2:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.