[R-sig-Geo] Overlay of SpatialGrid/Raster and polygons

[Robert J. Hijmans]

Hi Ned,

Here is an approach to get values from a RasterStack to all cells in
each polygon:

library(raster)
# a polygon
data(meuse.riv)
pol <- SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)), "x")))
# a raster
r1 <- raster(system.file("external/test.ag", package="sp"))
r2 <- sqrt(r1)
# a stack
s <- stack(r1, r2)

# convert the polygon to a RasterLayer
rr <- polygonsToRaster(pol, s)

par(mfrow=c(1,2))
plot(r1)
plot(pol, add=TRUE)
plot(rr)

# extract points that are not NA
pts <- rasterToPoints(rr) # see additional arguments to select a
subset (useful for very large rasters)
# perhaps subsample your points
sampl <- sample(1:length(pts[,1]), min(100, length(pts[,1])))
pts <- pts[sampl,]

v <- xyValues(s, pts[,1:2])
# if you have mutiple polygons, bind the polygon ID to the raster values:
v <- cbind(pts[,3], v)

# You could also sample points with spsample
pts <- spsample(pol, 100,  "random")
# remove duplicate cells
cells <- unique(cellFromXY(s, pts))
v2 <- cellValues(s, cells)

# or sample from a SpGDF if you can create it from the RasterLayer (if
it is not too big)
spdf <- as(rr, 'SpatialGridDataFrame')
pts <- spsample(spdf, 100,  "random")
v3 <- xyValues(s, pts)


For extremely large rasters, rasterToPoints could fail. If so, first
create a RasterLayer with random values with calc and "fun=runif" and
make values F for e.g x > 0.01, overlay that with rr, and try
rasterToPoints again...

Hth, Robert




On Mon, Nov 23, 2009 at 12:16 PM, Ned Horning <horning at amnh.org> wrote:
> Greetings,
>
> I am looking for a way to use R to extract pixel values (currently in a
> large RasterStack object) that fall under polygons (currently a
> SpatialPolygonsDataFrame object). I seem to recall a discussion about using
> overlay to do this but I can't find a method that would work.
>
> Any insight would be appreciated.
>
> Ned
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>