How to extract percentage of land cover from a raster according to multiple shapefiles?

Question

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.

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

Answer 1

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

Answer 2

I think the clearest way to do this is to use raster::extract in R, applied to files that contain binary indicators ("masks") of land cover created by raster::reclassify, for each of the land types. For example:

library(raster)

#Used only to generate random raster, not needed for your application
library(spatialEco)

#For creating sf shapefile
library(sf)

#Values used to encode land cover in the raster: 1 for grassland, 2 for forest, 3 for swamp
lcvals=c(1,2,3)

#Generate a random raster for this example
cover_raster=random.raster(n.row=10,n.col=10,n.layers=1,x=c(1,2,3),distribution="sample")
plot(cover_raster)

#Generate dummy shapefile:
shapepoints1=matrix(c(0,0,
                      0.5,0,
                      0.5,0.5,
                      0,0.5,
                      0,0),nrow=5,byrow=TRUE)
shape1=st_polygon(list(shapepoints1))
shapepoints2=matrix(c(1,0,
                      0.5,0,
                      0.5,0.5,
                      1,0.5,
                      1,0),nrow=5,byrow=TRUE)
shape2=st_polygon(list(shapepoints2))

my_shapefile=st_sf(st_sfc(list(shape1,shape2)))


#We will need a reclassification matrix to modify:
data=c(lcvals,0,0,0)
rcmat_base=matrix(data,nrow=3,ncol=2)

for (lcval in lcvals) {
  #Create a binary indicator of whether a pixel is land of the type lcval
  
  #Tell reclassify to set raster cells containing this land cover value to 1, all others to 0. This creates a binary mask for that lcval.
  rcmat=rcmat_base
  rcmat[which(rcmat[,1]==lcval),2]=1
  this_cover_mask=raster::reclassify(cover_raster, rcmat)

  #"Extract" can then calculate an average from the binary mask for that value.
  #The average of a 0/1 indicator is the fraction covered by that type of cell.
  #This will output a matrix containing the fraction with this value of cover 
  #for each shape in your shapefile.
  coverage_frac=raster::extract(this_cover_mask,my_shapefile,fun=mean,na.rm=T)

  #Append a column to your shapefile containing the fraction covered by that lcval.
  my_shapefile[,as.character(lcval)]=coverage_frac
}

The reclassify and extract operations may be parallelizable for speed using raster::clusterR.