3

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

5
  • 2
    Have you checked gis.stackexchange.com/questions/121116/…
    – underdark
    Jan 7, 2015 at 17:27
  • @underdark The table() function is helpful. I think I can already do it. Thanks.
    – D.MerBet
    Jan 7, 2015 at 19:54
  • and don't use readShapeSpatial because it will ignore the coordinate reference system. Use readOGR from the rgdal package.
    – Spacedman
    Jan 10, 2015 at 22:41
  • The r.stats tool in GRASS would do this. You'd need to convert your polygons to raster to mask off the areas you wanted to run it on. If you're working with rasters I wrote a how to here, writing the results to a database: scottishsnow.wordpress.com/2014/08/24/many-rastered-beast Jan 29, 2015 at 20:23
  • Do each of your 120 shapefiles contain multiple polygon features or just one? The way you code this problem is dependent on this answer. One point of confusion is that it is sounding like you want a global measure of landcover proportions for each "shapefile". Mar 3, 2016 at 0:09

2 Answers 2

3

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
2
  • I don't understand the list.shape[[k]] syntax. I have a shapefile loaded with multiple polygons, would an equivalent be file[1,] or file[k,] within a loop? Mar 24, 2018 at 0:42
  • 1
    Yeah 'k' would just be an index in the loop, assuming polygons are provided as a list. Depending on how your object is structured you should be able to use something like MyPolygonLayer[k,] instead Mar 25, 2018 at 2:16
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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