5

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

Community
  • 1
user20353
  • 51
  • 1
  • 1
  • 2

1 Answers1

6

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:

enter image description here

Which results in the mean pixel value for each polygon.

enter image description here

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)
Aaron
  • 51,658
  • 28
  • 154
  • 317
  • Thanks. That is pretty much what I got from http://gis.stackexchange.com/questions/23614/get-raster-values-from-a-polygon-overlay-in-opensource-gis-solutions. I guess the information I want is probably in there and I just have to figure out how to get it out. Does it exist as a table with the polygon attributes and the raster values? Then I could summarize or use a formula to pull out information by different attributes. Thanks. – user20353 Jul 23 '13 at 18:44
  • 3
    In this case the object resulting from extract is a list where each element in the list represents a polygon and the associated underlying pixel values. This is why the end if the code uses lapply (list apply). @Aaron just got fancy with a custom mean function that ignores empty polygon sets. Although, this is not necessary because lapply will recycle NA's. This does, however illustrate how you can pass a custom function to lapply. Nesting lapply in unlist returns an ordered vector that can be joined back to the data slot in the polygon data. – Jeffrey Evans Jul 23 '13 at 18:51
  • Okay, I think that last sentence is my key. So if I can figure out how to implement it, I can get the table I am imagining! – user20353 Jul 23 '13 at 19:08