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?

  • You should provide some sample data in order to get proper answers. – fdetsch Feb 18 '16 at 7:35

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.

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.

  • Thanks Jeffrey. I just returned to post some example code, and here is your suggestion. My examples include solutions similar to those you have outlined; using factors from SpatialPolygonDataFrames to query rasters. To answer your question, I was hoping however to use factor data (land use) from a rasterized SPDF layer within a RasterStack to query NDVI layers, also in the stack, as potentially be a neater way to do things - especially if there were a specific function for the purpose. Your method allows me to perform the stats I was having trouble with however, so thanks again. – Daniel Feb 19 '16 at 4:08

Perhaps zonal?

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


boxplot(r, z)

which is equivalent to

x <- cbind(v=values(r), zone=values(z))
boxplot(v~zone, data=x)
  • Cool RobertH, thanks a bunch – zonal is getting close. Now, if I could get it to treat the levels of factor in a stack attribute as zones, I’d be really happy. Here, using 3 levels in ‘Cult.Hist’ attribute of WandStack as zones in which to sum values from raster aug14. Maybe it’s only failing because of my poor syntax? zonal(aug14, WandStack@layers[[1]]@data@attributes[[1]]$Cult.Hist, 'sum', na.rm=TRUE ) – Daniel Feb 21 '16 at 3:05
  • Yes, your syntax looks confused. Reading slots (@) directly is rarely a good idea. Please provide a reproducible example. – Robert Hijmans Feb 21 '16 at 5:37

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.