Función Zonal Histogram de ArcGIS en R

En este tutorial explicare una aplicación de la función “zonal histogram” de ArcGIS en R, el objetivo, será calcular el área de las coberturas para los municipios del Valle del Cauca a partir de la capa de cobertura de la NASA y la división político administrativa del SIGOT.

Los datos a utilizar, en coordenadas geográficas WGS 1984, son los siguientes (disponibles aquí):

1.      Cobertura de la tierra reclasificada de la NASA.

2.      Municipios del Valle del Cauca

A continuación los pasos:

1.      Cargamos las librerías a utilizar, en caso que usted no las tenga instalada podría utilizar install.packages(“nameLibrary”)

# Cargamos librerias
require(raster)
require(rgdal)
require(foreign)
require(ggplot2)
require(dplyr)

2.      Cargamos los archivos, que son el raster, la tabla de la leyenda del raster, y el shapefile de municipios.

path <- "E:/R/_tutoriales/_tabulateArea/_blog/_data"
cob  <- raster(paste0(path, "/_raster/cov_eluglcde.tif"))
adm1 <- shapefile(paste0(path, "/_shp/mpiosValle.shp"))
leg  <- read.dbf(paste0(path, "/_raster/cov_eluglcde.tif.vat.dbf"))

3.      Ploteamos el shapefile para observar los límites de los municipios.

ggplot() + geom_polygon(data=adm1,  aes(x=long, y=lat, group=group), fill="cadetblue", color="grey") + coord_equal()

4.      Creamos la función de zonal histogram, pues esta función aún no existe como tal en alguna librería.

tabFunc                            <- function(indx, extracted, region, regname) {
  dat                              <- as.data.frame(table(extracted[[indx]]))
  dat$name                         <- region[[regname]][[indx]]
  return(dat)
}

5.      Ejecutamos la función creada anteriormente

ext <- extract(cob, adm1, method='simple')#extract valores
tabs <- lapply(seq(ext), tabFunc, ext, adm1, "NOM_MUNICI")
df <- do.call(rbind, tabs)#unimos todas las tablas en una sola

6.      Calculamos el porcentaje de cobertura por municipio y las hectáreas, teniendo en cuenta que cada píxel tiene un área aproximada de 0.20833 hectáreas

df  <-  df%>%
            group_by(name)%>%
            mutate(Porcentaje = (Freq/sum(Freq))*100)%>%
            ungroup()%>%
            mutate(HA = Freq * 0.2083333)

7.      Editamos los encabezados y unimos con la tabla de leyenda para conocer cada valor de celda a que categoría corresponde ej. El valor 9 corresponde a áreas urbanas.

colnames(df)      <- c("Category", "Count", "Municipio", "Porcentaje", "Has")
df_all            <- merge(df, leg, by.x = "Category", by.y = "Value")
df_all            <- df_all[,c("Category", "Count.x", "Municipio", "Porcentaje", "Has", "ELU_GLC_De")]
df_all            <- df_all[order(df_all$Municipio, decreasing = FALSE), ]
colnames(df_all)  <- c("Category", "Count", "Municipio", "Porcentaje", "Has", "Count2", "Cobertura")

      8.       Creamos gráfica para tres municipios del Valle del Cauca, esta gráfica representa el porcentaje de cada cobertura en el municipio según el raster de cobertura de la NASA para el año 2015.

mpios <- filter(df_all, Municipio %in% c("YUMBO", "CALI", "PALMIRA"))
gg    <- ggplot(mpios, aes(x = factor(Municipio), y = Porcentaje, group = Cobertura)) +
          geom_bar(aes(fill = Cobertura), stat = "identity", position = "dodge") +
          theme(legend.position = "bottom") + xlab("Municipio")
ggsave(plot = gg, paste0(path, "/cob_mpios_porc.png"), width = 9, height = 8, units = "in")
gg
dev.off()

9.      Guardamos la tabla como un archivo csv

write.csv(df_all, paste0(path, "/_table/zonalHistogram.csv"))

Gracias a Jeison Mesa por su colaboración en la realización de este tutorial!