[GIS] How to extract percentage of land cover from a raster according to multiple shapefiles

grassqgisrrastershapefile

I have a raster map that represents land cover of a region, and different shapefiles (total 120) in the same region. I want to know the percentage of every land cover for each shapefile.enter image description here

For example, I want to know % of forest, % of cultivated area, % of grassland, and so on, for each shapefile.

I already tried clipping the raster with the shapefiles in Qgis, but is a slow process. I think with R is possible too, but I have not been able to do it yet:

setwd("/home/shapefiles")
shps <- dir(getwd(), "*.shp") #looking for shapefiles

library(maptools)
list.shape<-list()
for (shp in shps){
  list.shape[[shp]]<- readShapeSpatial(shp) #Reading shapefiles
} 

#Opening raster
library(raster); library(rgeos); library(rgdal)
r.landcover <- raster('/home/raster')
shape.integrated<-c()
for (k in 1:length(list.shape)){
  shape.integrated<-gUnaryUnion(list.shape[[k]])  # Joining shapefiles
  r.vals[[k]]<- extract(r.landcover, shape.integrated)  # Extract values of the raster according to the shapefile 
  percent[[k]]<- lapply(r.vals[[k]], FUN=I NEED AN OPTION THAT CALCULATE PERCENTAGE OR SOMETHING SIMILAR HERE)  # Percentage
}

I can use Qgis, R, GRASS (Open source).

Best Answer

One solution would be to convert each of your polygons to rasters. Then take the portion of cell counts in each landcover class by the total cell count from the clip. Here is an example function: 

    PortionClassInPoly <- function(MySingleShape, MainRaster) {

    # Start with one shape from your list
    # Rasterize this shape according to the parameters of your background raster
       mini_rast <- rasterize(MySingleShape, MainRaster, fun='last')  
    # divide by itself - to make new set new raster values to 1 (NA values stays as NA). 
       mini_rast <- mini_rast/mini_rast
    # get the total number of cells for this shape (save for below)
       total_cells <- cellStats(mini_rast, 'sum')

    # multiply by original - to make a mask layer from your original raster
       my_cutout <- MainRaster*mini_rast
    # Count the number of cells for each discrete landcover class
       in_this_poly_unit <- freq(my_cutout)
    # Divide the cell count for each of these classes by the total 
    # number of cells in your current shape
       class_percents <- in_this_poly_unit[,2]/total_cells

    # make the function return a dataframe with
    # % landcover in each class for each shape
       output <- data.frame(myclass=in_this_poly_unit[,1], portion=class_percents)

       return(output)
    }

With your data, copy and paste the above & then try this function in your loop:

    PortionClassInPoly(MySingleShape=list.shape[[k]], MainRaster=r.landcover)

    # you should get something like this:
    # percentage of every land cover for a shapefile 

      myclass   portion
    1       1 0.1726333
    2       2 0.1662689
    3       3 0.1718377
    4       4 0.1789976
    5       5 0.1662689
    6       6 0.1439936
Related Question