Calculating block statistics from land cover raster using R?

Question

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

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)

Answer 1

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

Answer 2

I agree that raster::aggregate() works well in most cases, especially when a user-defined function is involved as in the answer given by @dbaston. However, you might want to have a look at the velox package as an alternative approach that comes in particularly handy when dealing with larger rasters.

The only caveat here is that VeloxRaster_aggregate() accepts "standard" aggregate functions (ie. c("sum", "mean", "min", "max", "median") only – any work steps beyond that level need to be carried out beforehand. Nevertheless, the resulting speed gain surely justifies the need for a little more code. In the example below, I increased the ncell() of your input raster by a factor of 100 and compared runtimes. As seen from the benchmark results, velox outperforms the raster::aggregate() based approach by a factor of 13, on average, when fed with the larger raster image.

library(raster)

r = raster(ncols = 1e3, nrows = 1e3
           , 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
set.seed(123) # make sampling reproducible
r[] <- sample(x, 1e6, replace = TRUE)


### benchmark ----

library(microbenchmark)
library(velox)

microbenchmark(

  ## raster approach
  rasterAgg = {
    lst1 = 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)
      })
    })
  }

  ## velox approach
  , veloxAgg = {
    vls = r[]

    lst2 = lapply(unique(r), function(i) {
      # set the current land cover class to 1, everything else to 0
      ids = (vls == i)
      vls[ids] = 1; vls[!ids] = 0
      rst = setValues(r, vls)

      # create velox raster and aggregate
      vlx = velox(rst)
      vlx$aggregate(factor = 2, aggtype = "mean")
      vlx$as.RasterLayer()
    })
  }

  ## repeat every chunk n times
  , times = 10)

And here are the benchmark results:

Unit: milliseconds
expr       min        lq      mean    median        uq       max neval cld
rasterAgg 5844.2713 5950.7560 6140.8879 6186.1065 6354.8678 6370.6968    10   b
veloxAgg   319.2344  336.5698  472.3171  375.2693  622.2361  744.8573    10  a

You can also plot() or subtract the single output layers in lst1 and lst2 from each other in order to verify that the two approaches create identical results.

As a reference: using your original 1e4 pixels raster as input, VeloxRaster_aggregate() performs only 4 times faster on my machine (at least on my machine), which leads me to the conclusion that the raster based approach is absolutely suitable for smaller images. However, velox is definitely worth a deeper look when dealing with larger datasets that still fit well into memory.