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

|improve this question
  • 2
    Have you checked gis.stackexchange.com/questions/121116/… – underdark♦ Jan 7 '15 at 17:27
  • @underdark The table() function is helpful. I think I can already do it. Thanks. – damb Jan 7 '15 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 '15 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 – MikeRSpencer Jan 29 '15 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". – Jeffrey Evans Mar 3 '16 at 0:09
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
|improve this answer
  • 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? – user3386170 Mar 24 '18 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 – Matthew Bayly Mar 25 '18 at 2:16

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.