Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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)
share|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. Winchester Feb 22 '11 at 20:32
add comment

2 Answers

up vote 6 down vote accepted

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).

share|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. Winchester 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
    
Thanks. That's amazingly slow. –  whuber Feb 22 '11 at 22:58
    
@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
show 1 more comment

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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