Cálculo de estadísticas de bloques a partir de un ráster de cobertura del suelo mediante R?

Question

Tengo un raster (con una resolución de 5 x 5 m) que representa la cobertura del suelo. Cada celda contiene un valor entero (de 1 a 10) que representa un total de diez clases de cobertura del suelo (por ejemplo, 1-bosque; 2-áreas agrícolas; 3-áreas urbanas; (...)10-arbustos).

A partir de este ráster de cobertura del suelo, me gustaría generar 10 rásters con una resolución de 10 x 10 m que representen la proporción de cada clase de cobertura del suelo por celda. En otras palabras, me gustaría producir un ráster que represente la proporción de bosque por celda, otro ráster que represente la proporción de áreas agrícolas por celda, etc... (véase la figura siguiente que muestra lo que quiero conseguir).

Sé que esta tarea puede llevarse a cabo en ArcGIS Spatial Analyst generando 10 rásteres binarios (1-presencia de una clase de uso del suelo específica; 0-ausencia de una clase de uso del suelo específica) y utilizando después la opción "Block statistics". Sin embargo, como no tengo licencia para ArcGIS, esta opción no es válida.

En cambio, me gustaría saber cómo hacer esto en R (versión 3.40). He revisado posts anteriores y he conseguido encontrar uno en el que un usuario pretendía hacer algo similar ( Cálculo de la proporción de clases de cobertura del suelo con una ventana móvil alrededor de un punto en R? ). Sin embargo, no quiero calcular la proporción de la cubierta terrestre con una ventana móvil.

A continuación encontrará el código para generar un objeto RasterLayer ficticio con cinco clases de cobertura del suelo codificadas como enteros.

    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

Puede utilizar el aggregate en la función raster paquete para hacerlo. Por ejemplo, para producir un ráster con la mitad de la resolución del original, mostrando la proporción de celdas cubiertas por la clase de tierra 1 :

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

Para hacer esto para todas las clases de tierra:

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

Estoy de acuerdo en que raster::aggregate() funciona bien en la mayoría de los casos, especialmente cuando se trata de una función definida por el usuario como en la respuesta dada por @dbaston. Sin embargo, es posible que desee echar un vistazo a la velox como un enfoque alternativo que resulta especialmente útil cuando se trata de rásters de gran tamaño.

La única advertencia es que VeloxRaster_aggregate() acepta funciones agregadas "estándar" (por ejemplo c("sum", "mean", "min", "max", "median") sólo - cualquier paso de trabajo más allá de ese nivel debe realizarse previamente. No obstante, la ganancia de velocidad resultante justifica sin duda la necesidad de un poco más de código. En el ejemplo siguiente, he aumentado el ncell() de su trama de entrada por un factor de 100 y los tiempos de ejecución comparados. Como se puede ver en los resultados del benchmark, velox supera al raster::aggregate() en un factor de 13, por término medio, cuando se alimenta con la imagen de trama más grande.

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)

Y aquí están los resultados de la evaluación comparativa:

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

También puede plot() o restar las capas de salida individuales en lst1 y lst2 para verificar que los dos enfoques crean resultados idénticos.

Como referencia: utilizando su trama original de 1e4 píxeles como entrada, VeloxRaster_aggregate() realiza sólo 4 veces más rápido en mi máquina (al menos en mi máquina), lo que me lleva a la conclusión de que el raster es absolutamente adecuado para las imágenes más pequeñas. Sin embargo, velox definitivamente vale la pena un vistazo más profundo cuando se trata de conjuntos de datos más grandes que todavía caben bien en la memoria.