Workflow escalonado para modelamiento de distribución de especies
Nivel regional, nivel local, comparación de modelos y ensamble final
Autor/a
Diego Ocampo Melgar
Fecha de publicación
22 de junio de 2026
1 Objetivo
Este documento presenta un workflow reproducible para modelar la distribución potencial de una especie mediante un enfoque escalonado. El análisis separa dos preguntas distintas: primero identifica la señal ambiental regional con predictores bioclimáticos de escala gruesa; luego refina la predicción en un área local usando predictores de mayor resolución espacial.
Nivel regional: estima la distribución potencial amplia usando ocurrencias de GBIF, un límite regional definido por shapefile y variables bioclimáticas. Este nivel también cumple una función de filtrado: selecciona un subconjunto de variables bioclimáticas informativas y reduce colinealidad mediante correlación de Pearson.
Nivel local: usa ocurrencias locales y predictores de mayor resolución, incluyendo capas satelitales, DEM y variables bioclimáticas seleccionadas en el Nivel 1. El objetivo no es repetir el análisis regional, sino mejorar la representación de gradientes ambientales finos dentro del área de estudio.
El flujo compara GLM, Random Forest y MaxEnt en el nivel regional. En el nivel local incorpora además CART, GAM y BRT, y genera ensambles por media, mediana, media ponderada por desempeño y committee averaging. La evaluación incluye matriz de confusión, ROC-AUC, PR-AUC, TSS, sensibilidad, especificidad y comparación visual de mapas.
Cómo leer este informe
El HTML está diseñado para lectura y auditoría. Los resultados, tablas y figuras quedan visibles; el código queda plegado por defecto y puede abrirse en cada sección con Mostrar código. La evaluación debe interpretarse junto con los mapas y las curvas de respuesta, no solo a partir de una métrica aislada.
Nota técnica
maxnet implementa una formulación tipo MaxEnt en R sin depender de Java. Si se requiere el MaxEnt clásico de Phillips et al. con maxent.jar, se puede reemplazar el bloque de maxnet por dismo::maxent().
2 Configuración
Esta sección define parámetros, carga paquetes y prepara carpetas de salida. El bloque se mantiene oculto en el HTML porque es infraestructura de ejecución, no un resultado interpretativo.
3 Funciones auxiliares
Las funciones auxiliares concentran operaciones repetidas del workflow: lectura de límites, extracción de valores raster, muestreo de background, evaluación de modelos, predicción espacial y generación de ensambles. Mantenerlas agrupadas facilita revisar el análisis sin duplicar código en cada sección.
4 Nivel 1: distribución regional
El nivel regional estima la distribución potencial amplia de la especie. Trabaja con predictores climáticos de escala gruesa y ocurrencias públicas, por lo que se usa principalmente para capturar restricciones ambientales generales y seleccionar variables bioclimáticas que luego alimentan el nivel local.
4.1 Límite regional y datos GBIF
Se carga el límite regional desde un shapefile y se descargan ocurrencias georreferenciadas desde GBIF. Los registros se depuran para mantener coordenadas válidas, remover duplicados espaciales y reducir sobre-representación dentro de una misma celda raster.
Reading layer `regiones' from data source
`G:\My Drive\Programming\R projects\ColoColo_SDM\input\shp\regiones.shp'
using driver `ESRI Shapefile'
Simple feature collection with 13 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -109.4549 ymin: -40.68094 xmax: -66.9905 ymax: -17.4984
Geodetic CRS: SIRGAS-Chile 2002
4.2 Predictores bioclimáticos gruesos
Se descargan las 19 variables bioclimáticas y se recortan al área regional. Estas variables resumen gradientes de temperatura y precipitación, por lo que son adecuadas para representar restricciones climáticas de primer orden.
4.3 Contribución de variables y reducción por colinealidad
La selección se hace en dos pasos. Primero se ajusta un Random Forest exploratorio con las 19 variables bioclimáticas para estimar importancia relativa. Luego se ordenan las variables por importancia y se retiene un subconjunto de hasta 10 predictores, evitando pares con correlación de Pearson absoluta mayor o igual al umbral definido (params$pearson_threshold). Este criterio busca reducir redundancia sin eliminar de antemano variables ecológicamente relevantes.
4.4 Entrenamiento regional de GLM, Random Forest y MaxEnt
Con el subconjunto seleccionado de variables bioclimáticas se entrenan tres algoritmos con supuestos distintos. GLM entrega una referencia paramétrica interpretable, Random Forest captura relaciones no lineales e interacciones, y MaxEnt/maxnet representa un enfoque habitual para datos de presencia-background.
4.5 Curvas ROC, precision-recall y matriz de confusión
La evaluación combina métricas de ranking y clasificación. ROC-AUC resume la separación general entre presencias y background; PR-AUC es más exigente cuando las presencias son escasas; TSS, sensibilidad y especificidad permiten revisar el balance entre omisiones y falsas presencias.
4.6 Predicción espacial regional y comparación visual
Los modelos regionales se proyectan sobre el stack ambiental regional. La comparación visual permite detectar diferencias en extensión, continuidad espacial y zonas de extrapolación potencial que no siempre son evidentes en las métricas.
5 Nivel 2: distribución local
El nivel local refina la predicción dentro de un área acotada. Usa ocurrencias locales y predictores de mayor resolución para capturar heterogeneidad ambiental fina, manteniendo coherencia con la señal bioclimática seleccionada en el nivel regional.
5.1 Límite local y ocurrencias locales
Se carga un límite local contenido dentro del área regional y una base de ocurrencias construida para el estudio. A diferencia del Nivel 1, aquí no se consulta GBIF: el objetivo es usar información local controlada y consistente con el diseño de muestreo.
5.2 Predictor de referencia Sentinel 2
Se usa un raster Landsat 8 como referencia de resolución, extensión, CRS y grilla. Todas las capas locales se alinean a esta referencia para evitar errores por diferencias de resolución, desplazamiento de celdas o extensión espacial.
Los predictores deben estar en params$sentinel_predictors_dir y pueden corresponder a bandas, índices espectrales o texturas ya calculadas.
5.3 DEM 30 m y derivados topográficos
El DEM local aporta variables topográficas relevantes para la distribución de hábitat. A partir de la elevación se derivan pendiente, rugosidad, TRI, curvatura aproximada y variables asociadas a aspecto. Estas capas representan condiciones microtopográficas que no están capturadas por las bioclimáticas gruesas.
[1] TRUE
5.4 Verdor asociado a oasis de niebla
La capa generada corresponde a un mapa de frecuencia temporal de verdor, construido con imágenes Sentinel 2 L2A entre 2017 y 2025. Representa una síntesis temporal de cuántas veces cada pixel mostró una señal espectral compatible con vegetación escasa, efímera o vegetación más persistente.
El criterio está basado en el enfoque de Moat et al. 2021, quienes usaron NDVI en composiciones de 16 días para distinguir vegetación efímera y núcleos verdantes de oasis de niebla. En el estudio original se usó MODIS a 250 m; aquí se adapta la lógica a Landsat 8, con resolución espacial de 30 m. El paper define NDVI entre 0.10 y 0.15 como vegetación escasa o efímera, y NDVI > 0.15 como vegetación herbácea/leñosa de oasis de niebla.
El procesamiento realiza los siguientes pasos:
Excluye coberturas no deseadas usando ESA WorldCover: urbano, cuerpos de agua, agricultura, nieve/hielo.
Filtra pixeles problemáticos de Landsat 8 usando las bandas de calidad.
Calcula NDVI para cada escena Landsat 8 válida.
Agrupa las imágenes en ventanas de 16 días, imitando la lógica temporal del producto MODIS usado en el paper.
Para cada ventana de 16 días, identifica si el pixel tuvo:
- NDVI entre 0.10 y 0.15 → vegetación escasa/efímera.
- NDVI mayor a 0.15 → vegetación verdante más clara.
Cuenta cuántas veces ocurrió cada condición durante todo el periodo.
Convierte esa frecuencia en clases finales:
0: Sin clasificación (No mostró señal suficiente o quedó fuera de la máscara)
1: Vegetación efímera (Señal débil y muy ocasional, NDVI 0.10–0.15)
2: Oasis verdante muy ocasional (NDVI > 0.15 entre 2–4% de los periodos válidos)
3: Oasis verdante ocasional (NDVI > 0.15 entre 4–8% de los periodos válidos)
4: Oasis verdante estacional/intermitente (NDVI > 0.15 entre 8–25% de los periodos válidos)
5: Oasis verdante persistente (NDVI > 0.15 entre 25–50% de los periodos válidos)
6: Oasis verdante muy persistente (NDVI > 0.15 en más del 50% de los periodos válidos)
Posteriormente, el raster fue sometido a una limpieza espacial para reducir pixeles aislados y ruido tipo “sal y pimienta”. Las bandas categóricas fueron suavizadas mediante una ventana focal de moda, asignando a cada pixel la clase dominante en su vecindad local. Las bandas continuas fueron suavizadas mediante mediana focal, reduciendo valores extremos locales sin modificar fuertemente la estructura espacial general. Se mantuvieron bandas auxiliares para identificar los pixeles rellenados o modificados durante el suavizado.
[1] TRUE
5.5 Frecuencia de niebla y nubosidad
La capa de frecuencia de niebla costera fue generada como un proxy satelital de recurrencia de nubosidad baja o niebla potencial, usando imágenes Landsat 8 Collection 2 Level 2 entre 2013 y 2025. El análisis fue restringido a una franja costera dentro del área de estudio y excluyó áreas urbanas, agrícolas, cuerpos de agua y nieve/hielo. Para cada imagen se identificaron pixeles con señal tipo nube/fog mediante las banderas de calidad QA_PIXEL de Landsat 8, combinadas con criterios espectrales de brillo y blancura en las bandas visibles. Posteriormente, la serie fue agregada en ventanas de 16 días y se calculó el porcentaje de periodos válidos en que cada pixel presentó señal compatible con niebla o nube baja.
La banda principal anual es: all_fog_pct
Representa el porcentaje de periodos válidos de 16 días en que cada pixel tuvo señal compatible con niebla/nube baja costera.
La banda estacional es: season_fog_pct
Representa lo mismo, pero solo para los meses definidos como temporada de niebla. Por defecto mayo a noviembre.
Posteriormente, el raster fue sometido a una limpieza espacial para reducir pixeles aislados y ruido tipo “sal y pimienta”. Las bandas categóricas fueron suavizadas mediante una ventana focal de moda, asignando a cada pixel la clase dominante en su vecindad local. Las bandas continuas fueron suavizadas mediante mediana focal, reduciendo valores extremos locales sin modificar fuertemente la estructura espacial general. Se mantuvieron bandas auxiliares para identificar los pixeles rellenados o modificados durante el suavizado.
[1] TRUE
5.6 Frecuencia de humedad
Se generó una capa complementaria de humedad superficial y agua superficial potencial a partir de imágenes Landsat 8 Collection 2 Level 2, manteniendo la misma lógica temporal aplicada a las capas de verdor y nubosidad/niebla. El análisis se restringió a la franja costera del área de estudio y se aplicaron máscaras de calidad para excluir nubes, sombras, nieve, saturación radiométrica y coberturas no objetivo como urbano y agricultura. Para agua superficial se utilizaron índices espectrales orientados a realzar cuerpos de agua, principalmente MNDWI y AWEI, debido a su capacidad de reducir confusiones con suelo, vegetación, sombras y superficies construidas. Para humedad superficial se utilizó NDMI/LSWI, basado en la relación entre NIR y SWIR, sensible a contenido de agua en vegetación y superficies húmedas.
La serie temporal fue agregada en ventanas de 16 días, siguiendo la lógica de frecuencia temporal utilizada para mapear oasis de niebla. Para cada pixel se calculó la recurrencia de condiciones compatibles con agua superficial y la frecuencia de señales de humedad superficial relativa. El resultado debe interpretarse como un proxy óptico de humedad superficial y agua expuesta, no como una medición directa de humedad volumétrica del suelo. Por esta razón, la capa se utiliza como variable auxiliar para interpretar la relación entre niebla costera, respuesta vegetal y disponibilidad superficial de humedad.
Banda
Interpretación
surface_water_pct
Porcentaje de periodos válidos con agua superficial detectada
surface_water_count
Número de periodos de 16 días con agua superficial
water_valid_count
Número de periodos válidos para agua
surface_water_class
Clase ordinal de recurrencia de agua
surface_wetness_freq
Porcentaje de periodos con NDMI alto relativo
surface_wetness_count
Número de periodos con señal alta de humedad superficial
5.9 Variables bioclimáticas seleccionadas del Nivel 1 remuestreadas a Sentinel-2
Las variables bioclimáticas seleccionadas en el Nivel 1 se incorporan al Nivel 2 después de ser remuestreadas a la grilla de referencia. Este paso preserva la señal climática regional, pero evita reinterpretar la resolución original de las bioclimáticas como si fuera información climática local fina.
[1] TRUE
5.10 Ensamble de predictores locales
Aquí se integran los predictores locales en un único stack ambiental: capas satelitales, topográficas, proxies de humedad/niebla/verdor y variables bioclimáticas seleccionadas. Luego se remueven capas sin variación o con problemas de nombres para asegurar compatibilidad con los modelos.
El nivel local compara seis modelos: GLM, RF, MaxEnt, CART, GAM y BRT. La evaluación se realiza sobre datos de prueba separados del entrenamiento; si se activa el aumento de presencias por buffer, esos puntos adicionales se usan solo para entrenamiento y no para evaluar desempeño.
5.12 Predicción espacial local y comparación visual
Los modelos locales se proyectan sobre la grilla ambiental de alta resolución. La comparación visual es parte de la evaluación: un modelo puede tener buenas métricas y aun así producir artefactos espaciales, bordes abruptos o patrones ecológicamente poco plausibles.
5.13 Ensamble local de modelos
El ensamble se calcula usando todos los modelos disponibles en models_local: GLM, RF, MaxEnt, CART, GAM y BRT. Se generan cuatro salidas complementarias: media de probabilidades, mediana de probabilidades, media ponderada por TSS y committee averaging de predicciones binarias. La media y la mediana resumen consenso probabilístico; la media ponderada favorece modelos con mejor desempeño; el committee averaging expresa la proporción de modelos que clasifican cada celda como presencia.
6 Comparación global y decisión de modelo
Esta sección resume el desempeño de los modelos regionales y locales, identifica el mejor modelo individual y contrasta esa decisión con los ensambles. La selección final no debe basarse únicamente en la métrica máxima, sino en el equilibrio entre desempeño, estabilidad espacial e interpretabilidad ecológica.
La decisión de modelo se basa en tres criterios:
Desempeño predictivo: prioridad a modelos con AUC-ROC, AUC-PR y TSS altos en datos de prueba.
Balance de errores: revisión explícita de sensibilidad y especificidad para distinguir entre omisiones y falsas presencias.
Plausibilidad ecológica: revisión de mapas, curvas parciales y superficies de respuesta para descartar artefactos, extrapolaciones no justificadas o respuestas incompatibles con el conocimiento de la especie.
El ensamble no reemplaza automáticamente al mejor modelo individual. Debe interpretarse como una predicción compuesta que puede ser más estable cuando los modelos coinciden, pero también puede ocultar desacuerdos importantes entre algoritmos.
[1] "RF"
En Nivel 1 regional, el set de prueba tiene 26 presencias y 299 negativos aleatorios. La prevalencia es ~8%. Por eso la precisión parece baja incluso en buenos modelos.
Esto no significa necesariamente que el modelo sea malo. Significa que, al haber muchos negativos, incluso una cantidad moderada de falsos positivos baja mucho la precisión.
En Nivel 2 local, el problema es más severo: solo hay 11 presencias y 191 negativos. La prevalencia es ~5.4%. El AUC-PR base esperado por azar está cerca de 0.054.
Todos superan el azar, pero no por mucho. Esto indica que el modelo local detecta alguna estructura ambiental, pero la señal es débil o el test tiene pocas presencias.
Para el modelo regional: Random Forest. Tiene el mejor AUC_ROC = 0.9098, mejor AUC_PR = 0.5652, alta especificidad 0.8528, buena sensibilidad 0.8462 y menos falsos positivos que MaxEnt. MaxEnt tiene un TSS levemente mayor, pero RF es más equilibrado.
Para el modelo local: MaxEnt. Tiene el mejor TSS = 0.5669, sensibilidad alta 0.8182 y especificidad aceptable 0.7487. Es el mejor compromiso para presencia-background/pseudoausencia, especialmente si el objetivo es delimitar hábitat potencial.
7 Curvas de respuesta parcial e inflated response curves
Las curvas de respuesta parcial muestran la predicción esperada al variar un predictor y mantener los demás constantes. Las inflated response curves extienden esta lógica incorporando el rango observado de los datos de calibración, lo que permite visualizar zonas donde el modelo extrapola o depende de supuestos fuera del dominio ambiental observado.
7.1 Función para curvas parciales e inflated response curves
7.2 Curvas para predictores principales
7.3 Superficies de respuesta 3D para dos predictores
Estas superficies no son mapas espaciales. Representan la predicción en el espacio ambiental definido por dos variables continuas, manteniendo el resto de predictores en su mediana. La altura corresponde a la probabilidad predicha.
8 Exportación final
Los productos finales se exportan como rasters para uso posterior en SIG, validación de campo o priorización espacial. Se guarda el mapa del mejor modelo local y el ensamble ponderado por desempeño, que suele ser una alternativa más estable cuando varios algoritmos entregan señales complementarias.
Barbosa, A. M. modEvA: paquete para evaluación y análisis de modelos de distribución de especies.
Zurell, D. et al. Protocolo ODMAP para reporte de modelos de distribución de especies.
Zurell, D. SDM Intro: respuesta parcial, superficies de respuesta e inflated response curves.
Zurell, D. EEC-MGC: ensemble species distribution modelling, model averaging, weighted means and committee averaging.
Phillips, S. J. et al. MaxEnt para modelamiento de distribución con datos de presencia y background.
Elith, J., Leathwick, J. R., & Hastie, T. Boosted regression trees for ecological modelling.
Wood, S. N. Generalized additive models: an introduction with R.
Ejecutar el código
---title: "Workflow escalonado para modelamiento de distribución de especies"subtitle: "Nivel regional, nivel local, comparación de modelos y ensamble final"author: "Diego Ocampo Melgar"date: todaylang: esformat: html: theme: light: flatly dark: darkly css: styles.css toc: true toc-title: "Contenido" toc-depth: 3 toc-location: right number-sections: true anchor-sections: true smooth-scroll: true page-layout: article title-block-banner: true code-fold: true code-summary: "Mostrar código" code-tools: true code-copy: true code-overflow: wrap fig-align: center fig-width: 8 fig-height: 5 fig-cap-location: bottom tbl-cap-location: top df-print: paged link-external-newwindow: trueexecute: echo: true warning: false message: false freeze: autoparams: species: "Leopardus colocola" out_dir: "output_sdm" seed: 123 crs_work: "EPSG:4326" # Nivel 1: distribución regional regional_boundary: "limit0.shp" gbif_country: "CL" gbif_limit: 5000 bg_n_regional: 5000 bio_res_arcmin: 0.5 pearson_threshold: 0.70 n_selected_bioclim: 10 train_prop: 0.70 # Nivel 2: distribución local local_boundary: "limit1.shp" local_occurrences: "puntos.csv" # columnas: lon, lat, presence opcional sentinel_reference: "S2_surface_wetness_class_2017_2025_20m.tif" verdant_reference: "S2_fog_oasis_class_only_AOI_coastal_2017_2025_10m.tif" fog_reference: "L8_fog_frequency_AOI_coastal_2013_2025_30m_filled_light.tif" wetness_reference: "S2_surface_wetness_class_2017_2025_20m.tif" sentinel_predictors_dir: "S2_summer_mean_stack_2025_full.tif" dem_provider: "archivo_local" # geodata o archivo_local local_dem_file: "NASADEM_30m_Chinchilla.tif" bg_n_local: 2000 # Aumento opcional de presencias para entrenamiento. # Se usa solo en entrenamiento; la evaluación permanece con puntos originales. augment_presence: true augment_buffer_m: 250 augment_n_per_presence: 3 augment_crs_metric: "EPSG:32719" augment_apply_regional: false augment_apply_local: true---# ObjetivoEste documento presenta un workflow reproducible para modelar la distribución potencial de una especie mediante un enfoque escalonado. El análisis separa dos preguntas distintas: primero identifica la señal ambiental regional con predictores bioclimáticos de escala gruesa; luego refina la predicción en un área local usando predictores de mayor resolución espacial.1. **Nivel regional:** estima la distribución potencial amplia usando ocurrencias de GBIF, un límite regional definido por shapefile y variables bioclimáticas. Este nivel también cumple una función de filtrado: selecciona un subconjunto de variables bioclimáticas informativas y reduce colinealidad mediante correlación de Pearson.2. **Nivel local:** usa ocurrencias locales y predictores de mayor resolución, incluyendo capas satelitales, DEM y variables bioclimáticas seleccionadas en el Nivel 1. El objetivo no es repetir el análisis regional, sino mejorar la representación de gradientes ambientales finos dentro del área de estudio.El flujo compara **GLM**, **Random Forest** y **MaxEnt** en el nivel regional. En el nivel local incorpora además **CART**, **GAM** y **BRT**, y genera ensambles por media, mediana, media ponderada por desempeño y committee averaging. La evaluación incluye matriz de confusión, ROC-AUC, PR-AUC, TSS, sensibilidad, especificidad y comparación visual de mapas.::: callout-note## Cómo leer este informeEl HTML está diseñado para lectura y auditoría. Los resultados, tablas y figuras quedan visibles; el código queda plegado por defecto y puede abrirse en cada sección con **Mostrar código**. La evaluación debe interpretarse junto con los mapas y las curvas de respuesta, no solo a partir de una métrica aislada.:::::: callout-important## Nota técnica`maxnet` implementa una formulación tipo MaxEnt en R sin depender de Java. Si se requiere el MaxEnt clásico de Phillips et al. con `maxent.jar`, se puede reemplazar el bloque de `maxnet` por `dismo::maxent()`.:::# ConfiguraciónEsta sección define parámetros, carga paquetes y prepara carpetas de salida. El bloque se mantiene oculto en el HTML porque es infraestructura de ejecución, no un resultado interpretativo.```{r setup}#| warning: false#| include: falseset.seed(params$seed)pkgs <-c("sf", "terra", "dplyr", "tidyr", "purrr", "ggplot2", "stringr", "readr","rgbif", "geodata", "randomForest", "maxnet", "pROC", "PRROC","modEvA", "yardstick", "patchwork", "vip", "pdp", "here", "dismo","mgcv", "rpart", "gbm", "lattice", "RColorBrewer", "kernlab")missing_pkgs <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly =TRUE)]if (length(missing_pkgs) >0) {stop("Faltan paquetes: ", paste(missing_pkgs, collapse =", "),"\nInstálelos antes de ejecutar este Quarto. Para modEvA use: ","install.packages('modEvA', repos = 'http://R-Forge.R-project.org')" )}invisible(lapply(pkgs, library, character.only =TRUE))dir.create(params$out_dir, showWarnings =FALSE, recursive =TRUE)dir.create(file.path(params$out_dir, "figures"), showWarnings =FALSE, recursive =TRUE)dir.create(file.path(params$out_dir, "rasters"), showWarnings =FALSE, recursive =TRUE)dir.create(file.path(params$out_dir, "tables"), showWarnings =FALSE, recursive =TRUE)species_clean <- stringr::str_replace_all(params$species, "[^A-Za-z0-9]+", "_")```# Funciones auxiliaresLas funciones auxiliares concentran operaciones repetidas del workflow: lectura de límites, extracción de valores raster, muestreo de background, evaluación de modelos, predicción espacial y generación de ensambles. Mantenerlas agrupadas facilita revisar el análisis sin duplicar código en cada sección.```{r helper-functions}#| warning: false#| include: falseutils_env <-new.env()sys.source("utilsFixed.r", envir = utils_env)loaded_functions <-ls(utils_env)[vapply(mget(ls(utils_env), envir = utils_env), is.function,logical(1) )]source("utilsFixed.r")source("utilsNoCART.r")knitr::kable(data.frame(function_name = loaded_functions),caption ="Functions loaded from utils.r")`%||%`<-function(a, b) if (!is.null(a)) a else b```# Nivel 1: distribución regionalEl nivel regional estima la distribución potencial amplia de la especie. Trabaja con predictores climáticos de escala gruesa y ocurrencias públicas, por lo que se usa principalmente para capturar restricciones ambientales generales y seleccionar variables bioclimáticas que luego alimentan el nivel local.## Límite regional y datos GBIFSe carga el límite regional desde un shapefile y se descargan ocurrencias georreferenciadas desde GBIF. Los registros se depuran para mantener coordenadas válidas, remover duplicados espaciales y reducir sobre-representación dentro de una misma celda raster.```{r level1-boundary-gbif}#| echo: falseregional_boundary <-read_boundary(here("input", "shp", params$regional_boundary), params$crs_work)bb <- sf::st_bbox(regional_boundary)# Consulta GBIF. Ajuste 'hasCoordinate = TRUE' y 'hasGeospatialIssue = FALSE' para reducir errores.gbif_raw <- rgbif::occ_search(scientificName = params$species,country = params$gbif_country,hasCoordinate =TRUE,hasGeospatialIssue =FALSE,limit = params$gbif_limit)$dataocc_regional <- gbif_raw |> dplyr::transmute(lon = decimalLongitude,lat = decimalLatitude,basisOfRecord = basisOfRecord,year = year,source ="GBIF" ) |>extract_xy()occ_regional_sf <- sf::st_as_sf(occ_regional, coords =c("lon", "lat"), crs =4326, remove =FALSE) |> sf::st_filter(regional_boundary, .predicate = sf::st_intersects)occ_regional <- occ_regional_sf |> sf::st_drop_geometry()readr::write_csv(occ_regional, file.path(params$out_dir, "tables", paste0(species_clean, "_gbif_regional_clean.csv")))regiones<-st_read(here("input", "shp", "regiones.shp"))regiones <-st_transform(regiones, st_crs(regional_boundary))regiones<-st_crop(regiones, regional_boundary)plot(sf::st_geometry(regional_boundary), main ="Ocurrencias regionales GBIF dentro del área de estudio")plot(sf::st_geometry(regiones), add =TRUE, border ="blue")plot(sf::st_geometry(occ_regional_sf), add =TRUE, pch =20, col ="red")```## Predictores bioclimáticos gruesosSe descargan las 19 variables bioclimáticas y se recortan al área regional. Estas variables resumen gradientes de temperatura y precipitación, por lo que son adecuadas para representar restricciones climáticas de primer orden.```{r level1-bioclim}#| echo: falsebio_global <- geodata::worldclim_global(var ="bio",res = params$bio_res_arcmin,path =file.path(params$out_dir, "geodata"))bio_regional <-crop_mask_stack(bio_global, regional_boundary)names(bio_regional) <-paste0("bio", 1:19)# ==========================================================# Add topographic predictors at bioclimatic resolution# ==========================================================# Use the regional bioclimatic raster as the exact grid templatebio_template <- bio_regional[[1]]# Ensure boundary is a terra vectorregional_boundary_v <- regional_boundaryif (inherits(regional_boundary_v, "sf")) { regional_boundary_v <- terra::vect(regional_boundary_v)}# Match boundary CRS to bioclimatic predictorsregional_boundary_bio_crs <- terra::project( regional_boundary_v, terra::crs(bio_template))# Download global elevation at the same nominal resolution as WorldClim bioclimdem_global <- geodata::elevation_global(res = params$bio_res_arcmin,path =file.path(params$out_dir, "geodata"))# Crop first to reduce processingdem_crop <- terra::crop( dem_global, terra::ext(bio_template),snap ="out")# Force DEM to match exactly the bioclimatic griddem_regional <- terra::project( dem_crop, bio_template,method ="bilinear")# Mask to regional boundarydem_regional <- terra::mask( dem_regional, regional_boundary_bio_crs)names(dem_regional) <-"elevation"# Calculate slope and aspect from aligned DEMterrain_vars <- terra::terrain( dem_regional,v =c("slope", "aspect"),unit ="degrees",neighbors =8)names(terrain_vars) <-c("slope_deg", "aspect_deg")# Convert aspect to radiansaspect_rad <- terrain_vars$aspect_deg * pi /180# Northness = cos(aspect)northness <- terra::ifel(is.na(terrain_vars$slope_deg),NA, terra::ifel(is.na(terrain_vars$aspect_deg) | terrain_vars$slope_deg <=0.001,0,cos(aspect_rad) ))# Eastness = sin(aspect)eastness <- terra::ifel(is.na(terrain_vars$slope_deg),NA, terra::ifel(is.na(terrain_vars$aspect_deg) | terrain_vars$slope_deg <=0.001,0,sin(aspect_rad) ))names(northness) <-"northness"names(eastness) <-"eastness"# Final topographic stacktopo_regional <-c( dem_regional, terrain_vars$slope_deg, northness, eastness)# Check geometry before combiningterra::compareGeom( bio_regional, topo_regional,stopOnError =TRUE)# Add topographic predictors to bioclimatic stackpredictors_regional <-c( bio_regional, topo_regional)names(predictors_regional)#predictors_regionalocc_regional <-thin_by_raster_cell(occ_regional, predictors_regional)occ_regional$ID <-seq(1:nrow(occ_regional))bg_regional <-sample_background(predictors_regional, params$bg_n_regional)bg_regional$ID <-seq(1:nrow(bg_regional))pa_regional_all <-make_pa_table(occ_regional, bg_regional, predictors_regional)spl_regional <-split_train_test(pa_regional_all, params$train_prop)train_regional_all <- spl_regional$traintest_regional_all <- spl_regional$testtrain_regional_all <-add_augmented_presences_to_train(train_dat = train_regional_all,boundary_sf = regional_boundary,predictors = predictors_regional,augment_presence =isTRUE(params$augment_presence) &&isTRUE(params$augment_apply_regional),buffer_m = params$augment_buffer_m,n_per_presence = params$augment_n_per_presence,metric_crs = params$augment_crs_metric,source_crs = params$crs_work)bio_vars <-names(predictors_regional)```## Contribución de variables y reducción por colinealidadLa selección se hace en dos pasos. Primero se ajusta un Random Forest exploratorio con las 19 variables bioclimáticas para estimar importancia relativa. Luego se ordenan las variables por importancia y se retiene un subconjunto de hasta 10 predictores, evitando pares con correlación de Pearson absoluta mayor o igual al umbral definido (`params$pearson_threshold`). Este criterio busca reducir redundancia sin eliminar de antemano variables ecológicamente relevantes.```{r level1-variable-selection}#| echo: falserf_screen <- randomForest::randomForest(factor(presence) ~ .,data = train_regional_all |> dplyr::select(presence, dplyr::all_of(bio_vars)),ntree =500,importance =TRUE)imp <- randomForest::importance(rf_screen, type =1)imp_tbl <- tibble::tibble(variable =rownames(imp),importance =as.numeric(imp[, 1])) |> dplyr::arrange(dplyr::desc(importance))selected_bio <-cor_select(dat = train_regional_all,vars = bio_vars,threshold = params$pearson_threshold,priority = imp_tbl$variable) |>head(params$n_selected_bioclim)cor_matrix <- stats::cor(train_regional_all[, selected_bio], use ="pairwise.complete.obs", method ="pearson")readr::write_csv(imp_tbl, file.path(params$out_dir, "tables", "level1_bioclim_variable_importance_rf_screen.csv"))readr::write_csv(tibble::tibble(selected_bio = selected_bio), file.path(params$out_dir, "tables", "level1_selected_bioclim.csv"))imp_tbl |> dplyr::slice_head(n =19)selected_bioround(cor_matrix, 2)vip::vip(rf_screen, num_features =19) +theme_bw() +labs(title ="Contribución inicial de variables bioclimáticas")```## Entrenamiento regional de GLM, Random Forest y MaxEntCon el subconjunto seleccionado de variables bioclimáticas se entrenan tres algoritmos con supuestos distintos. GLM entrega una referencia paramétrica interpretable, Random Forest captura relaciones no lineales e interacciones, y MaxEnt/maxnet representa un enfoque habitual para datos de presencia-background.```{r level1-fit-models}#| echo: false#| warning: falsemodels_regional <-fit_models(train_regional_all, selected_bio, extra_models =FALSE)pred_regional_test <- purrr::imap_dfr(models_regional, function(mod, nm) { tibble::tibble(level ="Nivel 1 regional",model = nm,presence = test_regional_all$presence,pred =predict_model_df(mod, test_regional_all, nm, selected_bio) )})eval_regional <- pred_regional_test |> dplyr::group_by(model) |> tidyr::nest() |> dplyr::mutate(eval = purrr::map2( data, model,~calc_eval(.x$presence, .x$pred, .y, "Nivel 1 regional") |> dplyr::select(-model) ) ) |> dplyr::select(-data) |> tidyr::unnest(eval)```## Curvas ROC, precision-recall y matriz de confusiónLa evaluación combina métricas de ranking y clasificación. ROC-AUC resume la separación general entre presencias y background; PR-AUC es más exigente cuando las presencias son escasas; TSS, sensibilidad y especificidad permiten revisar el balance entre omisiones y falsas presencias.```{r level1-evaluation-plots, fig.height=5, fig.width=12}#| echo: false#| warning: falseplots_regional <-plot_roc_pr( eval_regional, pred_regional_test,"Nivel 1 regional")plots_regional$rocplots_regional$pr``````{r level1-confusion-matrices}#| echo: false#| warning: falseconf_regional <- pred_regional_test |> dplyr::left_join(eval_regional |> dplyr::select(model, threshold_youden), by ="model") |> dplyr::mutate(pred_class =ifelse(pred >= threshold_youden, 1L, 0L)) |> dplyr::count(model, presence, pred_class, name ="n")conf_regional```## Predicción espacial regional y comparación visualLos modelos regionales se proyectan sobre el stack ambiental regional. La comparación visual permite detectar diferencias en extensión, continuidad espacial y zonas de extrapolación potencial que no siempre son evidentes en las métricas.```{r level1-spatial-prediction, fig.height=5, fig.width=12}#| echo: false#| warning: falsepred_rasters_regional <- purrr::imap(models_regional, function(mod, nm) { r <-predict_model_raster(mod, predictors_regional, nm, selected_bio)names(r) <- nm terra::writeRaster(r, file.path(params$out_dir, "rasters", paste0("level1_regional_", nm, ".tif")), overwrite =TRUE) r})pred_stack_regional <- terra::rast(pred_rasters_regional)plot_prediction_maps(pred_stack_regional, "Predicción regional")terra::writeRaster( pred_stack_regional,file.path(params$out_dir, "rasters", paste0("FINAL_level1.tif")),overwrite =TRUE)```# Nivel 2: distribución localEl nivel local refina la predicción dentro de un área acotada. Usa ocurrencias locales y predictores de mayor resolución para capturar heterogeneidad ambiental fina, manteniendo coherencia con la señal bioclimática seleccionada en el nivel regional.## Límite local y ocurrencias localesSe carga un límite local contenido dentro del área regional y una base de ocurrencias construida para el estudio. A diferencia del Nivel 1, aquí no se consulta GBIF: el objetivo es usar información local controlada y consistente con el diseño de muestreo.```{r level2-boundary-occurrences}#| echo: false#| warning: falselocal_boundary <-read_boundary(here("input", "shp", params$local_boundary), params$crs_work)puntos<-here("input", "table", params$local_occurrences)if (!file.exists(puntos)) stop("No existe la base local de ocurrencias: ", params$local_occurrences)occ_local <- readr::read_csv(puntos, show_col_types =FALSE) |>extract_xy()occ_local_sf <- sf::st_as_sf(occ_local, coords =c("lon", "lat"), crs =4326, remove =FALSE) |> sf::st_filter(local_boundary, .predicate = sf::st_intersects)occ_local <- occ_local_sf |> sf::st_drop_geometry()#regiones <- st_transform(regiones, st_crs(local_boundary))regiones<-st_crop(regiones, local_boundary)plot(sf::st_geometry(local_boundary), main ="Ocurrencias locales dentro del área local")plot(sf::st_geometry(regiones), add =TRUE, border ="blue")plot(sf::st_geometry(occ_local_sf), add =TRUE, pch =20, col ="red")```## Predictor de referencia Sentinel 2Se usa un raster Landsat 8 como referencia de resolución, extensión, CRS y grilla. Todas las capas locales se alinean a esta referencia para evitar errores por diferencias de resolución, desplazamiento de celdas o extensión espacial.Los predictores deben estar en `params$sentinel_predictors_dir` y pueden corresponder a bandas, índices espectrales o texturas ya calculadas.```{r level2-sentinel}#| echo: false#| warning: falseif (!file.exists(here("input","ras","V0", params$sentinel_reference))) stop("No existe el raster de referencia: ", params$sentinel_reference)s2_ref <- terra::rast(here("input","ras","V0", params$sentinel_reference))```## DEM 30 m y derivados topográficosEl DEM local aporta variables topográficas relevantes para la distribución de hábitat. A partir de la elevación se derivan pendiente, rugosidad, TRI, curvatura aproximada y variables asociadas a aspecto. Estas capas representan condiciones microtopográficas que no están capturadas por las bioclimáticas gruesas.```{r level2-dem-terrain}#| echo: false#| warning: falsedem <- terra::rast(here("input","ras", params$local_dem_file))# Ensure boundary is a terra vectorlocal_boundary_v <- local_boundaryif (inherits(local_boundary_v, "sf")) { local_boundary_v <- terra::vect(local_boundary_v)}# Match boundary CRS to bioclimatic predictorslocal_boundary_bio_crs <- terra::project( local_boundary_v, terra::crs(s2_ref))# Force DEM to match exactly the bioclimatic griddem_local <- terra::project( dem, s2_ref,method ="bilinear")# Mask to regional boundarydem_local <- terra::mask( dem_local, local_boundary_bio_crs)# Calculate slope and aspect from aligned DEMterrain_vars <- terra::terrain( dem_local[["elevation"]],v =c("aspect"),unit ="degrees",neighbors =8)names(terrain_vars) <-c("aspect_deg")# Convert aspect to radiansaspect_rad <- terrain_vars$aspect_deg * pi /180terrain_vars$slope_deg<-dem_local$slope_deg# Northness = cos(aspect)northness <- terra::ifel(is.na(terrain_vars$slope_deg),NA, terra::ifel(is.na(terrain_vars$aspect_deg) | terrain_vars$slope_deg <=0.001,0,cos(aspect_rad) ))# Eastness = sin(aspect)eastness <- terra::ifel(is.na(terrain_vars$slope_deg),NA, terra::ifel(is.na(terrain_vars$aspect_deg) | terrain_vars$slope_deg <=0.001,0,sin(aspect_rad) ))names(northness) <-"northness"names(eastness) <-"eastness"# Final topographic stacktopo_local <-c( dem_local,#terrain_vars$slope_deg, northness, eastness)# Check geometry before combiningterra::compareGeom( s2_ref, topo_local,stopOnError =TRUE)roughness <- terra::terrain(topo_local[["elevation"]], v ="roughness")slope_rad <- topo_local[["slope_deg"]] * pi /180twi_approx <-log((roughness +1) /tan(slope_rad +0.001))names(roughness) <-"roughness"names(twi_approx) <-"twi_approx"topo_pred <-c(topo_local, roughness, twi_approx)```## Verdor asociado a oasis de nieblaLa capa generada corresponde a un **mapa de frecuencia temporal de verdor**, construido con imágenes Sentinel 2 L2A entre 2017 y 2025. Representa una síntesis temporal de cuántas veces cada pixel mostró una señal espectral compatible con vegetación escasa, efímera o vegetación más persistente.El criterio está basado en el enfoque de Moat et al. 2021, quienes usaron NDVI en composiciones de 16 días para distinguir vegetación efímera y núcleos verdantes de oasis de niebla. En el estudio original se usó MODIS a 250 m; aquí se adapta la lógica a Landsat 8, con resolución espacial de 30 m. El paper define NDVI entre 0.10 y 0.15 como vegetación escasa o efímera, y NDVI \> 0.15 como vegetación herbácea/leñosa de oasis de niebla.El procesamiento realiza los siguientes pasos:1. **Excluye coberturas no deseadas** usando ESA WorldCover: urbano, cuerpos de agua, agricultura, nieve/hielo.2. **Filtra pixeles problemáticos de Landsat 8** usando las bandas de calidad.3. Calcula **NDVI** para cada escena Landsat 8 válida.4. Agrupa las imágenes en ventanas de **16 días**, imitando la lógica temporal del producto MODIS usado en el paper.5. Para cada ventana de 16 días, identifica si el pixel tuvo:``` - NDVI entre 0.10 y 0.15 → vegetación escasa/efímera.- NDVI mayor a 0.15 → vegetación verdante más clara.```6. Cuenta cuántas veces ocurrió cada condición durante todo el periodo.7. Convierte esa frecuencia en clases finales:- 0: Sin clasificación (No mostró señal suficiente o quedó fuera de la máscara)- 1: Vegetación efímera (Señal débil y muy ocasional, NDVI 0.10–0.15)- 2: Oasis verdante muy ocasional (NDVI \> 0.15 entre 2–4% de los periodos válidos)- 3: Oasis verdante ocasional (NDVI \> 0.15 entre 4–8% de los periodos válidos)- 4: Oasis verdante estacional/intermitente (NDVI \> 0.15 entre 8–25% de los periodos válidos)- 5: Oasis verdante persistente (NDVI \> 0.15 entre 25–50% de los periodos válidos)- 6: Oasis verdante muy persistente (NDVI \> 0.15 en más del 50% de los periodos válidos)Posteriormente, el raster fue sometido a una limpieza espacial para reducir pixeles aislados y ruido tipo “sal y pimienta”. Las bandas categóricas fueron suavizadas mediante una ventana focal de moda, asignando a cada pixel la clase dominante en su vecindad local. Las bandas continuas fueron suavizadas mediante mediana focal, reduciendo valores extremos locales sin modificar fuertemente la estructura espacial general. Se mantuvieron bandas auxiliares para identificar los pixeles rellenados o modificados durante el suavizado.```{r}#| echo: false#| warning: falseif (!file.exists(here("input","ras","V0", params$verdant_reference))) stop("No existe el raster de referencia: ", params$verdant_reference)s2_ref <- terra::rast(here("input","ras","V0", params$verdant_reference))#s2_ref <- crop_mask_stack(s2_ref, local_boundary)s2_ref<-s2_ref[["fog_oasis_class_final"]]#s2_ref <- project(s2_ref, "EPSG:4326", method = "bilinear")# Force DEM to match exactly the bioclimatic grids2_ref <- terra::project( s2_ref, topo_pred,method ="bilinear")# Mask to regional boundarys2_ref <- terra::mask( s2_ref, local_boundary_bio_crs)# Check geometry before combiningterra::compareGeom( s2_ref, topo_pred,stopOnError =TRUE)predictors <-c(topo_pred, s2_ref)```## Frecuencia de niebla y nubosidadLa capa de frecuencia de niebla costera fue generada como un proxy satelital de recurrencia de nubosidad baja o niebla potencial, usando imágenes Landsat 8 Collection 2 Level 2 entre 2013 y 2025. El análisis fue restringido a una franja costera dentro del área de estudio y excluyó áreas urbanas, agrícolas, cuerpos de agua y nieve/hielo. Para cada imagen se identificaron pixeles con señal tipo nube/fog mediante las banderas de calidad `QA_PIXEL` de Landsat 8, combinadas con criterios espectrales de brillo y blancura en las bandas visibles. Posteriormente, la serie fue agregada en ventanas de 16 días y se calculó el porcentaje de periodos válidos en que cada pixel presentó señal compatible con niebla o nube baja.- La banda principal anual es: all_fog_pct Representa el **porcentaje de periodos válidos de 16 días** en que cada pixel tuvo señal compatible con niebla/nube baja costera. La banda estacional es: season_fog_pct Representa lo mismo, pero solo para los meses definidos como temporada de niebla. Por defecto mayo a noviembre.Posteriormente, el raster fue sometido a una limpieza espacial para reducir pixeles aislados y ruido tipo “sal y pimienta”. Las bandas categóricas fueron suavizadas mediante una ventana focal de moda, asignando a cada pixel la clase dominante en su vecindad local. Las bandas continuas fueron suavizadas mediante mediana focal, reduciendo valores extremos locales sin modificar fuertemente la estructura espacial general. Se mantuvieron bandas auxiliares para identificar los pixeles rellenados o modificados durante el suavizado.```{r}#| label: fog_raster#| echo: false#| warning: falseif (!file.exists(here("input","ras", params$fog_reference))) stop("No existe el raster de referencia: ", params$fog_reference)s2_ref <- terra::rast(here("input","ras", params$fog_reference))#s2_ref <- crop_mask_stack(s2_ref, local_boundary)s2_ref<-s2_ref[["all_fog_pct"]]# Force DEM to match exactly the bioclimatic grids2_ref <- terra::project( s2_ref, topo_pred,method ="bilinear")# Mask to regional boundarys2_ref <- terra::mask( s2_ref, local_boundary_bio_crs)# Check geometry before combiningterra::compareGeom( s2_ref, topo_pred,stopOnError =TRUE)predictors <-c(predictors, s2_ref)```## Frecuencia de humedadSe generó una capa complementaria de humedad superficial y agua superficial potencial a partir de imágenes Landsat 8 Collection 2 Level 2, manteniendo la misma lógica temporal aplicada a las capas de verdor y nubosidad/niebla. El análisis se restringió a la franja costera del área de estudio y se aplicaron máscaras de calidad para excluir nubes, sombras, nieve, saturación radiométrica y coberturas no objetivo como urbano y agricultura. Para agua superficial se utilizaron índices espectrales orientados a realzar cuerpos de agua, principalmente MNDWI y AWEI, debido a su capacidad de reducir confusiones con suelo, vegetación, sombras y superficies construidas. Para humedad superficial se utilizó NDMI/LSWI, basado en la relación entre NIR y SWIR, sensible a contenido de agua en vegetación y superficies húmedas.La serie temporal fue agregada en ventanas de 16 días, siguiendo la lógica de frecuencia temporal utilizada para mapear oasis de niebla. Para cada pixel se calculó la recurrencia de condiciones compatibles con agua superficial y la frecuencia de señales de humedad superficial relativa. El resultado debe interpretarse como un proxy óptico de humedad superficial y agua expuesta, no como una medición directa de humedad volumétrica del suelo. Por esta razón, la capa se utiliza como variable auxiliar para interpretar la relación entre niebla costera, respuesta vegetal y disponibilidad superficial de humedad.| Banda | Interpretación ||----|----||`surface_water_pct`| Porcentaje de periodos válidos con agua superficial detectada ||`surface_water_count`| Número de periodos de 16 días con agua superficial ||`water_valid_count`| Número de periodos válidos para agua ||`surface_water_class`| Clase ordinal de recurrencia de agua ||`surface_wetness_freq`| Porcentaje de periodos con NDMI alto relativo ||`surface_wetness_count`| Número de periodos con señal alta de humedad superficial ||`wetness_valid_count`| Número de periodos válidos para humedad ||`surface_wetness_class`| Clase ordinal de humedad superficial relativa ||`surface_wetness_mean`| NDMI medio temporal ||`surface_wetness_p75`| Percentil 75 temporal de NDMI ||`surface_wetness_p90`| Percentil 90 temporal de NDMI ||`NDMI_median`| NDMI mediano del periodo completo ||`elevation`| Elevación |```{r}#| label: wetness_raster#| echo: false#| warning: falseif (!file.exists(here("input","ras","V0", params$wetness_reference))) stop("No existe el raster de referencia: ", params$wetness_reference)s2_ref <- terra::rast(here("input","ras","V0", params$wetness_reference))s2_ref<-s2_ref[["surface_wetness_class"]]# Force DEM to match exactly the bioclimatic grids2_ref <- terra::project( s2_ref, topo_pred,method ="bilinear")# Mask to regional boundarys2_ref <- terra::mask( s2_ref, local_boundary_bio_crs)# Check geometry before combiningterra::compareGeom( s2_ref, topo_pred,stopOnError =TRUE)predictors <-c(predictors, s2_ref)```## Sentinel 2 Indices for summer 2025```{r}#| label: indices_raster#| echo: false#| warning: falseif (!file.exists(here("input","ras","V0", params$sentinel_predictors_dir))) stop("No existe el raster de referencia: ", params$sentinel_predictors_dir)s2_ref <- terra::rast(here("input","ras","V0", params$sentinel_predictors_dir))s2_selected <-c("SAVI", # vegetación ajustada por suelo"NDMI", # humedad / estrés hídrico"NBR2", # señal SWIR independiente"NDRE"# condición de vegetación red-edge)s2_ref<-s2_ref[[s2_selected]]#terra::pairs(s2_ref[[s2_selected]])#cor_df <- terra::as.data.frame(# s2_ref[[s2_selected]],# na.rm = TRUE#)#round(cor(cor_df, use = "complete.obs"), 2)# Force DEM to match exactly the bioclimatic grids2_ref <- terra::project( s2_ref, topo_pred,method ="bilinear")# Mask to regional boundarys2_ref <- terra::mask( s2_ref, local_boundary_bio_crs)# Check geometry before combiningterra::compareGeom( s2_ref, topo_pred,stopOnError =TRUE)predictors <-c(predictors, s2_ref)```## Predictores de proximidad```{r}#| label: distance_raster#| echo: false#| warning: falses2_ref <- terra::rast(here("input","ras","proximity","costa_prox_0_1.tif"))s2_ref <-project(s2_ref, "EPSG:4326", method ="bilinear")s2_ref <-crop_mask_stack(s2_ref, local_boundary)ras1<-s2_refs2_ref <- terra::rast(here("input","ras","proximity","quebradas_prox_0_1.tif"))s2_ref <-project(s2_ref, "EPSG:4326", method ="bilinear")s2_ref <-crop_mask_stack(s2_ref, local_boundary)ras1<-c(ras1, s2_ref)s2_ref <- terra::rast(here("input","ras","proximity","rios_prox_0_1.tif"))s2_ref <-project(s2_ref, "EPSG:4326", method ="bilinear")s2_ref <-crop_mask_stack(s2_ref, local_boundary)ras1<-c(ras1, s2_ref)# Template: misma grilla que predictorstemplate <- predictors[[1]]# Reproyectar + remuestrear ras1 a la grilla de predictors# Como son variables continuas 0-1, usa bilinearras1_aligned <- terra::project( ras1, template,method ="bilinear")# Opcional: asegurar rango 0-1 después de interpolaciónras1_aligned <- terra::clamp( ras1_aligned,lower =0,upper =1,values =TRUE)# Verificar geometríaterra::compareGeom( predictors, ras1_aligned,stopOnError =FALSE)# Unir al stackpredictors <-c(predictors, ras1_aligned)names(predictors)```## Variables bioclimáticas seleccionadas del Nivel 1 remuestreadas a Sentinel-2Las variables bioclimáticas seleccionadas en el Nivel 1 se incorporan al Nivel 2 después de ser remuestreadas a la grilla de referencia. Este paso preserva la señal climática regional, pero evita reinterpretar la resolución original de las bioclimáticas como si fuera información climática local fina.```{r level2-bioclim-resample}#| echo: false#| warning: falsetags<-gsub(pattern ="_",replacement ="",x =substr(names(bio_global), start =11, stop =20))names(bio_global)<-tagssel_bio<-selected_bio[which(selected_bio %in% tags)]bio_local_selected <- bio_global[[sel_bio]]bio_local_selected <- terra::project(bio_local_selected, predictors, method ="bilinear") |> terra::resample(predictors, method ="bilinear") |>crop_mask_stack(local_boundary)names(bio_local_selected) <- sel_bio# Force DEM to match exactly the bioclimatic gridbio_local_selected <- terra::project( bio_local_selected, topo_pred,method ="bilinear")# Mask to regional boundarybio_local_selected <- terra::mask( bio_local_selected, local_boundary_bio_crs)# Check geometry before combiningterra::compareGeom( bio_local_selected, topo_pred,stopOnError =TRUE)local_predictors_all <-c(predictors, bio_local_selected)local_predictors_all <-crop_mask_stack(local_predictors_all, local_boundary)```## Ensamble de predictores localesAquí se integran los predictores locales en un único stack ambiental: capas satelitales, topográficas, proxies de humedad/niebla/verdor y variables bioclimáticas seleccionadas. Luego se remueven capas sin variación o con problemas de nombres para asegurar compatibilidad con los modelos.```{r level2-predictor-stack}#| echo: false#| warning: false# Remueve capas sin variación o con nombres problemáticos.local_sd <- terra::global(local_predictors_all, "sd", na.rm =TRUE)[, 1]local_predictors_all <- local_predictors_all[[which(is.finite(local_sd) & local_sd >0)]]names(local_predictors_all) <-make.names(names(local_predictors_all), unique =TRUE)occ_local <-thin_by_raster_cell(occ_local, local_predictors_all)bg_local <-sample_background(local_predictors_all, params$bg_n_local)pa_local_all <-make_pa_table(occ_local, bg_local, local_predictors_all)spl_local <-split_train_test(pa_local_all, params$train_prop)train_local <- spl_local$traintest_local <- spl_local$testtrain_local <-add_augmented_presences_to_train(train_dat = train_local,boundary_sf = local_boundary,predictors = local_predictors_all,augment_presence =isTRUE(params$augment_presence) &&isTRUE(params$augment_apply_local),buffer_m = params$augment_buffer_m,n_per_presence = params$augment_n_per_presence,metric_crs = params$augment_crs_metric,source_crs = params$crs_work)local_vars_all <-names(local_predictors_all)# Selección local: RF exploratorio + filtro Pearson. Permite más predictores que nivel 1, pero controla redundancia.# Si augment_presence = TRUE, la selección usa las presencias aumentadas sólo en entrenamiento.rf_local_screen <- randomForest::randomForest(factor(presence) ~ .,data = train_local |> dplyr::select(presence, dplyr::all_of(local_vars_all)),ntree =500,importance =TRUE)local_imp <- randomForest::importance(rf_local_screen, type =1)local_imp_tbl <- tibble::tibble(variable =rownames(local_imp), importance =as.numeric(local_imp[, 1])) |> dplyr::arrange(dplyr::desc(importance))selected_local <-cor_select(dat = train_local,vars = local_vars_all,threshold = params$pearson_threshold,priority = local_imp_tbl$variable) |>head(15)readr::write_csv(local_imp_tbl, file.path(params$out_dir, "tables", "level2_variable_importance_rf_screen.csv"))readr::write_csv(tibble::tibble(selected_local = selected_local), file.path(params$out_dir, "tables", "level2_selected_predictors.csv"))selected_localvip::vip(rf_local_screen, num_features =min(25, length(local_vars_all))) +theme_bw() +labs(title ="Contribución inicial de predictores locales")```## Entrenamiento y evaluación localEl nivel local compara seis modelos: GLM, RF, MaxEnt, CART, GAM y BRT. La evaluación se realiza sobre datos de prueba separados del entrenamiento; si se activa el aumento de presencias por buffer, esos puntos adicionales se usan solo para entrenamiento y no para evaluar desempeño.```{r level2-fit-evaluate}#| echo: false#| warning: falsemodels_local <-fit_models(train_local, selected_local, extra_models =TRUE)# Verificación: el nivel local debe incluir seis modelos.names(models_local)pred_local_test <- purrr::imap_dfr(models_local, function(mod, nm) { tibble::tibble(level ="Nivel 2 local",model = nm,presence = test_local$presence,pred =predict_model_df(mod, test_local, nm, selected_local) )})pred_local_test |> dplyr::group_by(model) |> dplyr::summarise(n = dplyr::n(),n_presence =sum(presence ==1, na.rm =TRUE),n_absence =sum(presence ==0, na.rm =TRUE),pred_min =min(pred, na.rm =TRUE),pred_max =max(pred, na.rm =TRUE),pred_mean =mean(pred, na.rm =TRUE),pred_na =sum(is.na(pred)),.groups ="drop" )eval_local <- pred_local_test |> dplyr::filter(!is.na(presence),!is.na(pred) ) |> (\(x) split(x, x$model))() |> purrr::imap_dfr(function(df, nm) {calc_eval(obs = df$presence,pred = df$pred,model_name = nm,level_name ="Nivel 2 local" ) })readr::write_csv(eval_local, file.path(params$out_dir, "tables", "level2_model_evaluation.csv"))eval_local``````{r level2-evaluation-plots, fig.height=5, fig.width=12}#| echo: false#| warning: falseplot_roc_pr(eval_local, pred_local_test, "Nivel 2 local")``````{r level2-confusion-matrices}#| echo: false#| warning: false#| conf_local <- pred_local_test |> dplyr::left_join(eval_local |> dplyr::select(model, threshold_youden), by ="model") |> dplyr::mutate(pred_class =ifelse(pred >= threshold_youden, 1L, 0L)) |> dplyr::count(model, presence, pred_class, name ="n")conf_local```## Predicción espacial local y comparación visualLos modelos locales se proyectan sobre la grilla ambiental de alta resolución. La comparación visual es parte de la evaluación: un modelo puede tener buenas métricas y aun así producir artefactos espaciales, bordes abruptos o patrones ecológicamente poco plausibles.```{r level2-spatial-prediction, fig.height=5, fig.width=12}#| echo: false#| warning: falsepred_rasters_local <- purrr::imap(models_local, function(mod, nm) { r <-predict_model_raster(mod, local_predictors_all, nm, selected_local)names(r) <- nm terra::writeRaster(r, file.path(params$out_dir, "rasters", paste0("level2_local_", nm, ".tif")), overwrite =TRUE) r})pred_stack_local <- terra::rast(pred_rasters_local)plot_prediction_maps(pred_stack_local, "Predicción local")```## Ensamble local de modelosEl ensamble se calcula usando todos los modelos disponibles en `models_local`: GLM, RF, MaxEnt, CART, GAM y BRT. Se generan cuatro salidas complementarias: media de probabilidades, mediana de probabilidades, media ponderada por TSS y committee averaging de predicciones binarias. La media y la mediana resumen consenso probabilístico; la media ponderada favorece modelos con mejor desempeño; el committee averaging expresa la proporción de modelos que clasifican cada celda como presencia.```{r level2-ensemble-prediction, fig.height=5, fig.width=12}#| echo: false#| warning: falsemodel_names_local <-names(pred_stack_local)threshold_vec <- eval_local |> dplyr::filter(model %in% model_names_local) |> dplyr::select(model, threshold_youden) |> tibble::deframe()threshold_vec <- threshold_vec[model_names_local]threshold_vec[!is.finite(threshold_vec)] <-0.5weight_vec <- eval_local |> dplyr::filter(model %in% model_names_local) |> dplyr::select(model, TSS) |> tibble::deframe()weight_vec <- weight_vec[model_names_local]weight_vec[!is.finite(weight_vec)] <-0weight_vec[weight_vec <0] <-0if (sum(weight_vec) ==0) { weight_vec[] <-1}ensemble_mean <- terra::app( pred_stack_local,fun =function(x) mean(x, na.rm =TRUE))names(ensemble_mean) <-"ensemble_mean_probability"ensemble_median <- terra::app( pred_stack_local,fun =function(x) stats::median(x, na.rm =TRUE))names(ensemble_median) <-"ensemble_median_probability"ensemble_weighted <- terra::app( pred_stack_local,fun =function(x) stats::weighted.mean(x, w = weight_vec, na.rm =TRUE))names(ensemble_weighted) <-"ensemble_weighted_mean_TSS"binary_layers <- purrr::map(model_names_local, function(nm) { r <- pred_stack_local[[nm]] >= threshold_vec[[nm]]names(r) <- nm r}) |> terra::rast()ensemble_committee <- terra::app( binary_layers,fun =function(x) mean(x, na.rm =TRUE))names(ensemble_committee) <-"ensemble_committee_averaging"ensemble_sd <- terra::app( pred_stack_local,fun =function(x) stats::sd(x, na.rm =TRUE))names(ensemble_sd) <-"ensemble_model_disagreement_sd"ensemble_stack_local <-c( ensemble_mean, ensemble_median, ensemble_weighted, ensemble_committee, ensemble_sd)terra::writeRaster( ensemble_stack_local,file.path(params$out_dir, "rasters", "level2_local_ensemble_stack.tif"),overwrite =TRUE)terra::writeRaster( ensemble_weighted,file.path(params$out_dir, "rasters", "FINAL_level2_weighted_ensemble_TSS.tif"),overwrite =TRUE)plot_prediction_maps(ensemble_stack_local, "Ensambles locales")```# Comparación global y decisión de modeloEsta sección resume el desempeño de los modelos regionales y locales, identifica el mejor modelo individual y contrasta esa decisión con los ensambles. La selección final no debe basarse únicamente en la métrica máxima, sino en el equilibrio entre desempeño, estabilidad espacial e interpretabilidad ecológica.La decisión de modelo se basa en tres criterios:1. **Desempeño predictivo:** prioridad a modelos con AUC-ROC, AUC-PR y TSS altos en datos de prueba.2. **Balance de errores:** revisión explícita de sensibilidad y especificidad para distinguir entre omisiones y falsas presencias.3. **Plausibilidad ecológica:** revisión de mapas, curvas parciales y superficies de respuesta para descartar artefactos, extrapolaciones no justificadas o respuestas incompatibles con el conocimiento de la especie.El ensamble no reemplaza automáticamente al mejor modelo individual. Debe interpretarse como una predicción compuesta que puede ser más estable cuando los modelos coinciden, pero también puede ocultar desacuerdos importantes entre algoritmos.```{r model-decision}#| echo: false#| warning: falseeval_all <- dplyr::bind_rows(eval_regional, eval_local) |> dplyr::mutate(score_rank = dplyr::min_rank(dplyr::desc(AUC_PR)) + dplyr::min_rank(dplyr::desc(TSS)) + dplyr::min_rank(dplyr::desc(AUC_ROC)) ) |> dplyr::arrange(level, score_rank)readr::write_csv(eval_all, file.path(params$out_dir, "tables", "model_comparison_all_levels.csv"))eval_allbest_local_name <- eval_local |> dplyr::arrange(dplyr::desc(AUC_PR), dplyr::desc(TSS), dplyr::desc(AUC_ROC)) |> dplyr::slice(1) |> dplyr::pull(model)best_local_model <- models_local[[best_local_name]]best_local_name``````{r model-comparison-plot, fig.height=5, fig.width=9}#| echo: false#| warning: false#| eval_all |> dplyr::select(level, model, AUC_ROC, AUC_PR, TSS, sensitivity, specificity) |> tidyr::pivot_longer(cols =c(AUC_ROC, AUC_PR, TSS, sensitivity, specificity), names_to ="metric", values_to ="value") |>ggplot(aes(model, value, fill = model)) +geom_col(show.legend =FALSE) +facet_grid(level ~ metric) +coord_cartesian(ylim =c(0, 1)) +theme_bw() +labs(title ="Comparación de desempeño por nivel y modelo", x =NULL, y ="Valor")```En Nivel 1 regional, el set de prueba tiene `26` presencias y `299` negativos aleatorios. La prevalencia es \~`8%`. Por eso la precisión parece baja incluso en buenos modelos.Esto no significa necesariamente que el modelo sea malo. Significa que, al haber muchos negativos, incluso una cantidad moderada de falsos positivos baja mucho la precisión.En Nivel 2 local, el problema es más severo: solo hay `11` presencias y `191` negativos. La prevalencia es \~`5.4%`. El AUC-PR base esperado por azar está cerca de `0.054`.Todos superan el azar, pero **no por mucho**. Esto indica que el modelo local detecta alguna estructura ambiental, pero la señal es débil o el test tiene pocas presencias.- Para el modelo regional: Random Forest. Tiene el mejor `AUC_ROC = 0.9098`, mejor `AUC_PR = 0.5652`, alta especificidad `0.8528`, buena sensibilidad `0.8462` y menos falsos positivos que MaxEnt. MaxEnt tiene un TSS levemente mayor, pero RF es más equilibrado.- Para el modelo local: MaxEnt. Tiene el mejor `TSS = 0.5669`, sensibilidad alta `0.8182` y especificidad aceptable `0.7487`. Es el mejor compromiso para presencia-background/pseudoausencia, especialmente si el objetivo es delimitar hábitat potencial.# Curvas de respuesta parcial e inflated response curvesLas curvas de respuesta parcial muestran la predicción esperada al variar un predictor y mantener los demás constantes. Las **inflated response curves** extienden esta lógica incorporando el rango observado de los datos de calibración, lo que permite visualizar zonas donde el modelo extrapola o depende de supuestos fuera del dominio ambiental observado.## Función para curvas parciales e inflated response curves```{r response-functions}#| echo: false#| warning: falsepredict_best_local <-function(object, newdata) {predict_model_df(object, newdata, best_local_name, selected_local)}partial_response_data <-function(model, train_dat, vars, focal_var, n =100) { base <- train_dat[, vars, drop =FALSE] med <- base |> dplyr::summarise(dplyr::across(dplyr::everything(), ~ stats::median(.x, na.rm =TRUE))) grid <- med[rep(1, n), , drop =FALSE] grid[[focal_var]] <-seq(min(base[[focal_var]], na.rm =TRUE), max(base[[focal_var]], na.rm =TRUE), length.out = n) grid$prediction <-predict_best_local(model, grid) grid$focal_var <- focal_var grid$x <- grid[[focal_var]] grid$type <-"partial_median_context" grid}inflated_response_data <-function(model, train_dat, vars, focal_var, n =100, sample_n =500) { base <- train_dat[, vars, drop =FALSE] sampled <- base |> dplyr::slice_sample(n =min(sample_n, nrow(base))) xs <-seq(min(base[[focal_var]], na.rm =TRUE), max(base[[focal_var]], na.rm =TRUE), length.out = n) out <- purrr::map_dfr(xs, function(xv) { nd <- sampled nd[[focal_var]] <- xv preds <-predict_best_local(model, nd) tibble::tibble(x = xv,prediction = preds,focal_var = focal_var ) }) |> dplyr::group_by(focal_var, x) |> dplyr::summarise(pred_mean =mean(prediction, na.rm =TRUE),pred_lwr = stats::quantile(prediction, 0.025, na.rm =TRUE),pred_upr = stats::quantile(prediction, 0.975, na.rm =TRUE),.groups ="drop" ) out}plot_inflated_response <-function(model, train_dat, vars, focal_var) { pr <-partial_response_data(model, train_dat, vars, focal_var) ir <-inflated_response_data(model, train_dat, vars, focal_var)ggplot() +geom_ribbon(data = ir, aes(x = x, ymin = pred_lwr, ymax = pred_upr), alpha =0.25) +geom_line(data = ir, aes(x = x, y = pred_mean), linewidth =0.8) +geom_line(data = pr, aes(x = x, y = prediction), linetype ="dashed", linewidth =0.7) +geom_rug(data = train_dat, aes(x = .data[[focal_var]], color =factor(presence)), sides ="b", alpha =0.4) +theme_bw() +labs(title =paste("Inflated response curve:", focal_var),subtitle =paste("Modelo seleccionado:", best_local_name),x = focal_var,y ="Probabilidad / adecuabilidad predicha",color ="Presencia" )}```## Curvas para predictores principales```{r response-curves, fig.height=7, fig.width=10}#| echo: false#| warning: falsetop_response_vars <- local_imp_tbl |> dplyr::filter(variable %in% selected_local) |> dplyr::slice_head(n =min(6, length(selected_local))) |> dplyr::pull(variable)response_plots <- purrr::map(top_response_vars, ~plot_inflated_response(best_local_model, train_local, selected_local, .x))patchwork::wrap_plots(response_plots, ncol =2)```## Superficies de respuesta 3D para dos predictoresEstas superficies no son mapas espaciales. Representan la predicción en el espacio ambiental definido por dos variables continuas, manteniendo el resto de predictores en su mediana. La altura corresponde a la probabilidad predicha.```{r response-surface-3d, fig.height=7, fig.width=9}#| echo: false#| warning: falsetags <-as.character(tags)bio_tags <- tags[ stringr::str_detect(tags, "^bio[0-9]+$")]bio_surface_pool <- selected_local[ selected_local %in% bio_tags]is_valid_surface_var <-function(v) { v %in%names(train_local) &&is.numeric(train_local[[v]]) && dplyr::n_distinct(train_local[[v]], na.rm =TRUE) >=10}# First priority: top response variables, but only if they are bioclimatic tagssurface_candidate_vars <- top_response_vars[ top_response_vars %in% bio_surface_pool]surface_candidate_vars <- surface_candidate_vars[vapply( surface_candidate_vars, is_valid_surface_var,logical(1) )]# Fallback: any selected local bioclimatic variable present in tagsif (length(surface_candidate_vars) <2) { surface_candidate_vars <- bio_surface_pool[vapply( bio_surface_pool, is_valid_surface_var,logical(1) ) ]}if (length(surface_candidate_vars) <2) {stop("No hay al menos dos predictores bioclimáticos válidos en `tags` ","para construir la superficie 3D. Revisa `tags`, `selected_local` ","y que las variables existan en `train_local`." )}if (length(surface_candidate_vars) >=2) { v1 <- surface_candidate_vars[1] v2 <- surface_candidate_vars[2] base_env <- train_local[, selected_local, drop =FALSE] |> dplyr::summarise( dplyr::across( dplyr::everything(),~ stats::median(.x, na.rm =TRUE) ) ) surface_grid <-expand.grid(x1 =seq(min(train_local[[v1]], na.rm =TRUE),max(train_local[[v1]], na.rm =TRUE),length.out =60 ),x2 =seq(min(train_local[[v2]], na.rm =TRUE),max(train_local[[v2]], na.rm =TRUE),length.out =60 ) ) nd <- base_env[rep(1, nrow(surface_grid)), , drop =FALSE] nd[[v1]] <- surface_grid$x1 nd[[v2]] <- surface_grid$x2 cls <- grDevices::colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdYlBu")) )(100) plot_surface_3d <-function(df, z_col, title, z_label ="Occurrence probability") { plot_df <- df plot_df$z <- plot_df[[z_col]] lattice::wireframe( z ~ x1 + x2,data = plot_df,zlab =list(z_label, rot =90),xlab = v1,ylab = v2,main = title,drape =TRUE,col.regions = cls,scales =list(arrows =FALSE),zlim =c(0, 1),screen =list(z =-120, x =-70, y =3) ) }# Best local model response surface surface_grid$best_model <-predict_best_local(best_local_model, nd) p_best_3d <-plot_surface_3d( surface_grid,"best_model",paste("Response surface:", best_local_name) )print(p_best_3d) predict_surface_model <-function(model, model_name, nd, selected_vars) { model_vars <-attr(model, "predictor_vars") %||% selected_vars model_vars <- model_vars[model_vars %in%names(nd)]predict_model_df(model = model,newdata = nd,model_name = model_name,vars = model_vars ) } surface_model_preds <- purrr::imap_dfc( models_local,function(mod, nm) { tibble::tibble(!!nm :=predict_surface_model(model = mod,model_name = nm,nd = nd,selected_vars = selected_local ) ) } ) |> dplyr::mutate( dplyr::across( dplyr::everything(),~pmin(pmax(as.numeric(.x), 0), 1) ) ) model_names_surface <-names(surface_model_preds)# ==========================================================# Align thresholds and weights to surface model names# ==========================================================model_names_surface <-names(surface_model_preds)threshold_surface <-rep(0.5, length(model_names_surface))names(threshold_surface) <- model_names_surfacecommon_threshold_models <-intersect( model_names_surface,names(threshold_vec))threshold_surface[common_threshold_models] <- threshold_vec[common_threshold_models]threshold_surface[!is.finite(threshold_surface)] <-0.5weight_surface <-rep(1, length(model_names_surface))names(weight_surface) <- model_names_surfacecommon_weight_models <-intersect( model_names_surface,names(weight_vec))weight_surface[common_weight_models] <- weight_vec[common_weight_models]weight_surface[!is.finite(weight_surface)] <-0weight_surface[weight_surface <0] <-0if (sum(weight_surface, na.rm =TRUE) ==0) { weight_surface[] <-1}threshold_vecnames(threshold_vec)weight_vecnames(weight_vec)names(surface_model_preds)surface_bin_preds <- surface_model_predsfor (nm in model_names_surface) { surface_bin_preds[[nm]] <-ifelse( surface_model_preds[[nm]] >= threshold_surface[nm],1,0 )} surface_plot_df <- dplyr::bind_cols( surface_grid[, c("x1", "x2")], tibble::tibble(mean_prob =rowMeans(surface_model_preds, na.rm =TRUE),median_prob =apply(surface_model_preds, 1, stats::median, na.rm =TRUE),weighted_mean_prob =apply( surface_model_preds,1, stats::weighted.mean,w = weight_surface,na.rm =TRUE ),committee_averaging =rowMeans(surface_bin_preds, na.rm =TRUE),sd_prob =apply(surface_model_preds, 1, stats::sd, na.rm =TRUE) ) )print(plot_surface_3d(surface_plot_df, "mean_prob", "Ensemble response surface: mean probability"))print(plot_surface_3d(surface_plot_df, "median_prob", "Ensemble response surface: median probability"))print(plot_surface_3d(surface_plot_df, "weighted_mean_prob", "Ensemble response surface: weighted mean probability"))print(plot_surface_3d(surface_plot_df, "committee_averaging", "Ensemble response surface: committee averaging", "Proportion of models predicting presence"))print(plot_surface_3d(surface_plot_df, "sd_prob", "Ensemble response surface: model disagreement", "SD of predicted probability"))}```# Exportación finalLos productos finales se exportan como rasters para uso posterior en SIG, validación de campo o priorización espacial. Se guarda el mapa del mejor modelo local y el ensamble ponderado por desempeño, que suele ser una alternativa más estable cuando varios algoritmos entregan señales complementarias.```{r export-final}#| echo: false#| warning: falsefinal_local_prediction <- pred_stack_local[[best_local_name]]terra::writeRaster( final_local_prediction,file.path(params$out_dir, "rasters", paste0("FINAL_level2_selected_model_", best_local_name, ".tif")),overwrite =TRUE)# Ensamble recomendado: media ponderada por TSS.terra::writeRaster( ensemble_weighted,file.path(params$out_dir, "rasters", "FINAL_level2_ensemble_weighted_TSS.tif"),overwrite =TRUE)list(selected_model_local = best_local_name,models_local =names(models_local),selected_bioclim_level1 = selected_bio,selected_predictors_level2 = selected_local,final_selected_model_raster =file.path(params$out_dir, "rasters", paste0("FINAL_level2_selected_model_", best_local_name, ".tif")),final_weighted_ensemble_raster =file.path(params$out_dir, "rasters", "FINAL_level2_ensemble_weighted_TSS.tif"),output_directory =normalizePath(params$out_dir))```# Referencias operativas- Barbosa, A. M. `modEvA`: paquete para evaluación y análisis de modelos de distribución de especies.- Zurell, D. et al. Protocolo ODMAP para reporte de modelos de distribución de especies.- Zurell, D. SDM Intro: respuesta parcial, superficies de respuesta e inflated response curves.- Zurell, D. EEC-MGC: ensemble species distribution modelling, model averaging, weighted means and committee averaging.- Phillips, S. J. et al. MaxEnt para modelamiento de distribución con datos de presencia y background.- Elith, J., Leathwick, J. R., & Hastie, T. Boosted regression trees for ecological modelling.- Wood, S. N. Generalized additive models: an introduction with R.```{r}```
