17

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.

5
  • 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, 2012 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, 2012 at 13:18
  • Check r.univar in GRASS, as see grasswiki.osgeo.org/wiki/Zonal_statistics
    – markusN
    Dec 14, 2012 at 10:32
  • Hi! Thanks for making a QGIS plugin! I just wanted to mention, that the plugin crashes (>13000 polygons). It would be great if it would split up the task as to not crash. And it would be wonderful to have an option to add all classes at once (e.g. so the attribute table gets 2 new fields LandcoverID and Landcover% where both hold a CSV-list with the values) :)
    – Mfbaer
    Mar 18, 2017 at 16:33
  • @Joran : If you think this is a bug, raise a bug report rather than writing this in a comment ( github.com/Martin-Jung/LecoS/issues ). Furthermore 1) it is not the plugins job to serialize or batch process your tasks. Run it on smaller subsets then. 2) Sure. There would be many wonderful things to add. Code is open source, Feel free to code it :)
    – Curlew
    Mar 19, 2017 at 18:52

7 Answers 7

15

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.

0
6

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)
}
4

I know this post is quite old but I recently hat to carry out the same sort of analysis but downloading programs such as R is a bit of a hassle on my work computer and needs approval. After many hours of researching a method that I could use with only QGis and Excel I found this method worked the best for me and wanted to offer it to people in the same sort of situation.

  1. Clip polygon to raster layer (Raster → extraction → clipper : input file = raster layer, choose you output name & location, click on mask layer, choose you polygon → ok)

  2. Polygonise the clipper layer (Raster → Conversion → polygonise : input file = your clip layer, save output → ok)

  3. Calculating the number of pixels (Click on the shape file you’ve just created → open field calculator: tick “create new field” and add field name, Function = geometry → area → ok). You should now have a new column in your attribute table showing the number of number of pixels.

  4. Save polygonise layer (Right click polygonise layer, save as : format = DBF file, save as → ok)

  5. Summarizing number of pixels for each habitat (start excel, open file, if you don’t have titles add one now for each column, click on an empty cell, got to DATA tab, consolidate, make sure it’s on sum, click on the red arrow next to “browse…” and select you two columns (titles included), click “add” and tick both the “Top row” and “left column” boxes → ok)

  6. If, like me, you have lots of polygons to analyze and need to compare them in the same table, this step will come in helpful. In a new excel workbook list your habitats numbers in column A (for me 1 to 48) and place the two columns you’ve just consolidated in column B and C (habitat in B and number of pixels in C). In cell D1 write the following formula: =IFNA(INDEX(C:C; MATCH(A2;B:B; 0));"") and drag (or double click bottom right corner) down to your last value ( so if you have 48 habitats down to the cell D48). The number of pixels are now in the corresponding cells to you habitat and you can repeat this process for all your polygons.

3

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

2

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.

0
0

Opposed to most answers above, I would argue that the better option is to rasterize your polygons and than work with two raster data-sets instead of two polygon datasets. This is much less processing intensive and is consequently the only easy to implement solution to process large rasters and large polygon files in R.

After rasterizing your polygons to the exact same extent and resolution of the raster data, you can tabulate summary statistics as explained here, which is appropriate if your raster fits into memory (small/medium raster layers) or you can binarize each category with the reclass function and than calculate zonal statistics for each class. Here is a solution that incorporates the rasterization and zonal statistics into one function and works nice with very large datasets.

-1

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?

1
  • thanks for comment, but Zonal Statistics only eats raster without categories. Iam using QGIS Trunk 1.9
    – Curlew
    Apr 18, 2012 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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