4

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

|improve this question
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)
|improve this answer
  • Thanks. That is pretty much what I got from gis.stackexchange.com/questions/23614/…. 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

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.