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