[GIS] Step-by-step: How to extract Raster values from Polygon overlay with Q-GIS or R

polygonrraster

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!

Best Answer

One problem I see is that if you do not have all classes present in each polygon then the resulting data does not match the dimension of the data frame in the polygon @data slot. You can avoid a looping function entirely. Here is a quick fix that will add NA values where you do not have a class present. 

require(raster)

# Create interger class raster
r <- raster(ncol=36, nrow=18)
  r[] <- round(runif(ncell(r),1,10),digits=0)

# Create two polygons
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                  Polygons(list(Polygon(cds2)), 2))),data.frame(ID=c(1,2)))

# Extract raster values to polygons                             
( v <- extract(r, polys) )

# Get class counts for each polygon
v.counts <- lapply(v,table)

# Calculate class percentages for each polygon
( v.pct <- lapply(v.counts, FUN=function(x){ x / sum(x) } ) )

# Create a data.frame where missing classes are NA
class.df <- as.data.frame(t(sapply(v.pct,'[',1:length(unique(r)))))  

# Replace NA's with 0 and add names
class.df[is.na(class.df)] <- 0   
  names(class.df) <- paste("class", names(class.df),sep="")

# Add back to polygon data
polys@data <- data.frame(polys@data, class.df)
  head(polys@data)
Related Question