1 Introducción.

Los incendios forestales causan daños económicos, afectan el suelo, alteran los ciclos hidrológicos, provocan serios deterioros a los ecosistemas forestales y a la biodiversidad y aportan carbono a la atmósfera, contribuyendo en el calentamiento global de la tierra. Aunque en condiciones naturales los incendios forman parte del proceso dinámico de los ecosistemas (Ferreras, J. et al, 2001).

En las últimas décadas, los incendios forestales han sido una preocupación en distintas regiones del mundo, especialmente en Chile, por el incremento en su ocurrencia producto de actividades humanas y cambio climático. En este trabajo se presenta una aplicación práctica de Análisis Multicriterio desarrollada en el software R, utilizando datos vectoriales y ráster obtenidos de la base de datos del Centro de Información de Recursos Naturales (CIREN), IDE-Chile, e Imágenes del satélite Landsat-8.

Para plantear este trabajo se hizo una breve revisión de los elementos y factores geográficos que configuran la zona de estudio, lo cual permitió establecer la base de trabajo para el análisis espaciales, posteriormente se explican los conceptos teóricos y metodológicos desarrollados en R, y por último se analizan los resultados.

2 Área de estudio.

El área de estudio está ubicada en la región del bío bío, a unos 90 kilómetros al Sur-Este de la ciudad de Concepción, y comprende una superficie aproximada de 6.750 km2. Desde el punto de vista morfológico, se distinguen las unidades tradicionales del relieve chileno como lo son la cordillera de los Andes, depresión intermedia, cordillera de la Costa y planicies litorales.

La región marca la transición entre los climas templados secos de la zona central de Chile y los climas templados lluviosos que se desarrollan inmediatamente al sur del río Bío Bío. En la franja costera y en los sectores altos y laderas occidental de la cordillera de la Costa se presenta un clima templado húmedo, con una humedad constante con precipitaciones que fluctúan entre 1.200 mm y 2.000 mm anuales de norte a sur de la región. Hacia el interior el clima templado costero húmedo posee también temperaturas menos extremas donde las precipitaciones alcanzan 1.330 mm anuales con un período seco de cuatro meses. En el valle longitudinal las temperaturas presentan un mayor contraste entre día y noche. En la parte norte de la región predomina el clima templado mediterráneo abarcando toda la zona intermedia, bordes orientales de la cordillera de la Costa y los sectores más bajos de la precordillera. En la cordillera de los Andes por sobre los 1.500 metros de altura se desarrolló el clima frío de altura con abundantes precipitaciones, más de 2.000 mm anuales y las bajas temperaturas que permiten la presencia de nieves permanentes en las alturas de la cordillera. Las características climáticas que presenta la región en su extremo norte permiten la existencia del espino, asociado con boldo, peumo y quillay. En cambio hacia el sur se encuentra el bosque esclerófilo, en donde hoy es posible apreciar el cambio que se ha producido en la vegetación natural por las plantaciones forestales de pinos y por cultivos agrícolas.

Al sur del Bío Bío se ubica el bosque templado higromórfico, principalmente en la cordillera de la Costa y en la precordillera andina, donde predominan especies como roble, ciprés, coigüe, lenga y ñirre y en los sectores con mayores alturas es posible encontrar alerce y mañío. Además está acompañado por un denso sotobosque formado por canelo, olivillo, avellano y especies menores como el copihue, quila y ulmo.

Mapa del área de estudio.

Mapa del área de estudio.

3 Objetivo del análisis.

El objetivo fue obtener un modelo de áreas susceptibles a los incendios forestales en la región del Bio Bio. Este estudio en esta zona es importante ya que la industria forestal chilena dispone de más de un millón de hectáreas plantadas de eucaliptus y pino radiata, lo que sustenta aserraderos, fábricas de paneles, enchapados e industrias de celulosa y papel.

4 Materiales y métodos.

Para el desarrollo de este trabajo se utilizó el software R, y la librería GDAL, el entorno Google Earth Engine y el software QGis. A continuación señalamos los paquetes y herramientas usados.

  1. Paquetes R:
    • base
    • rmarkdown
    • sp
    • rgdal
    • raster
    • sf
    • landscapemetrics
  2. Herramientas GDAL:
    • gdalbuildvrt
    • gdalinfo
    • gdal_translate
    • gdalwarp
    • gdal_calc.py
    • val_repl.py
    • ogr2ogr
    • ogrinfo

Los datos vectoriales fueron obtenidos de las bases de datos del Centro de Información de Recursos Naturales (CIREN) https://www.ciren.cl/ , IDE-Chile http://www.ide.cl/ y las Imágenes del satélite Landsat-8 se obtuvieron utilizando Google Earth Engine https://code.earthengine.google.com/

Lista de capas de información

  1. Uso del suelo.
  2. Bandas Landsat:
    • Banda B2
    • Banda B3
    • Banda B4
    • Banda B5
    • Banda B6
    • Banda B7
  3. Imagen raster NDVI.
  4. Imagen raster DEM.
  5. Imagen raster Aspect.
  6. Imagen raster Slope.

Instalamos paquetes necesarios para ejecutar la sesión en código R.

if(!require(rgdal))  suppressMessages( suppressWarnings( install.packages("rgdal") ))
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
if(!require(sf)) suppressMessages( suppressWarnings( install.packages("sf") ))

Definimos una estructura de directorio con una carpeta de datos para insumos y resultados que sea común para cada uno de los talleres que se hicieron previamente.

ruta_proyecto <- dirname(getwd())
ruta_clase_01 <- file.path(ruta_proyecto, "Clase_01_GEE")
ruta_clase_02 <- file.path(ruta_proyecto, "Clase_02_R")
ruta_clase_03 <- file.path(ruta_proyecto, "Clase_03_Modelamiento")
ruta_clase_04 <- file.path(ruta_proyecto, "Clase_04_Maxent")
ruta_clase_05 <- file.path(ruta_proyecto, "Clase_05_Metricas")
ruta_clase_06 <- file.path(ruta_proyecto, "Clase_06_amc")
ruta_datos    <- file.path(ruta_proyecto, "datos")

4.1 Cargar y preparar las capas que utilizaremos.

Vamos a cargar los archivos correspondientes a la nueva zona de estudio. Comenzamos con el shapefile que contiene el límite de la zona de interés.

zona_estudio_amc <- sf::st_read(file.path(ruta_datos,"raster_amc/area_de_interes_nm.shp"),quiet=TRUE)

Este shapefile está en formato UTM18S, lo que nos permitirá configurar la resolución en metros. Revisamos la proyección del mismo para ver que se encuentre en metros.

raster::crs(zona_estudio_amc)
## CRS arguments:
##  +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

4.2 Preparar las entidades geográficas con coberturas.

Hemos generado un nuevo conjunto de variables dependientes para la zona de estudio. Las presentaremos a continuación. Comenzaremos con las coberturas separadas (uni clase) que se utilizan para MaxEnt y luego las coberturas juntas (multi clase) para RandomForest. Cargamos las coberturas y comprobamos su proyección que esté en metros.

puntos <- read.csv(file.path(ruta_datos,"raster_amc/cobertura_veg_points2.csv") , header=TRUE, sep=",")
puntos <- sf::st_as_sf(puntos, wkt='WKT')
sf::st_crs(puntos) <- 32718
raster::crs(puntos)
## CRS arguments:
##  +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Hacemos un recuento inicial de las coberturas existentes en la zona de interés.

setNames(aggregate( as.character(WKT) ~ USO, puntos, length),c('USO','Recuento'))
##          USO Recuento
## 1    Abierto       34
## 2   Agricola      290
## 3       Agua       40
## 4 Construido       68
## 5     Nativo      943
## 6 Plantacion      977

Realizamos una selección aleatoria al listado de coberturas presentes según cada tipo de uso. Definimos un tamaño muestral de 34, porque consideramos que es suficiente para realizar la etapa de entrenamiento del MaxEnt.

idx01 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Abierto"    , ] )) , 34 )
idx02 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Agricola"   , ] )) , 34 )
idx03 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Agua"       , ] )) , 34 )
idx04 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Construido" , ] )) , 34 )
idx05 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Nativo"     , ] )) , 34 )
idx06 <- sample( as.numeric(rownames( puntos [ puntos$USO == "Plantacion" , ] )) , 34 )

Aplicamos el índice generado en la selección aleatoria para separar en coberturas independientes desde el data frame original.

puntos_amc_abierto    <- puntos [ idx01 , ]
puntos_amc_agricola   <- puntos [ idx02 , ]
puntos_amc_agua       <- puntos [ idx03 , ]
puntos_amc_construido <- puntos [ idx04 , ]
puntos_amc_nativo     <- puntos [ idx05 , ]
puntos_amc_plantacion <- puntos [ idx06 , ]

Comprobamos las distribución posterior a la selección aleatoria para ver que cumpla con el tamaño muestral definido para cada una de las nuevas coberturas.

rbind(
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_abierto   , length),c('USO','Recuento')),
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_agricola  , length),c('USO','Recuento')),
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_agua      , length),c('USO','Recuento')),
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_construido, length),c('USO','Recuento')),
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_nativo    , length),c('USO','Recuento')),
setNames(aggregate( as.character(WKT) ~ USO, puntos_amc_plantacion, length),c('USO','Recuento'))
)
##          USO Recuento
## 1    Abierto       34
## 2   Agricola       34
## 3       Agua       34
## 4 Construido       34
## 5     Nativo       34
## 6 Plantacion       34

Ahora usamos los mismos puntos pero esta vez con coordenadas en formato XYZ.

puntos_xyz <- read.csv(file.path(ruta_datos,"raster_amc/cobertura_veg_points2_xyz.csv") , header=TRUE, sep=",")

Aplicamos el índice generado en la selección aleatoria para separar en coberturas independientes desde el data frame original.

puntos_xyz_abierto    <- puntos_xyz [ idx01 , ]
puntos_xyz_agricola   <- puntos_xyz [ idx02 , ]
puntos_xyz_agua       <- puntos_xyz [ idx03 , ]
puntos_xyz_construido <- puntos_xyz [ idx04 , ]
puntos_xyz_nativo     <- puntos_xyz [ idx05 , ]
puntos_xyz_plantacion <- puntos_xyz [ idx06 , ]

Volvemos a concatenar las coberturas esta vez muestreadas en caso de necesitarlas más adelante.

puntos_xyz_muestreados <-
rbind(
puntos_xyz_abierto   ,
puntos_xyz_agricola  ,
puntos_xyz_agua      ,
puntos_xyz_construido,
puntos_xyz_nativo    ,
puntos_xyz_plantacion
)

Lo convertimos en objeto sf para exportarlo como Shapefile.

puntos_xyz_muestreados <- st_as_sf( puntos_xyz_muestreados, coords=c('X','Y'))
st_crs( puntos_xyz_muestreados) <- 32718

Grabamos un archivo en disco para poder inspeccionar con QGis de ser necesario.

st_write(puntos_xyz_muestreados,
         file.path(ruta_datos,"raster_amc/cobertura_veg_points2_sample.shp"),
         delete_layer=TRUE,
         quiet=TRUE)

4.3 Preparar nuevas imágenes para el Análisis Multicriterio.

Ahora cargaremos los archivos raster que utilizaremos para el ejercicio. Estas imágenes fueron generadas para la zona de interés con GEE usando código Javascript. Posteriormente se procesaron con la librería GDAL, pasos tales como, reproyección, máscara y remuestreo. Listamos las imágenes GeoTiff con las que vamos a trabajar.

archivos_raster_amc <- list.files(file.path(ruta_datos,"raster_amc/Bandas_Landsat"),
                       pattern="_crop_32718.tif$|aspect.tif$|slope.tif$|dem.tif$")
archivos_raster_amc <- archivos_raster_amc[c(2:7,9,1,8,10)]
as.data.frame(archivos_raster_amc)
##    archivos_raster_amc
## 1    B2_crop_32718.tif
## 2    B3_crop_32718.tif
## 3    B4_crop_32718.tif
## 4    B5_crop_32718.tif
## 5    B6_crop_32718.tif
## 6    B7_crop_32718.tif
## 7  NDVI_crop_32718.tif
## 8           aspect.tif
## 9              dem.tif
## 10           slope.tif

Vamos a trabajar con las bandas B2, a B6, con NDVI y las demás las dejaremos para pasos posteriores. Con esas bandas, creamos un archivo stack usando la librería GDAL y la cargamos directamente en R desde nuestro directorio.

stack_bandas_amc <- raster::stack( file.path(ruta_datos,"raster_amc/Bandas_Landsat/stack_bandas_amc_32718.tif" ))
names(stack_bandas_amc) <- c('B2','B3','B4','B5','B6','NDVI')

4.3.1 Mapas temáticos mediante método MaxEnt.

Vamos a aplicar de nuevo el procedimiento de MaxEnt a las nuevas imágenes. El propósito es generar con el algoritmo MaxEnt un tipo de mapa temático que se use como insumo para aplicar el análisis multicriterio.

Retiramos la variable USO de cada capa de cobertura para trabajar únicamente con las coordenadas.

puntos_xyz_abierto   $USO <- NULL
puntos_xyz_agricola  $USO <- NULL
puntos_xyz_agua      $USO <- NULL
puntos_xyz_construido$USO <- NULL
puntos_xyz_nativo    $USO <- NULL
puntos_xyz_plantacion$USO <- NULL

Creamos grupos sintéticos para particionar posteriormente cada conjunto de datos.

grupos_amc_abierto    <- dismo::kfold( puntos_xyz_abierto   , k=5)
grupos_amc_agricola   <- dismo::kfold( puntos_xyz_agricola  , k=5)
grupos_amc_agua       <- dismo::kfold( puntos_xyz_agua      , k=5)
grupos_amc_construido <- dismo::kfold( puntos_xyz_construido, k=5)
grupos_amc_nativo     <- dismo::kfold( puntos_xyz_nativo    , k=5)
grupos_amc_plantacion <- dismo::kfold( puntos_xyz_plantacion, k=5)

Particionamos cada conjunto de datos para crear una muestra de entrenamiento y una de prueba respectivamente.

puntos_xyz_train_abierto    <- puntos_xyz_abierto    [ grupos_amc_abierto    != 1 , ]
puntos_xyz_train_agricola   <- puntos_xyz_agricola   [ grupos_amc_agricola   != 1 , ]
puntos_xyz_train_agua       <- puntos_xyz_agua       [ grupos_amc_agua       != 1 , ]
puntos_xyz_train_construido <- puntos_xyz_construido [ grupos_amc_construido != 1 , ]
puntos_xyz_train_nativo     <- puntos_xyz_nativo     [ grupos_amc_nativo     != 1 , ]
puntos_xyz_train_plantacion <- puntos_xyz_plantacion [ grupos_amc_plantacion != 1 , ]
puntos_xyz_teste_abierto    <- puntos_xyz_abierto    [ grupos_amc_abierto    == 1 , ]
puntos_xyz_teste_agricola   <- puntos_xyz_agricola   [ grupos_amc_agricola   == 1 , ]
puntos_xyz_teste_agua       <- puntos_xyz_agua       [ grupos_amc_agua       == 1 , ]
puntos_xyz_teste_construido <- puntos_xyz_construido [ grupos_amc_construido == 1 , ]
puntos_xyz_teste_nativo     <- puntos_xyz_nativo     [ grupos_amc_nativo     == 1 , ]
puntos_xyz_teste_plantacion <- puntos_xyz_plantacion [ grupos_amc_plantacion == 1 , ]

Ya tenemos los necesario para entrenar cada cobertura usando MaxEnt. Usamos el argumento path para que la ejecución sea más rápida a partir de la primera corrida con lo cual se almacenan varios archivos en disco que el algoritmo crea cada vez que corre.

maxent_xyz_train_abierto    <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_abierto    , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_abierto"   ), silent=FALSE)
## Loading required namespace: rJava
## This is MaxEnt version 3.4.1
maxent_xyz_train_agricola   <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_agricola   , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_agricola"  ), silent=FALSE)
## This is MaxEnt version 3.4.1
maxent_xyz_train_agua       <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_agua       , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_agua"      ), silent=FALSE)
## This is MaxEnt version 3.4.1
maxent_xyz_train_construido <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_construido , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_construido"), silent=FALSE)
## Warning in .local(x, p, ...): 1 (3.7%) of the presence points have NA predictor values
## This is MaxEnt version 3.4.1
maxent_xyz_train_nativo     <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_nativo     , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_nativo"    ), silent=FALSE)
## This is MaxEnt version 3.4.1
maxent_xyz_train_plantacion <- dismo::maxent( x=stack_bandas_amc, p=puntos_xyz_train_plantacion , ngb=2000, cache=TRUE, threads=4, memory_assigned=3600, skipifexists=TRUE, verbose=FALSE, path=file.path(ruta_datos,"maxent_xyz/train_plantacion"), silent=FALSE)
## Warning in .local(x, p, ...): 1 (3.7%) of the presence points have NA predictor values
## This is MaxEnt version 3.4.1

Graficamos la contribución de variables por cada cobertura.

#par(mfrow=c(2,3))
#plot(maxent_xyz_train_abierto    , main="Variable contribution\nAbierto"    )
#plot(maxent_xyz_train_agricola   , main="Variable contribution\nAgricola"   )
#plot(maxent_xyz_train_agua       , main="Variable contribution\nAgua"       )
#plot(maxent_xyz_train_construido , main="Variable contribution\nConstruido" )
#plot(maxent_xyz_train_nativo     , main="Variable contribution\nNativo"     )
#plot(maxent_xyz_train_plantacion , main="Variable contribution\nPlantacion" )
Contribución de variables al modelo MaxEnt.

Contribución de variables al modelo MaxEnt.

Obtenemos raster con las predicciones para cada cobertura.

maxent_xyz_predict_abierto    <- dismo::predict(maxent_xyz_train_abierto    , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_abierto.tif"   ), overwrite=TRUE)
## Loading required namespace: rJava
## This is MaxEnt version 3.4.1
maxent_xyz_predict_agricola   <- dismo::predict(maxent_xyz_train_agricola   , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_agricola.tif"  ), overwrite=TRUE)
## This is MaxEnt version 3.4.1
maxent_xyz_predict_agua       <- dismo::predict(maxent_xyz_train_agua       , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_agua.tif"      ), overwrite=TRUE)
## This is MaxEnt version 3.4.1
maxent_xyz_predict_construido <- dismo::predict(maxent_xyz_train_construido , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_construido.tif"), overwrite=TRUE)
## This is MaxEnt version 3.4.1
maxent_xyz_predict_nativo     <- dismo::predict(maxent_xyz_train_nativo     , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_nativo.tif"    ), overwrite=TRUE)
## This is MaxEnt version 3.4.1
maxent_xyz_predict_plantacion <- dismo::predict(maxent_xyz_train_plantacion , stack_bandas_amc, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_predict_plantacion.tif"), overwrite=TRUE)
## This is MaxEnt version 3.4.1

Generamos puntos background de testeo para poder evaluar las predicciones.

puntos_xyz_bg_teste      <- dismo::randomPoints( stack_bandas_amc, 1000)

Evaluamos cada modelo usando los mismos puntos de background generados previamente.

maxent_xyz_eval_abierto    <- dismo::evaluate(maxent_xyz_train_abierto   , p=puntos_xyz_train_abierto   , a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1
maxent_xyz_eval_agricola   <- dismo::evaluate(maxent_xyz_train_agricola  , p=puntos_xyz_train_agricola  , a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1
maxent_xyz_eval_agua       <- dismo::evaluate(maxent_xyz_train_agua      , p=puntos_xyz_train_agua      , a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1
maxent_xyz_eval_construido <- dismo::evaluate(maxent_xyz_train_construido, p=puntos_xyz_train_construido, a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1
maxent_xyz_eval_nativo     <- dismo::evaluate(maxent_xyz_train_nativo    , p=puntos_xyz_train_nativo    , a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1
maxent_xyz_eval_plantacion <- dismo::evaluate(maxent_xyz_train_plantacion, p=puntos_xyz_train_plantacion, a=puntos_xyz_bg_teste, x=stack_bandas_amc)
## This is MaxEnt version 3.4.1 
## This is MaxEnt version 3.4.1

Podemos graficar las curvas ROC por cada cobertura.

#par(mfrow=c(2,3))
#plot(maxent_xyz_eval_abierto    , 'ROC')
#plot(maxent_xyz_eval_agricola   , 'ROC')
#plot(maxent_xyz_eval_agua       , 'ROC')
#plot(maxent_xyz_eval_construido , 'ROC')
#plot(maxent_xyz_eval_nativo     , 'ROC')
#plot(maxent_xyz_eval_plantacion , 'ROC')
Curvas ROC de cada variables en el modelo MaxEnt.

Curvas ROC de cada variables en el modelo MaxEnt.

Extraemos los valores del estadígrafo MaxKappa.

maxkpp_xyz_eval_abierto    <- maxent_xyz_eval_abierto   @t[which.max(maxent_xyz_eval_abierto   @kappa)]
maxkpp_xyz_eval_agricola   <- maxent_xyz_eval_agricola  @t[which.max(maxent_xyz_eval_agricola  @kappa)]
maxkpp_xyz_eval_agua       <- maxent_xyz_eval_agua      @t[which.max(maxent_xyz_eval_agua      @kappa)]
maxkpp_xyz_eval_construido <- maxent_xyz_eval_construido@t[which.max(maxent_xyz_eval_construido@kappa)]
maxkpp_xyz_eval_nativo     <- maxent_xyz_eval_nativo    @t[which.max(maxent_xyz_eval_nativo    @kappa)]
maxkpp_xyz_eval_plantacion <- maxent_xyz_eval_plantacion@t[which.max(maxent_xyz_eval_plantacion@kappa)]

Extraemos los valores del estadígrafo TSS.

maxtss_xyz_eval_abierto    <- maxent_xyz_eval_abierto   @t[which.max(maxent_xyz_eval_abierto   @TPR+(maxent_xyz_eval_abierto   @TNR-1))]
maxtss_xyz_eval_agricola   <- maxent_xyz_eval_agricola  @t[which.max(maxent_xyz_eval_agricola  @TPR+(maxent_xyz_eval_agricola  @TNR-1))]
maxtss_xyz_eval_agua       <- maxent_xyz_eval_agua      @t[which.max(maxent_xyz_eval_agua      @TPR+(maxent_xyz_eval_agua      @TNR-1))]
maxtss_xyz_eval_construido <- maxent_xyz_eval_construido@t[which.max(maxent_xyz_eval_construido@TPR+(maxent_xyz_eval_construido@TNR-1))]
maxtss_xyz_eval_nativo     <- maxent_xyz_eval_nativo    @t[which.max(maxent_xyz_eval_nativo    @TPR+(maxent_xyz_eval_nativo    @TNR-1))]
maxtss_xyz_eval_plantacion <- maxent_xyz_eval_plantacion@t[which.max(maxent_xyz_eval_plantacion@TPR+(maxent_xyz_eval_plantacion@TNR-1))]

Promediamos ambas métricas para reclasificar cada imagen en función de ese valor promedio.

maxprom_xyz_eval_abierto    <- (maxkpp_xyz_eval_abierto   +maxtss_xyz_eval_abierto   )/2
maxprom_xyz_eval_agricola   <- (maxkpp_xyz_eval_agricola  +maxtss_xyz_eval_agricola  )/2
maxprom_xyz_eval_agua       <- (maxkpp_xyz_eval_agua      +maxtss_xyz_eval_agua      )/2
maxprom_xyz_eval_construido <- (maxkpp_xyz_eval_construido+maxtss_xyz_eval_construido)/2
maxprom_xyz_eval_nativo     <- (maxkpp_xyz_eval_nativo    +maxtss_xyz_eval_nativo    )/2
maxprom_xyz_eval_plantacion <- (maxkpp_xyz_eval_plantacion+maxtss_xyz_eval_plantacion)/2

Definimos una función que reclasifique cada imagen según un valor de corte suministrado, en este caso, el promedio entre el MaxKappa y el MaxTSS.

funReclValProm <- function(...,strval) {
ifelse(... <= strval, 0,
ifelse(... >  strval, 1, NA)) }

Luego llamamos dicha función con los argumentos que corresponde en cada caso.

maxent_xyz_reclass_abierto    <- calc(maxent_xyz_predict_abierto   , function(x){ funReclValProm(x,strval=maxprom_xyz_eval_abierto   )}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_abierto.tif"   ), overwrite=TRUE)
maxent_xyz_reclass_agricola   <- calc(maxent_xyz_predict_agricola  , function(x){ funReclValProm(x,strval=maxprom_xyz_eval_agricola  )}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_agricola.tif"  ), overwrite=TRUE)
maxent_xyz_reclass_agua       <- calc(maxent_xyz_predict_agua      , function(x){ funReclValProm(x,strval=maxprom_xyz_eval_agua      )}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_agua.tif"      ), overwrite=TRUE)
maxent_xyz_reclass_construido <- calc(maxent_xyz_predict_construido, function(x){ funReclValProm(x,strval=maxprom_xyz_eval_construido)}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_construido.tif"), overwrite=TRUE)
maxent_xyz_reclass_nativo     <- calc(maxent_xyz_predict_nativo    , function(x){ funReclValProm(x,strval=maxprom_xyz_eval_nativo    )}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_nativo.tif"    ), overwrite=TRUE)
maxent_xyz_reclass_plantacion <- calc(maxent_xyz_predict_plantacion, function(x){ funReclValProm(x,strval=maxprom_xyz_eval_plantacion)}, filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_plantacion.tif"), overwrite=TRUE)

Graficamos las imágenes temáticas.

png(file.path(ruta_datos,"figs/graficar_reclass_valor_cobertura-1.png"), units="in", width=8.45, height=6, res=300)
par(mfrow=c(2,3))
plot(maxent_xyz_reclass_abierto   ,main="Reclass\nAbierto   ")
plot(maxent_xyz_reclass_agricola  ,main="Reclass\nAgricola  ")
plot(maxent_xyz_reclass_agua      ,main="Reclass\nAgua      ")
plot(maxent_xyz_reclass_construido,main="Reclass\nConstruido")
plot(maxent_xyz_reclass_nativo    ,main="Reclass\nNativo    ")
plot(maxent_xyz_reclass_plantacion,main="Reclass\nPlantacion")
box()
dev.off()
## png 
##   2
Coberturas reclasificadas por Max Kappa y Max TSS.

Coberturas reclasificadas por Max Kappa y Max TSS.

4.3.2 Generar grilla para la zona de estudio.

Generamos un grillado usando el software QGis para poder dividir la zona de estudio en parches de menor tamaño. Encontramos que nuestra zona de interés posee una grilla que no son cuadrados regulares, sino que irregulares, por lo que la grilla generada por nosotros no se ajustó exactamente al área, quedando los parches de los bordes de un tamaño inferior.

Rejilla con cuadrantes para la zona de estudio.

Rejilla con cuadrantes para la zona de estudio.

Vamos a generar algunas métricas de paisaje para obtener coberturas representativas para identificar zonas con mayor propensión a incendios forestales.

if(!require(landscapemetrics)) suppressMessages( suppressWarnings( install.packages("landscapemetrics") ))

Verificamos que las capas temáticas resultantes de MaxEnt posean las propiedades necesarias.

check_landscape(maxent_xyz_reclass_abierto   )
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
check_landscape(maxent_xyz_reclass_agricola  )
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
check_landscape(maxent_xyz_reclass_agua      )
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
check_landscape(maxent_xyz_reclass_construido)
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
check_landscape(maxent_xyz_reclass_nativo    )
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
check_landscape(maxent_xyz_reclass_plantacion)
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔

Resumimos las áreas totales.

rbind(
lsm_l_ta( maxent_xyz_reclass_abierto    ),
lsm_l_ta( maxent_xyz_reclass_agricola   ),
lsm_l_ta( maxent_xyz_reclass_agua       ),
lsm_l_ta( maxent_xyz_reclass_construido ),
lsm_l_ta( maxent_xyz_reclass_nativo     ),
lsm_l_ta( maxent_xyz_reclass_plantacion )
)
## # A tibble: 6 x 6
##   layer level     class    id metric   value
##   <int> <chr>     <int> <int> <chr>    <dbl>
## 1     1 landscape    NA    NA ta     671141.
## 2     1 landscape    NA    NA ta     671141.
## 3     1 landscape    NA    NA ta     671141.
## 4     1 landscape    NA    NA ta     671141.
## 5     1 landscape    NA    NA ta     671141.
## 6     1 landscape    NA    NA ta     671141.

Resumimos los bordes totales.

rbind(
lsm_l_te( maxent_xyz_reclass_abierto    ),
lsm_l_te( maxent_xyz_reclass_agricola   ),
lsm_l_te( maxent_xyz_reclass_agua       ),
lsm_l_te( maxent_xyz_reclass_construido ),
lsm_l_te( maxent_xyz_reclass_nativo     ),
lsm_l_te( maxent_xyz_reclass_plantacion )
)
## # A tibble: 6 x 6
##   layer level     class    id metric     value
##   <int> <chr>     <int> <int> <chr>      <dbl>
## 1     1 landscape    NA    NA te     41657375.
## 2     1 landscape    NA    NA te     35474272.
## 3     1 landscape    NA    NA te     15178427.
## 4     1 landscape    NA    NA te      7722284.
## 5     1 landscape    NA    NA te     30634735.
## 6     1 landscape    NA    NA te     47951068.

Resumimos índique de complejidad de forma, a mayor valor, mayor complejidad.

rbind(
lsm_l_np( maxent_xyz_reclass_abierto    ),
lsm_l_np( maxent_xyz_reclass_agricola   ),
lsm_l_np( maxent_xyz_reclass_agua       ),
lsm_l_np( maxent_xyz_reclass_construido ),
lsm_l_np( maxent_xyz_reclass_nativo     ),
lsm_l_np( maxent_xyz_reclass_plantacion )
)
## # A tibble: 6 x 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA np     43110
## 2     1 landscape    NA    NA np     43484
## 3     1 landscape    NA    NA np     37309
## 4     1 landscape    NA    NA np     16696
## 5     1 landscape    NA    NA np     54597
## 6     1 landscape    NA    NA np     58535

Resumimos el índice de Shannon que indica mayor diversidad a mayor valor del índice.

rbind(
lsm_l_shdi( maxent_xyz_reclass_abierto    ),
lsm_l_shdi( maxent_xyz_reclass_agricola   ),
lsm_l_shdi( maxent_xyz_reclass_agua       ),
lsm_l_shdi( maxent_xyz_reclass_construido ),
lsm_l_shdi( maxent_xyz_reclass_nativo     ),
lsm_l_shdi( maxent_xyz_reclass_plantacion )
)
## # A tibble: 6 x 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA shdi   0.616
## 2     1 landscape    NA    NA shdi   0.490
## 3     1 landscape    NA    NA shdi   0.210
## 4     1 landscape    NA    NA shdi   0.143
## 5     1 landscape    NA    NA shdi   0.369
## 6     1 landscape    NA    NA shdi   0.645

Resumimos a nivel de porcentaje del paisaje por cada clase.

rbind(
lsm_c_pland( maxent_xyz_reclass_abierto    ),
lsm_c_pland( maxent_xyz_reclass_agricola   ),
lsm_c_pland( maxent_xyz_reclass_agua       ),
lsm_c_pland( maxent_xyz_reclass_construido ),
lsm_c_pland( maxent_xyz_reclass_nativo     ),
lsm_c_pland( maxent_xyz_reclass_plantacion )
)
## # A tibble: 12 x 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 class     0    NA pland  69.4 
##  2     1 class     1    NA pland  30.6 
##  3     1 class     0    NA pland  80.7 
##  4     1 class     1    NA pland  19.3 
##  5     1 class     0    NA pland  94.6 
##  6     1 class     1    NA pland   5.39
##  7     1 class     0    NA pland  96.7 
##  8     1 class     1    NA pland   3.25
##  9     1 class     0    NA pland  87.9 
## 10     1 class     1    NA pland  12.1 
## 11     1 class     0    NA pland  65.4 
## 12     1 class     1    NA pland  34.6

Resumimos a nivel de número de parches por cada clase.

rbind(
lsm_c_np( maxent_xyz_reclass_abierto    ),
lsm_c_np( maxent_xyz_reclass_agricola   ),
lsm_c_np( maxent_xyz_reclass_agua       ),
lsm_c_np( maxent_xyz_reclass_construido ),
lsm_c_np( maxent_xyz_reclass_nativo     ),
lsm_c_np( maxent_xyz_reclass_plantacion )
)
## # A tibble: 12 x 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 class     0    NA np     10295
##  2     1 class     1    NA np     32815
##  3     1 class     0    NA np      7333
##  4     1 class     1    NA np     36151
##  5     1 class     0    NA np       881
##  6     1 class     1    NA np     36428
##  7     1 class     0    NA np       832
##  8     1 class     1    NA np     15864
##  9     1 class     0    NA np      3971
## 10     1 class     1    NA np     50626
## 11     1 class     0    NA np     17374
## 12     1 class     1    NA np     41161

Calcularemos algunas métricas por cada zona de la grilla generada para el estudio.

rejilla_nm <- st_read(file.path(ruta_datos,"zip/6x6_rectangle_angle.shp"),quiet=TRUE)
grilla_01_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[1,] )  ; dsp_01_xyz_reclass_abierto <- lsm_l_shdi( grilla_01_xyz_reclass_abierto )
grilla_02_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[2,] )  ; dsp_02_xyz_reclass_abierto <- lsm_l_shdi( grilla_02_xyz_reclass_abierto )
grilla_03_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[3,] )  ; dsp_03_xyz_reclass_abierto <- lsm_l_shdi( grilla_03_xyz_reclass_abierto )
grilla_04_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[4,] )  ; dsp_04_xyz_reclass_abierto <- lsm_l_shdi( grilla_04_xyz_reclass_abierto )
grilla_05_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[5,] )  ; dsp_05_xyz_reclass_abierto <- lsm_l_shdi( grilla_05_xyz_reclass_abierto )
grilla_06_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[6,] )  ; dsp_06_xyz_reclass_abierto <- lsm_l_shdi( grilla_06_xyz_reclass_abierto )
grilla_07_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[7,] )  ; dsp_07_xyz_reclass_abierto <- lsm_l_shdi( grilla_07_xyz_reclass_abierto )
grilla_08_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[8,] )  ; dsp_08_xyz_reclass_abierto <- lsm_l_shdi( grilla_08_xyz_reclass_abierto )
grilla_09_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[9,] )  ; dsp_09_xyz_reclass_abierto <- lsm_l_shdi( grilla_09_xyz_reclass_abierto )
grilla_10_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[10,] ) ; dsp_10_xyz_reclass_abierto <- lsm_l_shdi( grilla_10_xyz_reclass_abierto )
grilla_11_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[11,] ) ; dsp_11_xyz_reclass_abierto <- lsm_l_shdi( grilla_11_xyz_reclass_abierto )
grilla_12_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[12,] ) ; dsp_12_xyz_reclass_abierto <- lsm_l_shdi( grilla_12_xyz_reclass_abierto )
grilla_13_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[13,] ) ; dsp_13_xyz_reclass_abierto <- lsm_l_shdi( grilla_13_xyz_reclass_abierto )
grilla_14_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[14,] ) ; dsp_14_xyz_reclass_abierto <- lsm_l_shdi( grilla_14_xyz_reclass_abierto )
grilla_15_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[15,] ) ; dsp_15_xyz_reclass_abierto <- lsm_l_shdi( grilla_15_xyz_reclass_abierto )
grilla_16_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[16,] ) ; dsp_16_xyz_reclass_abierto <- lsm_l_shdi( grilla_16_xyz_reclass_abierto )
grilla_17_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[17,] ) ; dsp_17_xyz_reclass_abierto <- lsm_l_shdi( grilla_17_xyz_reclass_abierto )
grilla_18_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[18,] ) ; dsp_18_xyz_reclass_abierto <- lsm_l_shdi( grilla_18_xyz_reclass_abierto )
grilla_19_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[19,] ) ; dsp_19_xyz_reclass_abierto <- lsm_l_shdi( grilla_19_xyz_reclass_abierto )
grilla_20_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[20,] ) ; dsp_20_xyz_reclass_abierto <- lsm_l_shdi( grilla_20_xyz_reclass_abierto )
grilla_21_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[21,] ) ; dsp_21_xyz_reclass_abierto <- lsm_l_shdi( grilla_21_xyz_reclass_abierto )
grilla_22_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[22,] ) ; dsp_22_xyz_reclass_abierto <- lsm_l_shdi( grilla_22_xyz_reclass_abierto )
grilla_23_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[23,] ) ; dsp_23_xyz_reclass_abierto <- lsm_l_shdi( grilla_23_xyz_reclass_abierto )
grilla_24_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[24,] ) ; dsp_24_xyz_reclass_abierto <- lsm_l_shdi( grilla_24_xyz_reclass_abierto )
grilla_25_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[25,] ) ; dsp_25_xyz_reclass_abierto <- lsm_l_shdi( grilla_25_xyz_reclass_abierto )
grilla_26_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[26,] ) ; dsp_26_xyz_reclass_abierto <- lsm_l_shdi( grilla_26_xyz_reclass_abierto )
grilla_27_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[27,] ) ; dsp_27_xyz_reclass_abierto <- lsm_l_shdi( grilla_27_xyz_reclass_abierto )
grilla_28_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[28,] ) ; dsp_28_xyz_reclass_abierto <- lsm_l_shdi( grilla_28_xyz_reclass_abierto )
grilla_29_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[29,] ) ; dsp_29_xyz_reclass_abierto <- lsm_l_shdi( grilla_29_xyz_reclass_abierto )
grilla_30_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[30,] ) ; dsp_30_xyz_reclass_abierto <- lsm_l_shdi( grilla_30_xyz_reclass_abierto )
grilla_31_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[31,] ) ; dsp_31_xyz_reclass_abierto <- lsm_l_shdi( grilla_31_xyz_reclass_abierto )
grilla_32_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[32,] ) ; dsp_32_xyz_reclass_abierto <- lsm_l_shdi( grilla_32_xyz_reclass_abierto )
grilla_33_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[33,] ) ; dsp_33_xyz_reclass_abierto <- lsm_l_shdi( grilla_33_xyz_reclass_abierto )
grilla_34_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[34,] ) ; dsp_34_xyz_reclass_abierto <- lsm_l_shdi( grilla_34_xyz_reclass_abierto )
grilla_35_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[35,] ) ; dsp_35_xyz_reclass_abierto <- lsm_l_shdi( grilla_35_xyz_reclass_abierto )
grilla_36_xyz_reclass_abierto <- crop( maxent_xyz_reclass_abierto , rejilla_nm[36,] ) ; dsp_36_xyz_reclass_abierto <- lsm_l_shdi( grilla_36_xyz_reclass_abierto )

Consolidamos cada valor por rejilla en una sola capa.

shannon_xyz_reclass_abierto <- rbind(
dsp_01_xyz_reclass_abierto,
dsp_02_xyz_reclass_abierto,
dsp_03_xyz_reclass_abierto,
dsp_04_xyz_reclass_abierto,
dsp_05_xyz_reclass_abierto,
dsp_06_xyz_reclass_abierto,
dsp_07_xyz_reclass_abierto,
dsp_08_xyz_reclass_abierto,
dsp_09_xyz_reclass_abierto,
dsp_10_xyz_reclass_abierto,
dsp_11_xyz_reclass_abierto,
dsp_12_xyz_reclass_abierto,
dsp_13_xyz_reclass_abierto,
dsp_14_xyz_reclass_abierto,
dsp_15_xyz_reclass_abierto,
dsp_16_xyz_reclass_abierto,
dsp_17_xyz_reclass_abierto,
dsp_18_xyz_reclass_abierto,
dsp_19_xyz_reclass_abierto,
dsp_20_xyz_reclass_abierto,
dsp_21_xyz_reclass_abierto,
dsp_22_xyz_reclass_abierto,
dsp_23_xyz_reclass_abierto,
dsp_24_xyz_reclass_abierto,
dsp_25_xyz_reclass_abierto,
dsp_26_xyz_reclass_abierto,
dsp_27_xyz_reclass_abierto,
dsp_28_xyz_reclass_abierto,
dsp_29_xyz_reclass_abierto,
dsp_30_xyz_reclass_abierto,
dsp_31_xyz_reclass_abierto,
dsp_32_xyz_reclass_abierto,
dsp_33_xyz_reclass_abierto,
dsp_34_xyz_reclass_abierto,
dsp_35_xyz_reclass_abierto,
dsp_36_xyz_reclass_abierto)
id_rejilla <- c(
paste0('A',seq(1:6)),
paste0('B',seq(1:6)),
paste0('C',seq(1:6)),
paste0('D',seq(1:6)),
paste0('E',seq(1:6)),
paste0('F',seq(1:6))
)

Reordenamos las columnas para poder graficar el nuevo indicador de Shannon a nivel de cada polígono.

shannon_xyz_reclass_abierto <- cbind(NAME=id_rejilla,shannon_xyz_reclass_abierto[,6])
shannon_xyz_rejilla_abierto <- merge(rejilla_nm, shannon_xyz_reclass_abierto,by='NAME')
png(file.path(ruta_datos,"figs/plot_concatenar_dsp_xyz_cobertura.png"), units="in", width=6, height=6, res=300)
plot(shannon_xyz_rejilla_abierto[,3],main="Índice Shannon")
dev.off()
## png 
##   2
Índice de Shannon obtenido por cada segmento.

Índice de Shannon obtenido por cada segmento.

Para completar neustras capas, vamos a rasterizar esta capa como Diversidad.

maxent_xyz_reclass_diversidad <- rasterize( shannon_xyz_rejilla_abierto ,
filename=file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_diversidad.tif"),
maxent_xyz_reclass_abierto ,
field='value' ,
silent=TRUE,
overwrite=TRUE)

Vamos a cargar un raster stack que hemos preparado con las nuevas variables que usaremos para el análisis multicriterio.

stack_xyz_predict_amc <- raster::brick(file.path(ruta_datos,"raster_xyz/stack_xyz_predict_amc.tif"))
names(stack_xyz_predict_amc) <- c(
"aspect"    ,
"construido",
"diversidad",
"nativo"    ,
"plantacion",
"slope"     )

Revisamos las variables con las que vamos a trabajar en el análisis multicriterio.

variables_amc <- data.frame(variables_amc=ls(pattern="maxent_xyz_reclass_"))
variables_amc <- data.frame(variables_amc=variables_amc[c(4:9),])
variables_amc
##                   variables_amc
## 1     maxent_xyz_reclass_aspect
## 2 maxent_xyz_reclass_construido
## 3 maxent_xyz_reclass_diversidad
## 4     maxent_xyz_reclass_nativo
## 5 maxent_xyz_reclass_plantacion
## 6      maxent_xyz_reclass_slope

4.4 Generar el análisis multicriterio (AMC).

Primero veamos la escala de nuestros datos para ver la necesidad de transformarlos en indicadores compatibles. Veamos cuál es la escala de valores de nuestras variables.

¿Son compatibles/conmensurables nuestros datos? En otras palabras, ¿son operacionalmente compatibles? Por lo tanto: Lo primero dentro de un AMC es normalizar los datos de las variables en indicadores compatibles.

En este caso vamos a trabajar con las siguientes 6 variables:

  • Capa Maxent Construido
 ==> Mayor nivel de construcción, Mayor densidad poblacional y probabilidad de incendio.
  • Capa Maxent Nativo
 ==> Mayor nivel de bosque nativo, Mayor biomasa disponible para incendio.
  • Capa Maxent Plantación
 ==> Mayor nivel de plantación forestal, Mayor biomasa disponible para incendio.
  • Diversidad de Paisaje
 ==> Mayor nivel de diversidad de paisaje, Menor probabilidad de que el incendio se propague.
  • Pendiente terreno
 ==> Mayor nivel de Pendiente, Mayor probabilidad de vientos ascendentes y propagación de incendio.
  • Aspecto terreno
 ==> Aquellos lugares con exposición norte y este debiesen tener mayor radiación y temperatura.

En este práctico normalizaremos todas nuestras variables en el rango de 0 a 1 en base a distintas reglas. Para todos nuestros indicadores, valores mayores indicarán mayor probabilidad de incendio de gran magnitud.

Revisamos las estadísticas de cada cobertura.

summa_xyz_predict_amc <-
data.frame(rbind(
stack_xyz_predict_amc@data@max   ,
stack_xyz_predict_amc@data@min
))
colnames(summa_xyz_predict_amc) <-
stack_xyz_predict_amc@data@names
summa_xyz_predict_amc
##     aspect  construido diversidad       nativo  plantacion    slope
## 1 359.7582 1.000000000  0.6828299 0.9176571369 0.873118997 60.65885
## 2  -1.0000 0.003941121  0.3597561 0.0008257756 0.008149054  0.00000

4.5 Normalización de las variables.

Transformación de datos espaciales en indicadores compatibles.

CAPAS MAXENT: Las capas probabilísticas de las clases generadas en el práctico Maxent vienen en valores entre 0 a 1, por lo que no es necesario normalizar. Pero sólo como un ejemplo práctico, lo que haremos será utilizar el cuadrado del valor original para maximizar diferencias entre valores extremos.

#maxent_xyz_normals_aspect      <- subset(stack_xyz_predict_amc,1)^2 # aspect
 maxent_xyz_normals_construido  <- subset(stack_xyz_predict_amc,2)^2 # construido
#maxent_xyz_normals_diversidad  <- subset(stack_xyz_predict_amc,3)^2 # diversidad
 maxent_xyz_normals_nativo      <- subset(stack_xyz_predict_amc,4)^2 # nativo
 maxent_xyz_normals_plantacion  <- subset(stack_xyz_predict_amc,5)^2 # plantacion
#maxent_xyz_normals_slope       <- subset(stack_xyz_predict_amc,6)^2 # slope
#writeRaster(maxent_xyz_normals_aspect    ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"    ),overwrite=TRUE)
 writeRaster(maxent_xyz_normals_construido,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_construido.tif"),overwrite=TRUE)
#writeRaster(maxent_xyz_normals_diversidad,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_diversidad.tif"),overwrite=TRUE)
 writeRaster(maxent_xyz_normals_nativo    ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_nativo.tif"    ),overwrite=TRUE)
 writeRaster(maxent_xyz_normals_plantacion,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_plantacion.tif"),overwrite=TRUE)
#writeRaster(maxent_xyz_normals_slope     ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"     ),overwrite=TRUE)

NATIVO: Veamos como cambiaron nuestros datos, antes y después de normalizar.

maxent_xyz_normals_nativo@data@min
## [1] 6.819053e-07
maxent_xyz_normals_nativo@data@max
## [1] 0.8420946
maxent_xyz_normals_nativo <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_nativo.tif"))
png(file.path(ruta_datos,"figs/plot_maxent_xyz_normals_nativo.png"), units="px", width=1344, height=960, pointsize=24)
plot(maxent_xyz_normals_nativo, main="Variable nativo\nnormalizada")
dev.off()
## png 
##   2
Variable nativo normalizada en rango 0 a 1..

Variable nativo normalizada en rango 0 a 1..

png(file.path(ruta_datos,"figs/hist_minmax_xyz_normals_nativo.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
hist( maxent_xyz_predict_nativo, breaks=100, main = "Original")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
hist( maxent_xyz_normals_nativo, breaks=100, main = "Normalizado")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
dev.off()
## png 
##   2
Variable nativo original versus normalizada.

Variable nativo original versus normalizada.

DIVERSIDAD PAISAJE: Utilizaremos una aproximación de reescalamiento de tipo "max-min": Y = (x-min)/(max-min).

maxent_xyz_normals_diversidad <-
  (subset(stack_xyz_predict_amc,3) - summa_xyz_predict_amc[2,3]) /
  (summa_xyz_predict_amc[1,3]      - summa_xyz_predict_amc[2,3])
 writeRaster(maxent_xyz_normals_diversidad,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_diversidad.tif"),overwrite=TRUE)

Veamos si la variable diversidad quedó normalizada.

maxent_xyz_normals_diversidad@data@min
## [1] -2.920972e-15
maxent_xyz_normals_diversidad@data@max
## [1] 1
maxent_xyz_normals_diversidad <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_diversidad.tif"))
png(file.path(ruta_datos,"figs/plot_maxent_xyz_normals_diversidad.png"), units="px", width=1344, height=960, pointsize=24)
plot(maxent_xyz_normals_diversidad, main="Variable diversidad\nnormalizada")
dev.off()
## png 
##   2
Variable diversidad normalizada en rango 0 a 1..

Variable diversidad normalizada en rango 0 a 1..

png(file.path(ruta_datos,"figs/hist_minmax_xyz_normals_diversidad.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
hist( maxent_xyz_reclass_diversidad, breaks=100, main = "Original")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
hist( maxent_xyz_normals_diversidad, breaks=100, main = "Normalizado")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
dev.off()
## png 
##   2
Variable diversidad original versus normalizada.

Variable diversidad original versus normalizada.

CAPAS TERRENO: Para estas capas tomaremos dos aproximaciones diferentes para cada una.

PENDIENTE: Esta capa tiene valores entre 0 y 75. Para normalizar los datos generaremos 5 clases, con valores entre 0.5 y 1.

maxent_xyz_normals_slope     <-
  ((subset(stack_xyz_predict_amc,6) >= 50)*1 +
   (subset(stack_xyz_predict_amc,6) <  50 & subset(stack_xyz_predict_amc,6) >= 40 )*0.9 +
   (subset(stack_xyz_predict_amc,6) <  40 & subset(stack_xyz_predict_amc,6) >= 30 )*0.8 +
   (subset(stack_xyz_predict_amc,6) <  30 & subset(stack_xyz_predict_amc,6) >= 20 )*0.7 +
   (subset(stack_xyz_predict_amc,6) <  20 & subset(stack_xyz_predict_amc,6) >= 10 )*0.6 +
   (subset(stack_xyz_predict_amc,6) <  10)*0.5)
 writeRaster(maxent_xyz_normals_slope     ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"     ),overwrite=TRUE)

Veamos si la variable slope quedó normalizada.

maxent_xyz_normals_slope@data@min
## [1] 0.5
maxent_xyz_normals_slope@data@max
## [1] 1
maxent_xyz_normals_slope <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"))
png(file.path(ruta_datos,"figs/plot_maxent_xyz_normals_slope.png"), units="px", width=1344, height=960, pointsize=24)
plot(maxent_xyz_normals_slope, main="Variable slope\nnormalizada")
dev.off()
## png 
##   2
Variable slope normalizada en rango 0 a 1..

Variable slope normalizada en rango 0 a 1..

maxent_xyz_reclass_slope <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_slope.tif"))
maxent_xyz_normals_slope <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"))
png(file.path(ruta_datos,"figs/hist_minmax_xyz_normals_slope.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
hist( maxent_xyz_reclass_slope, breaks=100, main = "Original")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
hist( maxent_xyz_normals_slope, breaks=100, main = "Normalizado")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
dev.off()
## png 
##   2
Variable slope original versus normalizada.

Variable slope original versus normalizada.

ASPECTO: lo normalizaremos en base a una regla cuatrinaria: Noreste = 1, Sureste = 0.75, Suroeste = 0.50, Sureste = 0.25.

maxent_xyz_normals_aspect <-
  ((subset(stack_xyz_predict_amc,1) >= 0   & subset(stack_xyz_predict_amc,1) <   90)*1.00 +
   (subset(stack_xyz_predict_amc,1) >= 90  & subset(stack_xyz_predict_amc,1) <  180)*0.75 +
   (subset(stack_xyz_predict_amc,1) >= 180 & subset(stack_xyz_predict_amc,1) <  270)*0.50 +
   (subset(stack_xyz_predict_amc,1) >= 270 & subset(stack_xyz_predict_amc,1) <= 360)*0.25)
 writeRaster(maxent_xyz_normals_aspect    ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"    ),overwrite=TRUE)

Veamos si la variable aspect quedó normalizada.

maxent_xyz_normals_aspect@data@min
## [1] 0
maxent_xyz_normals_aspect@data@max
## [1] 1
maxent_xyz_normals_aspect <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"))
png(file.path(ruta_datos,"figs/plot_maxent_xyz_normals_aspect.png"), units="px", width=1344, height=960, pointsize=24)
plot(maxent_xyz_normals_aspect, main="Variable aspect\nnormalizada")
dev.off()
## png 
##   2
Variable aspect normalizada en rango 0 a 1..

Variable aspect normalizada en rango 0 a 1..

Veamos el resultado.

maxent_xyz_reclass_aspect <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_reclass_aspect.tif"))
maxent_xyz_normals_aspect <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"))
png(file.path(ruta_datos,"figs/hist_minmax_xyz_normals_aspect.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
hist( maxent_xyz_reclass_aspect, breaks=100, main = "Original")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
hist( maxent_xyz_normals_aspect, breaks=100, main = "Normalizado")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2% of the raster cells were used. 100000 values used.
dev.off()
## png 
##   2
Variable aspect original versus normalizada.

Variable aspect original versus normalizada.

Volvemos a grabar en disco las capas.

 writeRaster(maxent_xyz_normals_aspect    ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"    ),overwrite=TRUE)
#writeRaster(maxent_xyz_normals_construido,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_construido.tif"),overwrite=TRUE)
 writeRaster(maxent_xyz_normals_diversidad,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_diversidad.tif"),overwrite=TRUE)
#writeRaster(maxent_xyz_normals_nativo    ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_nativo.tif"    ),overwrite=TRUE)
#writeRaster(maxent_xyz_normals_plantacion,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_plantacion.tif"),overwrite=TRUE)
 writeRaster(maxent_xyz_normals_slope     ,file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"     ),overwrite=TRUE)
archivos_normals <- list.files(file.path(ruta_datos,"raster_xyz"),pattern="normals")
as.data.frame(archivos_normals)
##                    archivos_normals
## 1     maxent_xyz_normals_aspect.tif
## 2 maxent_xyz_normals_construido.tif
## 3 maxent_xyz_normals_diversidad.tif
## 4     maxent_xyz_normals_nativo.tif
## 5 maxent_xyz_normals_plantacion.tif
## 6      maxent_xyz_normals_slope.tif
stack_xyz_normals_amc <- raster::stack(file.path(ruta_datos,"raster_xyz",archivos_normals))

4.6 Integración de indicadores mediante ecuación multicriterio.

Existen 2 formas básicas de integrar indicadores normalizados: Media Aritmética, o Media Geométrica.

Media Aritmética: (Ind.1 + Ind.2 + Ind.3 + ... + Ind.n)/n.

maxent_xyz_medarit_amc <-
  ((subset(stack_xyz_normals_amc,1))
  +(subset(stack_xyz_normals_amc,2))
  +(subset(stack_xyz_normals_amc,3))
  +(subset(stack_xyz_normals_amc,4))
  +(subset(stack_xyz_normals_amc,5))
  +(subset(stack_xyz_normals_amc,6)))/6
 writeRaster(maxent_xyz_medarit_amc     ,file.path(ruta_datos,"raster_xyz/maxent_xyz_medarit_amc.tif"     ),overwrite=TRUE)
png(file.path(ruta_datos,"figs/plot_maxent_xyz_medarit_amc.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
plot(maxent_xyz_medarit_amc, main="M aritmética", col=bpy.colors(100))
hist(maxent_xyz_medarit_amc, main="M aritmética")
dev.off()
## png 
##   2
Mapa e histograma para media aritmética.

Mapa e histograma para media aritmética.

Media Geométrica: (Ind.1 x Ind.2 x Ind.3 x ... x Ind.n)^(1/n).

maxent_xyz_medgeom_amc <-
  ((subset(stack_xyz_normals_amc,1))
  *(subset(stack_xyz_normals_amc,2))
  *(subset(stack_xyz_normals_amc,3))
  *(subset(stack_xyz_normals_amc,4))
  *(subset(stack_xyz_normals_amc,5))
  *(subset(stack_xyz_normals_amc,6)))^(1/6)
 writeRaster(maxent_xyz_medgeom_amc     ,file.path(ruta_datos,"raster_xyz/maxent_xyz_medgeom_amc.tif"     ),overwrite=TRUE)
png(file.path(ruta_datos,"figs/plot_maxent_xyz_medgeom_amc.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
plot(maxent_xyz_medgeom_amc, main="M geométrica", col=bpy.colors(100))
hist(maxent_xyz_medgeom_amc, main="M geométrica")
dev.off()
## png 
##   2
Mapa e histograma para media geométrica.

Mapa e histograma para media geométrica.

4.7 Asignar y justificar el peso relativo de los indicadores.

¿NOS FALTA ALGO? / ¿Son todas los indicadores igual de importantes? Para darle distinta importancia a cada variable podemos utilizar "Ponderadores". Definamos nuestro indicadores y ponderadores para hacer más simple cambiar los parámetros.

Idx_A <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_aspect.tif"     ))
Idx_C <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_construido.tif" ))
Idx_B <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_diversidad.tif" ))
Idx_D <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_nativo.tif"     ))
Idx_E <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_plantacion.tif" ))
Idx_F <- raster(file.path(ruta_datos,"raster_xyz/maxent_xyz_normals_slope.tif"      ))
Pa <- c(1) # aspect
Pc <- c(1) # construido
Pb <- c(1) # diversidad
Pd <- c(1) # nativo
Pe <- c(1) # plantacion
Pf <- c(1) # slope

Calculamos las medias y graficamos.

PND_AMC_MA_111111 <- (Idx_A*Pa + Idx_B*Pb + Idx_C*Pc + Idx_D*Pd + Idx_E*Pe + Idx_F*Pf)/(Pa+Pb+Pc+Pd+Pe+Pf)
PND_AMC_MG_111111 <- (Idx_A^Pa * Idx_B^Pb * Idx_C^Pc * Idx_D^Pd * Idx_E^Pe * Idx_F^Pf)^(1/(Pa+Pb+Pc+Pd+Pe+Pf))
 writeRaster(PND_AMC_MA_111111 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MA_111111.tif"     ),overwrite=TRUE)
 writeRaster(PND_AMC_MG_111111 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MG_111111.tif"     ),overwrite=TRUE)
png(file.path(ruta_datos,"figs/plot_ponderacion_indicadores_0001.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
plot(PND_AMC_MA_111111, main = "MA_111111", col=bpy.colors(100))
plot(PND_AMC_MG_111111, main = "MG_111111", col=bpy.colors(100))
dev.off()
## png 
##   2
Medias aritmética y geométrica para fórmula 111111.

Medias aritmética y geométrica para fórmula 111111.

Probemos con otros valores de ponderación. Por ejemplo:

Pa <- c(1) # aspect
Pb <- c(2) # construido
Pc <- c(1) # diversidad
Pd <- c(2) # nativo
Pe <- c(3) # plantacion
Pf <- c(3) # slope

Calculamos las medias y graficamos.

PND_AMC_MA_121233 <- (Idx_A*Pa + Idx_B*Pb + Idx_C*Pc + Idx_D*Pd + Idx_E*Pe + Idx_F*Pf)/(Pa+Pb+Pc+Pd+Pe+Pf)
PND_AMC_MG_121233 <- (Idx_A^Pa * Idx_B^Pb * Idx_C^Pc * Idx_D^Pd * Idx_E^Pe * Idx_F^Pf)^(1/(Pa+Pb+Pc+Pd+Pe+Pf))
 writeRaster(PND_AMC_MA_121233 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MA_121233.tif"     ),overwrite=TRUE)
 writeRaster(PND_AMC_MG_121233 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MG_121233.tif"     ),overwrite=TRUE)
png(file.path(ruta_datos,"figs/plot_ponderacion_indicadores_0002.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
plot(PND_AMC_MA_121233, main = "MA_121233", col=bpy.colors(100))
plot(PND_AMC_MG_121233, main = "MG_121233", col=bpy.colors(100))
dev.off()
## png 
##   2
Medias aritmética y geométrica para fórmula 121233.

Medias aritmética y geométrica para fórmula 121233.

Probemos con un último ejemplo:

Pa <- c(2) # aspect
Pb <- c(3) # construido
Pc <- c(1) # diversidad
Pd <- c(3) # nativo
Pe <- c(1) # plantacion
Pf <- c(4) # slope

Calculamos las medias y graficamos.

PND_AMC_MA_231314 <- (Idx_A*Pa + Idx_B*Pb + Idx_C*Pc + Idx_D*Pd + Idx_E*Pe + Idx_F*Pf)/(Pa+Pb+Pc+Pd+Pe+Pf)
PND_AMC_MG_231314 <- (Idx_A^Pa * Idx_B^Pb * Idx_C^Pc * Idx_D^Pd * Idx_E^Pe * Idx_F^Pf)^(1/(Pa+Pb+Pc+Pd+Pe+Pf))
 writeRaster(PND_AMC_MA_231314 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MA_231314.tif"     ),overwrite=TRUE)
 writeRaster(PND_AMC_MG_231314 ,file.path(ruta_datos,"raster_xyz/PND_AMC_MG_231314.tif"     ),overwrite=TRUE)
PND_AMC_MA_231314 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MA_231314.tif"))
PND_AMC_MG_231314 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MG_231314.tif"))
png(file.path(ruta_datos,"figs/plot_ponderacion_indicadores_0003.png"), units="px", width=1344, height=960, pointsize=24)
par(mfrow=c(1,2))
plot(PND_AMC_MA_231314, main = "MA_231314", col=bpy.colors(100))
plot(PND_AMC_MG_231314, main = "MG_231314", col=bpy.colors(100))
dev.off()
## png 
##   2
Medias aritmética y geométrica para fórmula 231314.

Medias aritmética y geométrica para fórmula 231314.

Finalmente: Pongamos nuestros rasters en una sola capa para poder mapearlos y compararlos juntos

stack_ponderacion_indicadores_0001 <-
stack(PND_AMC_MA_111111,
      PND_AMC_MA_121233,
      PND_AMC_MA_231314,
      PND_AMC_MG_111111,
      PND_AMC_MG_121233,
      PND_AMC_MG_231314)
Nombres_Capas <- c("MA_111111", "MA_121233", "MA_231314", "MG_111111", "MG_121233", "MG_231314")
names(stack_ponderacion_indicadores_0001) <- Nombres_Capas
 writeRaster(stack_ponderacion_indicadores_0001 ,file.path(ruta_datos,"raster_xyz/stack_ponderacion_indicadores_0001.tif"     ),overwrite=TRUE)
stack_ponderacion_indicadores_0001 <- raster::stack(file.path(ruta_datos,"raster_xyz/stack_ponderacion_indicadores_0001.tif"))
Nombres_Capas <- c("MA_111111", "MA_121233", "MA_231314", "MG_111111", "MG_121233", "MG_231314")
names(stack_ponderacion_indicadores_0001) <- Nombres_Capas
png(file.path(ruta_datos,"figs/stack_ponderacion_indicadores_0001.png"), units="px", width=1344, height=960, pointsize=24)
plot(stack_ponderacion_indicadores_0001, col=bpy.colors(100))
dev.off()
## png 
##   2
Resumen de medias aritméticas y geométricas.

Resumen de medias aritméticas y geométricas.

Como último paso, ahora mapearemos sólamente las áreas correspondientes al 5% de mayor prioridad. Con esta función cortamos la distribución de nuestros raster en el 5% superior.

PND_AMC_MA_111111 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MA_111111.tif"))
PND_AMC_MA_121233 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MA_121233.tif"))
PND_AMC_MA_231314 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MA_231314.tif"))
PND_AMC_MG_111111 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MG_111111.tif"))
PND_AMC_MG_121233 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MG_121233.tif"))
PND_AMC_MG_231314 <- raster(file.path(ruta_datos,"raster_xyz/PND_AMC_MG_231314.tif"))
QNT_AMC_MA_111111_Perc_5 <- quantile(PND_AMC_MA_111111, probs = c(0.95), na.rm=TRUE)
QNT_AMC_MA_121233_Perc_5 <- quantile(PND_AMC_MA_121233, probs = c(0.95), na.rm=TRUE)
QNT_AMC_MA_231314_Perc_5 <- quantile(PND_AMC_MA_231314, probs = c(0.95), na.rm=TRUE)
QNT_AMC_MG_111111_Perc_5 <- quantile(PND_AMC_MG_111111, probs = c(0.95), na.rm=TRUE)
QNT_AMC_MG_121233_Perc_5 <- quantile(PND_AMC_MG_121233, probs = c(0.95), na.rm=TRUE)
QNT_AMC_MG_231314_Perc_5 <- quantile(PND_AMC_MG_231314, probs = c(0.95), na.rm=TRUE)
funReclValPerc_5 <- function(...,strval) {
ifelse(... <= strval, 0,
ifelse(... >  strval, 1, NA)) }
RECL_AMC_MA_111111_Perc_5 <- calc(PND_AMC_MA_111111, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MA_111111_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MA_111111_Perc_5.tif"), overwrite=TRUE)
RECL_AMC_MA_121233_Perc_5 <- calc(PND_AMC_MA_121233, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MA_121233_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MA_121233_Perc_5.tif"), overwrite=TRUE)
RECL_AMC_MA_231314_Perc_5 <- calc(PND_AMC_MA_231314, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MA_231314_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MA_231314_Perc_5.tif"), overwrite=TRUE)
RECL_AMC_MG_111111_Perc_5 <- calc(PND_AMC_MG_111111, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MG_111111_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MG_111111_Perc_5.tif"), overwrite=TRUE)
RECL_AMC_MG_121233_Perc_5 <- calc(PND_AMC_MG_121233, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MG_121233_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MG_121233_Perc_5.tif"), overwrite=TRUE)
RECL_AMC_MG_231314_Perc_5 <- calc(PND_AMC_MG_231314, function(x) { funReclValPerc_5(x, strval=QNT_AMC_MG_231314_Perc_5) }, filename=file.path(ruta_datos,"raster_xyz/RECL_AMC_MG_231314_Perc_5.tif"), overwrite=TRUE)

Ahora hacemos nuevamente un raster stack y vemos los resultados.

stack_quantile_indicadores_0001 <-
stack(RECL_AMC_MA_111111_Perc_5,
      RECL_AMC_MA_121233_Perc_5,
      RECL_AMC_MA_231314_Perc_5,
      RECL_AMC_MG_111111_Perc_5,
      RECL_AMC_MG_121233_Perc_5,
      RECL_AMC_MG_231314_Perc_5)
Nombres_Capas_Perc.5 <- c("MA_111111_5%", "MA_121233_5%", "MA_231314_5%", "MG_111111_5%", "MG_121233_5%", "MG_231314_5%")
names(stack_quantile_indicadores_0001) <- Nombres_Capas_Perc.5
 writeRaster(stack_quantile_indicadores_0001 ,file.path(ruta_datos,"raster_xyz/stack_quantile_indicadores_0001.tif"     ),overwrite=TRUE)
png(file.path(ruta_datos,"figs/plot_stack_quantile_indicadores_0001.png"), units="px", width=1344, height=960, pointsize=24)
plot(stack_quantile_indicadores_0001)
dev.off()
## png 
##   2
Resumen de coberturas con mayor propensión a incendios.

Resumen de coberturas con mayor propensión a incendios.

5 Resultados y conclusiones.

En este trabajo se tomó en cuenta variables explicativas como, tipos de usos de la tierra, insolación, pendientes, NDVI y altitud sobre el nivel de mar, modeladas con técnicas de Evaluación Multicriterio (ECM) en R. Dado que el concepto riesgo se refiere a una condición de probabilidad, se podría agregar para el desarrollo de mejores escenarios la capa de eventos ocurridos de incendios que se han presentado en años anteriores, como por ejemplos los datos que se pueden obtener en The Fire Information for Resource Management System (FIRMS) https://firms.modaps.eosdis.nasa.gov/

Se debe advertir que no se consideró la variable viento como factor dispersor del fuego, ya que el objetivo era determinar lugares con condición favorable para desarrollarse un incendio sin mediar los factores que contribuyen en su distribución.

El estudio ha demostrado el potencial de riesgo de incendio en la zona de estudio. Lo anterior permite analizar la bondad de ajuste, en donde se corrobora que el modelo calculado se ajusta efectivamente a los datos usados para estimarlo.

Se puede observar la coincidencia de las áreas susceptibles determinadas por el modelo desarrollado en R y los datos “hotspots” de los últimos 30 días Enero y Febrero del 2021 en el sistema The Fire Information for Resource Management System (FIRMS)

Herramienta de diagnóstico FIRMS.

Herramienta de diagnóstico FIRMS.

El modelo presenta buenos resultados ya que existe una relación directa de incendios con las diferentes categorías de riesgo establecidas. Por ello se recomienda su uso siempre y cuando exista la información disponible para poder modelar los datos según lo establecido en la metodología.

El uso de técnicas de teledetección y Análisis espacial en R demuestra ser una herramienta muy poderosa y confiable para llevar a cabo procesos de validación y contextualizar el desempeño del modelaje de los datos, además facilita el monitoreo de áreas con susceptibilidad a presentar estos eventos.