combine raster and polygon values in R

Question

My question is similar to this one Get Raster Values from a Polygon Overlay in Opensource GIS Solutions but I think I need another step.

I have a polygon layer of ecoregions which I brought into R with "readOGR". It has many attributes such as ecoregion and biome (biomes encompass multiple ecoregion polygons). I have a raster layer (I used "raster" on a tif) that is 0-2. I also have continuous raster layers, that I'll use later so I'm looking for a generalized solution.

I would like to be able to do a variety of evaluations like total area and proportion of each biome >0 or when I use the continuous raster, things like mean or range. Ideally I'd like to end up with some sort of attribute table that I can work with that combines the polygon attributes with the raster values. I'm not sure if it's most accurate & efficient to convert both to polygons, both to rasters, or to do summaries as is.

I am trying to work in R so that I can script the process; sample code is particularly helpful as I'm a novice. I appreciate any thoughts.
Thanks a lot

Answer 1

The following R example essentially performs what ArcGIS users call Zonal Statistics. This should be a good building block for your analysis. The main function performing this analysis is extract() from the raster package.

require(raster)

# Create some sample raster data
r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
plot(r)

#Create some sample polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                              Polygons(list(Polygon(cds2)), 2)))
plot(polys)

# Extract the raster values underlying the polygons
v <- extract(r, polys)
v

# simplify to display mean values
output = unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA ))
print(output)

Edit:

Here is a simple calculation of the zonal mean using your own shapefile and single band raster:

Which results in the mean pixel value for each polygon.

require(raster)
require(maptools)

# Read the polygon shapefile
poly = readShapePoly("C:/temp/poly.shp")
plot(poly)

# Read the single band raster
raster = raster("C:/temp/subset.tif")

# Extract the raster values underlying the polygons
v <- extract(raster, poly, fun = mean)
output = data.frame(v)
print(output)