Extract raster layer data by attributes in RasterStacks R

Question

I have a RasterStack with 2 raster layers of the same area, one from a SpatialPolygonsDataFrame (layer1) and another NDVI satellite image (NDVI). layer1 has a series of polygons, each with various land-use attributes.

I am interested in using attributes in layer1 (those partitioning the layer according to various Factors), to query the NDVI layer. I'd like to be able to compare different levels of the factor statistically, or graphically by boxplots or similar. I'd also like to be able to add additional NDVI images, and compare across NDVI layers for the same Factors.

Is it possible to do this - which functions / packages should I look at? extract() allows queries by separate Spatial objects, doesn't seem to allow query by attributes. I haven't been able to find a solution on this site.

Clearly these operations can be done by using separate Raster and Spatial objects, but I suspect using a RasterStack would be more efficient. Or is there another way altogether?

Answer 1

Your question is a good example of why site guidelines request that you provide an example of what you have already tried and ideally, some example data. There are many aspects of your question that leave one speculating. For example, I am not clear as to what you mean by "one from a SpatialPolygonsDataFrame". Is your land-use raster data rasterized from polygons?

If this is the case, then the answer is quite simple. Leave the land-use data as a SpatialPolygonsDataFrame, query it to create a desired subset and then use extract. This is clearly illustrated in the help for the extract function.

In this example we create some dummy polygon data with two classes (ag, urban) and a raster bounded -1 to 1 (like NDVI). We then query the polygons to extract values associated with a single class.

library(raster)
r <- raster(ncol=36, nrow=18)
  r[] <- runif(ncell(r), -1,1)
lu <- spPolygons( rbind(c(-180,-20), c(-160,5), c(-60, 0), 
                  c(-160,-60), c(-180,-20)),rbind(c(80,0), 
                  c(100,60), c(120,0), c(120,-55), c(80,0)))
lu <- SpatialPolygonsDataFrame(lu, data.frame(class=c("urban","ag")))

# Extract raster values associated with "urban"
extract(r, lu[lu$class == "urban",])

# Here we can pull a summary
summary(unlist(extract(r, lu[lu$class == "urban",])))

# Boxplot of multiple classes (this could be dangerous with large data)
boxplot( unlist(extract(r, lu[lu$class == "urban",])), 
         unlist(extract(r, lu[lu$class == "ag",])), pch=20)

If you have a raster stack (eg. representing multi-date NDVI) then the data structure, resulting from extract, will be different. Instead of a list with a vector of values for each polygon there will be a list containing a data.frame for each polygon where each column represents a raster in the stack. You can use do.call to create a single data.frame.

Answer 2

Perhaps zonal?

library(raster)
r <- raster(ncols=10, nrows=10)
r[] <- runif(ncell(r)) * 1:ncell(r)
z <- setValues(r, rep(1:5, each=20))
zonal(r, z, 'sum')

or

boxplot(r, z)

which is equivalent to

x <- cbind(v=values(r), zone=values(z))
boxplot(v~zone, data=x)