Saturday, 16 December 2017

Step-by-step: How do I extract Raster values from Polygon overlay with Q-GIS or R?


I am writing my master thesis at the moment and am having the exact same problem as user Curley in this question: I first need to overlay a Corine 2006 raster with a polygon layer (2638 polygons). Then I need to extract the proportion (%) for each class within a polygon layer so I get a table like


                1 (e.g. arable land)           2 (e.g. forest)          3 (e.g. water)

Polygon 1 11 % 0% 89%
Polygon 2 23 % 60% 17%

The statistics plugin Lecos that is being provided crashes my computer despite updating the libraries.


Since I am lacking experience am not yet able to completely reproduce the solution in R. I could follow the example as far as


> 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 ...

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

26 27
287 255

[[2]]


2 11 18
67 99 792

but when I try to reproduce Curlew's next step:


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

I'll receive the following Error message:


 [Error] length of 'dimnames' [2] not equal to array extent

So I created the list of tables with lapply. Now I need to attach it to the attribute table of my polygon layer.


The first row of each table in this list contains the landscape classes (their ID) that lie within the polygon. The second row is the number of raster cells of each class within the polygon. It is presence-only, so the tables never show a class with zero cells and have different sizes.


The landscape class ID should now be the header of the new columns I want to attach to the polygon layer.



Then I need to calculate the percentage of the classes for each polygon.


I would be very grateful if you could give me a more detailed step-by-step description of what happened in the final step by Curlew, or maybe provide another solution in R as I will be using it for my data analysis later. Thanks!




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...