[GIS] Calculating block statistics from land cover raster using R

land-coverland-userraster

I have a raster (5 x 5 m resolution) representing land cover. Each cell contains a integer value (from 1 to 10) representing a total of ten land cover classes (e.g. 1-forest; 2-agricultural areas; 3-urban areas;(…)10-scrublands).

From this land cover raster, I would like to generate 10 rasters with a resolution of 10 x 10 m representing the proportion of each land cover class per cell. In another words, I would like to produce a raster representing the proportion of forest per cell, another raster representing the proportion of agricultural areas per cell, etc… (please see the figure below that shows what I want to accomplish).

enter image description here

I know this task can be accomplished in ArcGIS Spatial Analyst by generating 10 binary rasters (1- presence of a specific land use class; 0-absence of a specific land use class) and then using the "Block statistics" option. However, because I do not hold a license for ArcGIS, this is not a valid option.

Instead, I would like to know how to do this in R (version 3.40). I have checked previous posts and I managed to find one where an user intended to do something similar (Calculating Proportion of Land Cover Classes with moving window around point in R?). However, I do not want to calculate proportion of land cover with a moving window.

Below find the code to generate a dummy RasterLayer object with five land cover classes coded as integers.

    r<-raster(ncols=100, nrows=100,xmn=100,xmx=1100,ymn=0,ymx=1000, crs="+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")

    x<-0:5

    r[]<-sample(x,10000,replace=T)

Best Answer

You can use the aggregate function in the raster package to do this. For example, to produce a raster at half the resolution of the original, showing the proportion of cells covered with land class 1: 

aggregate(r, fact=2, fun=function(vals, na.rm) {
  sum(vals==1, na.rm=na.rm)/length(vals)
})

To do this for all land classes:

cov_pct <- lapply(unique(r), function(land_class) {
             aggregate(r, fact=2, fun=function(vals, na.rm) {
               sum(vals==land_class, na.rm=na.rm)/length(vals)
             })
           })
Related Question