Get Raster Values from a Polygon Overlay in Opensource GIS Solutions

Use 'extract' to overlay polygon features from a SpatialPolygonsDataFrame (which can be read from a shapefile using maptools:readShapeSpatial) on a raster, then use 'table' to summarise. Example:

> c=raster("cumbria.tif") # this is my CORINE land use raster
> summary(spd)
Object of class SpatialPolygonsDataFrame
[...]
> nrow(spd)  # how many polygons:
[1] 2
> ovR = extract(c,spd)
> str(ovR)
List of 2
 $ : num [1:542] 26 26 26 26 26 26 26 26 26 26 ...
 $ : num [1:958] 18 18 18 18 18 18 18 18 18 18 ...

So my first polygon covers 542 pixels, and my second covers 958. I can summarise each of them:

> lapply(ovR,table)
[[1]]

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792 

So my first polygon is 287 pixels of class 26, and 255 pixels of class 27. Easy enough to sum and divide and multiply by 100 to get percentages.


I wanted to report back and here i am. Spacedman's solution worked great and i was able to export all information for every polygon in my shape. Just in case someone has the same problem, here is how i preceded:

...
tab <- apply(ovR,table)
# Calculate percentage of landcover types for each polygon-field.
# landcover is a datastream with the names of every polygon
for(i in 1:length(tab)){
 s <- sum(tab[[i]])
 mat <- as.matrix(tab[[i]])
 landcover[i,paste("X",row.names(mat),sep="")] <- as.numeric(tab[[i]]/s)
}

if I understand correctly what you want, and assuming you have the vector layer 'mypolygonlayer' and the raster layer 'corina' already in your GRASS GIS database:

First I would convert vector to raster. The cat will ensure you'll have one unique identifier per polygon. If you have a column with a unique numerical identifier, you can use that column instead. The labelcolumn is optional:

v.to.rast input=mypolygonlayer layer=1 output=mypolygons use=cat labelcolumn=NameMappingUnit

Then run r.stats to get your statistics:

r.stats -a -l input=mypolygons,corina separator=; output=/home/paulo/corinastats.csv

The last step is to open the corinastats.csv in e.g., LibreOffice and create pivot table or use R to create your cross table