Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have two layers. A polygon-shape-layer with many tiles and a raster-layer containing CORINE 2006 land cover with many categories in a colourmap. I want to obtain for each polygon in the shapelayer a sum of each land cover category of the raster-layer.

For example there is a polygon with id '2' and i want to Attributes like this for this polygon (in percent or square meters):

  • Arable land : 15 %
  • Forest: 11 %
  • Streets:2 % (... and so one)

I tried to do it in grass, qgis (no function), saga (just sums up every to a total value) r(total sum), but i still found no solution. Most plugins (zonal statistics in qgis) only support 0-1 raster layers. v.rast.stats didn't help either. Iam open to any good and smart solution!. Maybe i even used a wrong approach or made mistakes.

In Arcgis this task is quite easy, if am remember right, but i am still missing a good solution for your everyday linux user.

I am running a debian linux system and this why i can only use programs for this OS.


EDIT: Because this question still has so many views and visitors: I wrote a QGIS-plugin, which also is capable of calculating the landcover of raster layer. I have'nt coded a polygon overlay yet, but it definitely planed. Find the plugin here and install the Scipy library first.

share|improve this question
    
It can definitely be done in R, its just a question of working out which functions. You need to overlay each polygon with the raster, and then use table() to get a summary of the "cookie-cut" pixels. Packages raster, rgdal, and rgeos may be useful. Read the "R Spatial Task View" (google will find it) –  Spacedman Apr 17 '12 at 12:49
    
sure, but how can i get such a summary. You can easily overlay a polygon layer with a raster layer with !is.na(overlay(Poly, Raster)), but with commands like extract i can only calculate the total area in the cookie-cut pixel and not different categories of a colourmap. I didn't know rgeos, but i looked through the api and found no function to do this. –  Curlew Apr 17 '12 at 13:18
    
Check r.univar in GRASS, as see grasswiki.osgeo.org/wiki/Zonal_statistics –  markusN Dec 14 '12 at 10:32
add comment

5 Answers

up vote 8 down vote accepted

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.

share|improve this answer
    
+1 ... beat me by a few minutes! ;) –  Simbamangu Apr 17 '12 at 14:38
    
Great, thanks a lot for the effort. I will try that and report back :-) –  Curlew Apr 17 '12 at 16:21
add comment

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)
}
share|improve this answer
add comment

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

share|improve this answer
add comment

How about converting the CORINE data into a vector polygon dataset using QGIS (Raster > Conversion > Polygonize) and then using the Union function (Vector > Geoprocessing Tools > Union) to combine with the polygons. The resulting vector dataset would contain the areas of each CORINE class in each polygon.

share|improve this answer
    
thanks for this suggestion. Haven't thought about vector union yet. Maybe i will try that, if R-processing somehow fails. –  Curlew Apr 17 '12 at 16:23
add comment

QGIS.

In the QGIS trunk, there is another version of ZonalStats available, it is called Zonal Statistics.

This carries out the function you require.

As to the workflow, I am not clear as to how many rasters you have or are they just bands in a raster?

share|improve this answer
    
thanks for comment, but Zonal Statistics only eats raster without categories. Iam using QGIS Trunk 1.9 –  Curlew Apr 18 '12 at 14:43
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.