WSIM-GLDAS: Adquisición, Exploración e Integración

Authors

Josh Brinks

Elaine Famutimi

Published

July 30, 2024

Descripción General

En esta lección, adquirirás el conjunto de datos llamado Modelo de Indicadores de Seguridad Hídrica - Sistema Global de Asimilación de Datos Territoriales (WSIM-GLDAS) desde el sitio web del NASA Socioeconomic Data and Applications Center (SEDAC), realizarás varias tareas de preprocesamiento y explorarás el conjunto de datos con visualizaciones avanzadas. También integraremos WSIM-GLDAS con un conjunto de datos de límites administrativos globales llamado geoBoundaries y un conjunto de datos de población llamado Gridded Population of the World. Aprenderás cómo realizar numerosas manipulaciones de datos y visualizaciones a lo largo de esta lección.

Objetivos de Aprendizaje

Después de completar esta lección, deberías ser capaz de:

  • Recuperar datos de WSIM-GLDAS desde SEDAC.
  • Recuperar límites administrativos utilizando la API de geoBoundaries.
  • Subconjuntar datos de WSIM-GLDAS para una región y período de interés.
  • Visualizar datos geoespaciales para resaltar patrones de déficit de precipitación.
  • Escribir un archivo preprocesado en formato NetCDF al disco.
  • Realizar exploración visual con histogramas, coropletas y mapas de series temporales.
  • Integrar datos de población enrejada con datos de WSIM-GLDAS para realizar análisis y construir visualizaciones para entender cómo las personas son impactadas.
  • Resumir datos raster de WSIM-GLDAS y población utilizando estadísticas zonales.

Introducción

El ciclo del agua es el proceso de circulación del agua en, sobre y bajo la superficie de la Tierra (NOAA 2019). Las actividades humanas producen emisiones de gases de efecto invernadero, cambios en el uso del suelo, desarrollo de presas y embalses, y extracción de aguas subterráneas que han afectado el ciclo natural del agua en las últimas décadas (IPCC 2023). La influencia de estas actividades humanas en el ciclo del agua tiene impactos consecuentes en los procesos oceánicos, de aguas subterráneas y terrestres, influyendo en fenómenos como sequías e inundaciones (Zhou 2016).

Los déficits de precipitación, o períodos de lluvia por debajo del promedio, pueden llevar a la sequía, caracterizada por períodos prolongados de poca o ninguna lluvia y escasez de agua resultante. Las sequías a menudo desencadenan estrés ambiental y pueden crear ciclos de refuerzo que impactan ecosistemas y personas (Rodgers 2023). Por ejemplo, California frecuentemente experimenta sequías, pero la combinación de períodos secos prolongados y temperaturas sostenidas altas impidió el reabastecimiento de agua fresca y fría al río Klamath, lo que llevó a graves escaseces de agua en 2003 y nuevamente entre 2012 y 2014. Estas escaseces afectan áreas agrícolas como el Valle Central, que cultiva almendras, uno de los cultivos más importantes de California, con el estado produciendo el 80% de las almendras del mundo. Estas sequías severas, junto con la competencia por recursos limitados de agua dulce, resultaron en poblaciones decrecientes de salmón Chinook debido al estrés térmico y la enfermedad de putrefacción de las branquias, lo que interrumpió el suministro de alimentos para los grupos tribales de la cuenca de Klamath (Guillen 2002; Bland 2014).

Ciclo del Agua1

Revisión de Ciencia de Datos

Un raster es un tipo de datos geográficos en formato de imagen digital con información numérica almacenada en cada píxel. (Los rasteres a menudo se llaman cuadrículas debido a su estructura de datos matricial regularmente formada). Los rasteres pueden almacenar muchos tipos de información y pueden tener dimensiones que incluyen latitud, longitud y tiempo. NetCDF es un formato para datos raster; otros incluyen Geotiff, ASCII y muchos más. Varios formatos raster como NetCDF pueden almacenar múltiples capas raster, o una “pila raster”, lo cual puede ser útil para almacenar y analizar una serie de rasteres.

El conjunto de datos Water Security (WSIM-GLDAS) Monthly Grids, v1 (1948 - 2014) “identifica y caracteriza los excedentes y déficits de agua dulce, y los parámetros que determinan estas anomalías, en intervalos mensuales durante el período de enero de 1948 a diciembre de 2014” (ISciences and Center For International Earth Science Information Network-CIESIN-Columbia University 2022b). El conjunto de datos puede descargarse desde el sitio web de NASA SEDAC. Las descargas de datos de WSIM-GLDAS están organizadas por una combinación de variables temáticas (excedente/déficit compuesto, temperatura, PETmE, escorrentía, humedad del suelo, precipitación) y períodos de integración (una agregación temporal) (1, 3, 6, 12 meses). Cada combinación de variable-integración consta de un archivo raster NetCDF (.nc) (con una dimensión de tiempo que contiene una capa raster para cada uno de los 804 meses entre enero de 1948 y diciembre de 2014. Algunas variables también contienen múltiples atributos, cada uno con su propia serie temporal. Por lo tanto, este es un archivo grande que puede tardar mucho tiempo en descargarse y puede causar problemas de memoria en ciertos sistemas. Esto se considera data GRANDE.

Revisión de Conocimientos
  1. ¿Cómo describirías mejor el ciclo del agua?
    1. Un período prolongado de poca o ninguna lluvia.
    2. Baja precipitación combinada con temperaturas récord.
    3. La circulación del agua en y sobre la superficie de la Tierra.
    4. Un ciclo que ocurre debido a la sequía.
  2. ¿Qué intervenciones humanas afectan el ciclo del agua (selecciona todas las que apliquen)?
    1. Emisiones de gases de efecto invernadero.
    2. Cambios en el uso del suelo.
    3. Desarrollo de presas y embalses.
    4. Sobreexplotación de aguas subterráneas.
  3. ¿Qué es un déficit de precipitación?
    1. Un período de precipitación por debajo del promedio.
    2. Un período prolongado de poca o ninguna lluvia.
    3. Un período de reacciones en cadena.
    4. Un período de precipitación por encima del promedio.

Adquisición de los Datos

Revisión de Ciencia de Datos

El conjunto de datos Water Security (WSIM-GLDAS) Monthly Grids utilizado en esta lección está alojado por el Centro de Datos y Aplicaciones Socioeconómicas de la NASA (SEDAC), uno de los varios Centros de Archivo Activo Distribuidos (DAACs). SEDAC alberga una variedad de productos de datos que incluyen datos de población geoespacial, asentamientos humanos e infraestructura, exposición y vulnerabilidad al cambio climático, y datos de calidad del aire, agua y uso del suelo basados en satélites. Para descargar datos alojados por SEDAC, se requiere tener una cuenta gratuita de NASA EarthData. Puedes crear una cuenta aquí: NASA EarthData.

Para esta lección, trabajaremos con el archivo NetCDF WSIM-GLDAS dataset Composite Anomaly Twelve-Month Return Period. Este contiene variables de agua, excedente y “anomalía compuesta” combinada con un período de integración de 12 meses. El período de integración representa el tiempo durante el cual se promedian los valores de anomalía. La integración de 12 meses promedia déficits de agua (sequías), excedentes (inundaciones) y combinados (presencia de sequías e inundaciones) durante un período de 12 meses que finaliza con el mes especificado. Comenzamos con la agregación de 12 meses porque es una instantánea de anomalías para todo el año, lo que la hace útil para obtener una comprensión de un año completo; una vez que identificamos períodos de tiempo de interés en los datos, podemos examinar más de cerca utilizando los períodos de integración de 3 meses o 1 mes.

Cuando trabajas en tu sistema local puedes descargar directamente WSIM-GLDAS desde el sitio web de SEDAC, sin embargo, estamos trabajando en la nube, por lo que nos conectaremos al cubo de datos TOPS-SCHOOL S3. La documentación del conjunto de datos describe las variables compuestas como características clave de WSIM-GLDAS que combinan “los períodos de retorno de múltiples parámetros de agua en índices compuestos de excedentes y déficits de agua en general (ISciences and Center For International Earth Science Information Network-CIESIN-Columbia University 2022a)”. Los archivos de anomalía compuesta presentan los datos en términos de la frecuencia con la que ocurren; o un “período de retorno”. Por ejemplo, un período de retorno de déficit de 25 significa una sequía tan severa que solo esperaríamos que ocurra una vez cada 25 años. Por favor, procede a descargar el archivo.

Lectura de los Datos

Revisión de Ciencia de Datos

Esta lección utiliza los paquetes stars, terra, sf, lubridate, httr, cubelyr, exactextractr, ggplot2, data.table, y kableExtra. Si deseas aprender más sobre las funciones utilizadas en esta lección, puedes utilizar las guías de ayuda en los sitios web de sus paquetes.

Los paquetes stars y terra están diseñados para trabajar con datos espaciales raster grandes y complejos, mientras que el paquete sf te permite manejar y analizar datos espaciales vectoriales (más sobre esto más adelante). El paquete cubelyr te ayuda a crear y analizar cubos de datos multidimensionales, lo que facilita explorar conjuntos de datos complejos y descubrir patrones. El paquete lubridate facilita trabajar con fechas y horas en R, y httr proporciona métodos para que los usuarios interactúen con bases de datos y sitios web en línea para descargar datos directamente. exactextractr realiza estadísticas zonales rápidas y eficientes que resumen datos raster y vectoriales. data.table facilita el procesamiento estándar de datos, y finalmente, ggplot2 y kableExtra se utilizan para crear gráficos, mapas y tablas.

Asegúrate de que estén instalados antes de comenzar a trabajar con el código en este documento. Si deseas aprender más sobre las funciones utilizadas en esta lección, puedes utilizar las guías de ayuda en los sitios web de sus paquetes.

A continuación, instala y carga los paquetes de R necesarios para este ejercicio. Los paquetes que estamos utilizando hoy se enumeran a continuación, pero estamos utilizando la imagen Docker de TOPS-SCHOOL para nuestro análisis, por lo que el entorno en la nube que estamos utilizando hoy ya está configurado con todo lo que necesitamos.

packages_to_check <- c("stars", "terra", "sf", "cubelyr", "lubridate", "httr", "data.table", "exactextractr", "ggplot2", "kableExtra", "jsonlite", "tmap", "basemaps", "sp")

# Check and install packages
for (package_name in packages_to_check) {
  if (!package_name %in% rownames(installed.packages())) {
    install.packages(package_name)
    cat(paste("Package", package_name, "installed.\n"))
  } else {
    cat(paste("Package", package_name, "is already installed.\n"))
  }
  library(package_name, character.only = TRUE)
}

Cargar Datos y Selección de Atributos

Primero nos conectaremos a nuestro “cubo” de almacenamiento en la nube de TOPS-SCHOOL para descargar el archivo WSIM-GLDAS. Si estuvieras trabajando en un sistema local, utilizarías el sitio web de SEDAC para descargar los archivos, pero eso no es factible en un taller en la nube en vivo. Leer el archivo con el argumento proxy = TRUE te permite inspeccionar los elementos básicos del archivo sin leer el archivo completo en la memoria. Los conjuntos de datos raster multidimensionales pueden ser extremadamente grandes y detener tu entorno de computación si tienes limitaciones de memoria.

# leer el archivo de integración de 12 meses WSIM-GLDAS en nuestro cubo s3
wsim_gldas <- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_12mo.nc", proxy = TRUE)
deficit, deficit_cause, surplus, surplus_cause, both, 
print(wsim_gldas)
stars_proxy object with 5 attributes in 5 file(s):
$deficit
[1] "[...]/composite_12mo.nc:deficit"

$deficit_cause
[1] "[...]/composite_12mo.nc:deficit_cause"

$surplus
[1] "[...]/composite_12mo.nc:surplus"

$surplus_cause
[1] "[...]/composite_12mo.nc:surplus_cause"

$both
[1] "[...]/composite_12mo.nc:both"

dimension(s):
     from   to offset delta  refsys                    values x/y
x       1 1440   -180  0.25  WGS 84                      NULL [x]
y       1  600     90 -0.25  WGS 84                      NULL [y]
time    1  793     NA    NA POSIXct 1948-12-01,...,2014-12-01    

Ahora podemos usar el comando print para ver la información básica. La salida nos dice que tenemos 5 atributos (déficit, causa del déficit, excedente, causa del excedente, ambos) y 3 dimensiones. Las primeras 2 dimensiones son las extensiones espaciales (x/y–longitud/latitud) y el tiempo es la 3ª dimensión.

Esto significa que el número total de capas raster individuales en este NetCDF es 3965 (5 atributos x 793 pasos de tiempo/meses). De nuevo, data grande.

Podemos gestionar este archivo grande seleccionando una sola variable; en este caso “both” (ambos). Leeremos los datos de nuevo; esta vez con proxy = FALSE para leer los datos en memoria real, pero solo seleccionando la capa de déficit con sub = ‘surplus’ para que sea más manejable. Esto toma aproximadamente 45-60 segundos en la plataforma de 2i2c.

# leer solo la capa de déficit usando proxy = false
# wsim_gldas <- stars::read_stars("data/composite_1mo.nc", sub = 'deficit')

wsim_gldas <- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_12mo.nc", proxy = FALSE, sub = 'surplus')
surplus, 
print(wsim_gldas)
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
         Min.   1st Qu.   Median      Mean  3rd Qu. Max.  NA's
surplus   -60 -7.275392 -2.70342 -4.822883 2.431405   60 94088
dimension(s):
     from   to offset delta  refsys                    values x/y
x       1 1440   -180  0.25  WGS 84                      NULL [x]
y       1  600     90 -0.25  WGS 84                      NULL [y]
time    1  793     NA    NA POSIXct 1948-12-01,...,2014-12-01    

Selección de Tiempo

Especificar un rango temporal de interés hará que el tamaño del archivo sea más pequeño y manejable. Seleccionaremos todos los años para el rango 2000-2014. Esto se puede lograr generando una secuencia para cada año entre diciembre de 2000 y diciembre de 2014, y luego pasando esa lista de fechas, en forma de vector, usando la función “filtro”. Recuerde, estamos utilizando la integración de 12 meses de WSIM-GLDAS. Esto significa que cada paso de tiempo enumerado promedia el déficit de los 12 meses anteriores. Por lo tanto, si seleccionamos solo una secuencia de meses de diciembre que abarca el período 2000-2014, cada capa resultante representará el déficit promedio de ese año.

# generar un vector de fechas para subconjuntos
keeps<-seq(lubridate::ymd("2000-12-01"),
           lubridate::ymd("2014-12-01"),
           by = "year")
#cambiar tipo de datos a POSIXct
keeps <- as.POSIXct(keeps)
# filtrar el conjunto de datos usando el vector de fechas
wsim_gldas_2000_2014 <- dplyr::filter(wsim_gldas, time %in% keeps)
# volver a verificar la información básica
print(wsim_gldas_2000_2014)
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
         Min.   1st Qu.   Median      Mean  3rd Qu. Max.  NA's
surplus   -60 -2.767541 2.250129 -2.196594 3.180827   60 94088
dimension(s):
     from   to offset delta  refsys                    values x/y
x       1 1440   -180  0.25  WGS 84                      NULL [x]
y       1  600     90 -0.25  WGS 84                      NULL [y]
time    1   15     NA    NA POSIXct 2000-12-01,...,2014-12-01    

Ahora nos hemos reducido a un único atributo (“déficit”) con 15 pasos de tiempo. Podemos echar un vistazo a los primeros 6 años pasando el objeto por “slice” y luego por “trama”.

wsim_gldas_2000_2014 |>
  # corta los primeros 6 pasos de tiempo
  dplyr::slice(index = 1:6, along = "time") |>
  # trazarlos con algunas quiebras de prueba para la leyenda
  plot(key.pos = 1,  breaks = c(-60, -30, -5,  0, 5, 30, 70), key.lab = "Pariodos de retorno de ambos (excedentes y deficits)") 

# breaks = c(0, -1, -5, -10, -20, -30, -50, -100, -200, -500),
# breaks = c(0, 1, 5, 10, 20, 30, 50),

Aunque ahora hemos reducido los datos a un solo atributo con un tiempo de interés restringido, podemos ir un paso más allá y limitar la extensión espacial a un país o estado de interés.

Verificación de conocimientos
  1. ¿Cuál de estos describe mejor un conjunto de datos ráster?
    1. Un tipo de datos geográficos formados por celdas cuadradas.

    2. Una tabla o lista de números geográficos.

    3. Una región geográfica de interés.

    4. Un tipo de datos geográficos formados por líneas y puntos

  2. ¿Cuál de las siguientes afirmaciones es cierta sobre la información que pueden almacenar los rásteres? (seleccione todas las que correspondan)
    1. Atributos (contenido temático)
    2. Dimensiones (información que expresa información de extensión espacial o temporal)
    3. Coordenadas geográficas
    4. Una lista de números

Selección espacial

Revisión de ciencia de datos

GeoJSON es un formato para codificar, almacenar y compartir datos geográficos como datos vectoriales, es decir, puntos, líneas y polígonos. Se utiliza comúnmente en aplicaciones de mapas web para representar características geográficas y sus atributos asociados. Los archivos GeoJSON son fácilmente legibles tanto por humanos como por computadoras, lo que los convierte en una opción popular para compartir datos geográficos a través de Internet.

Podemos recortar espacialmente la pila de ráster con algunos métodos diferentes. Las opciones incluyen el uso de un cuadro delimitador en el que se especifican las coordenadas geográficas exteriores (xmin, ymin, xmax, ymax), el uso de otro objeto ráster o el uso de un límite vectorial como un archivo de forma o GeoJSON para recortar la extensión de los datos ráster originales.

En este ejemplo utilizamos un límite vectorial para realizar la tarea de geoprocesamiento de recortar los datos a una unidad administrativa o política. Primero, adquirimos los datos en formato GeoJSON para Ecuador desde la API geoBoundaries. (Tenga en cuenta que también es posible descargar los límites vectorizados directamente desde https://www.geoboundaries.org/ en lugar de utilizar la API).

Para utilizar la API de geoBoundaries, la URL raíz a continuación se modifica para incluir un código de 3 letras de la Organización Internacional de Estándares utilizado para identificar países (ISO3) y un nivel administrativo para la solicitud de datos. Los niveles administrativos corresponden a unidades geográficas como el País (nivel administrativo 0), el Estado/Provincia (nivel administrativo 1), el Condado/Distrito (nivel administrativo 2), etc.:

https://www.geoboundaries.org/api/current/gbOpen/ISO3/NIVEL/

Revisión de ciencia de datos

Creada por la comunidad y William & Mary geoLab, la base de datos global de límites administrativos políticos geoBoundaries es un recurso en línea con licencia abierta (CC BY 4.0 / ODbL) para límites administrativos (es decir, , estado, condado, provincia) para todos los países del mundo. Desde 2016, este proyecto ha rastreado aproximadamente 1 millón de unidades espaciales dentro de más de 200 entidades; incluidos todos los estados miembros de la ONU.

Para este ejemplo, ajustamos los componentes en negrita de la dirección URL de muestra a continuación para especificar el país que queremos usando el código de país de caracteres ISO3 para Ecuador (ECU) y el nivel administrativo del estado deseado (ADM2 ).

# realizar la solicitud al sitio web de geoboundarie para los límites de Ecuador al nivel de administracion 2.
ecuador_ad2 <- httr::GET("https://www.geoboundaries.org/api/current/gbOpen/ECU/ADM2/")

En la línea de código anterior, usamos httr:GET para obtener metadatos de la URL. Asignamos el resultado a una nueva variable llamada “ecuador_ad2”. A continuación examinaremos el “contenido”.

# analizar el contenido en un formato legible
ecuador_ad2<- httr::content(ecuador_ad2)
# ecuador_ad3 <- fromJSON(ecuador_ad3)
# mire las etiquetas para obtener información disponible
names(ecuador_ad2)
 [1] "boundaryID"                "boundaryName"             
 [3] "boundaryISO"               "boundaryYearRepresented"  
 [5] "boundaryType"              "boundaryCanonical"        
 [7] "boundarySource"            "boundaryLicense"          
 [9] "licenseDetail"             "licenseSource"            
[11] "boundarySourceURL"         "sourceDataUpdateDate"     
[13] "buildDate"                 "Continent"                
[15] "UNSDG-region"              "UNSDG-subregion"          
[17] "worldBankIncomeGroup"      "admUnitCount"             
[19] "meanVertices"              "minVertices"              
[21] "maxVertices"               "meanPerimeterLengthKM"    
[23] "minPerimeterLengthKM"      "maxPerimeterLengthKM"     
[25] "meanAreaSqKM"              "minAreaSqKM"              
[27] "maxAreaSqKM"               "staticDownloadLink"       
[29] "gjDownloadURL"             "tjDownloadURL"            
[31] "imagePreview"              "simplifiedGeometryGeoJSON"

El contenido analizado contiene 32 componentes. El elemento 29 es un enlace directo al archivo GeoJSON (gjDownloadURL) donde se encuentran los datos del límite del vector. A continuación obtendremos el GeoJSon y comprobaremos los resultados.

# leer directamente en el geojson con sf desde el servidor de geoserver

ecuador_ad2 <- sf::st_read(ecuador_ad2$gjDownloadURL)
Reading layer `geoBoundaries-ECU-ADM2' from data source 
  `https://github.com/wmgeolab/geoBoundaries/raw/9469f09/releaseData/gbOpen/ECU/ADM2/geoBoundaries-ECU-ADM2.geojson' 
  using driver `GeoJSON'
Simple feature collection with 223 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -92.00897 ymin: -5.016157 xmax: -75.18715 ymax: 1.681835
Geodetic CRS:  WGS 84
# revisa las imágenes
plot(sf::st_geometry(ecuador_ad2))

# ver el nombre de las areas administrativas 2
ecuador_ad2$shapeName
  [1] "24 de Mayo"                  "Aguarico"                   
  [3] "Alausi"                      "Alfredo Baquerizo Moreno"   
  [5] "Ambato"                      "Antonio Ante"               
  [7] "Arajuno"                     "Archidona"                  
  [9] "Arenillas"                   "Atacames"                   
 [11] "Atahualpa"                   "Azogues"                    
 [13] "Baba"                        "Babahoyo"                   
 [15] "Balao"                       "Balsas"                     
 [17] "Balzar"                      "Baños de Agua Santa"        
 [19] "Biblian"                     "Bolivar"                    
 [21] "Bolivar"                     "Buena Fe"                   
 [23] "Caluma"                      "Calvas"                     
 [25] "Camilo Ponce Enriquez"       "Cañar"                      
 [27] "Carlos Julio Arosemena Tola" "Cascales"                   
 [29] "Catamayo"                    "Cayambe"                    
 [31] "Celica"                      "Centinela del Condor"       
 [33] "Cevallos"                    "Chaguarpamba"               
 [35] "Chambo"                      "Chilla"                     
 [37] "Chillanes"                   "Chimbo"                     
 [39] "Chinchipe"                   "Chone"                      
 [41] "Chordeleg"                   "Chunchi"                    
 [43] "Colimes"                     "Colta"                      
 [45] "Cotacachi"                   "Crnel. Marcelino Maridueña" 
 [47] "Cuenca"                      "Cumanda"                    
 [49] "Cuyabeno"                    "Daule"                      
 [51] "Deleg"                       "Duran"                      
 [53] "Echeandia"                   "El Carmen"                  
 [55] "El Chaco"                    "El Guabo"                   
 [57] "El Pan"                      "El Pangui"                  
 [59] "El Piedrero"                 "El Tambo"                   
 [61] "El Triunfo"                  "Eloy Alfaro"                
 [63] "Empalme"                     "Esmeraldas"                 
 [65] "Espejo"                      "Espindola"                  
 [67] "Flavio Alfaro"               "Giron"                      
 [69] "Gnral. Antonio Elizalde"     "Gonzalo Pizarro"            
 [71] "Gonzanama"                   "Guachapala"                 
 [73] "Gualaceo"                    "Gualaquiza"                 
 [75] "Guamote"                     "Guano"                      
 [77] "Guaranda"                    "Guayaquil"                  
 [79] "Huamboya"                    "Huaquillas"                 
 [81] "Ibarra"                      "Isabela"                    
 [83] "Isidro Ayora"                "Jama"                       
 [85] "Jaramijo"                    "Jipijapa"                   
 [87] "Junin"                       "La Concordia"               
 [89] "La Joya de los Sachas"       "La Libertad"                
 [91] "La Mana"                     "La Troncal"                 
 [93] "Lago Agrio"                  "Las Golondrinas"            
 [95] "Las Lajas"                   "Las Naves"                  
 [97] "Latacunga"                   "Limon Indanza"              
 [99] "Logroño"                     "Loja"                       
[101] "Lomas de Sargentillo"        "Loreto"                     
[103] "Macara"                      "Machala"                    
[105] "Manta"                       "Marcabeli"                  
[107] "Mejia"                       "Mera"                       
[109] "Milagro"                     "Mira"                       
[111] "Mocache"                     "Mocha"                      
[113] "Montalvo"                    "Montecristi"                
[115] "Montufar"                    "Morona"                     
[117] "Muisne"                      "Nabon"                      
[119] "Nangaritza"                  "Naranjal"                   
[121] "Naranjito"                   "Nobol"                      
[123] "Olmedo"                      "Olmedo"                     
[125] "Oña"                         "Orellana"                   
[127] "Otavalo"                     "Pablo Sexto"                
[129] "Pajan"                       "Palanda"                    
[131] "Palenque"                    "Palestina"                  
[133] "Pallatanga"                  "Palora"                     
[135] "Paltas"                      "Pangua"                     
[137] "Paquisha"                    "Pasaje"                     
[139] "Pastaza"                     "Patate"                     
[141] "Paute"                       "Pedernales"                 
[143] "Pedro Carbo"                 "Pedro Moncayo"              
[145] "Pedro Vicente Maldonado"     "Penipe"                     
[147] "Pichincha"                   "Pimampiro"                  
[149] "Piñas"                       "Pindal"                     
[151] "Playas"                      "Portovelo"                  
[153] "Portoviejo"                  "Pucara"                     
[155] "Puebloviejo"                 "Puerto Lopez"               
[157] "Puerto Quito"                "Pujili"                     
[159] "Putumayo"                    "Puyango"                    
[161] "Quero"                       "Quevedo"                    
[163] "Quijos"                      "Quilanga"                   
[165] "Quininde"                    "Quinsaloma"                 
[167] "Quito"                       "Riobamba"                   
[169] "Rioverde"                    "Rocafuerte"                 
[171] "Rumiñahui"                   "Salcedo"                    
[173] "Salinas"                     "Salitre"                    
[175] "Samborondon"                 "San Cristobal"              
[177] "San Fernando"                "San Jacinto de Yaguachi"    
[179] "San Juan Bosco"              "San Lorenzo"                
[181] "San Miguel"                  "San Miguel de Los Bancos"   
[183] "San Miguel de Urcuqui"       "San Pedro de Huaca"         
[185] "San Pedro de Pelileo"        "San Vicente"                
[187] "Santa Ana"                   "Santa Clara"                
[189] "Santa Cruz"                  "Santa Elena"                
[191] "Santa Isabel"                "Santa Lucia"                
[193] "Santa Rosa"                  "Santiago"                   
[195] "Santiago de Pillaro"         "Santo Domingo"              
[197] "Saquisili"                   "Saraguro"                   
[199] "Sevilla de Oro"              "Shushufindi"                
[201] "Sigchos"                     "Sigsig"                     
[203] "Simon Bolivar"               "Sozoranga"                  
[205] "Sucre"                       "Sucua"                      
[207] "Sucumbios"                   "Suscal"                     
[209] "Taisha"                      "Tena"                       
[211] "Tisaleo"                     "Tiwintza"                   
[213] "Tosagua"                     "Tulcan"                     
[215] "Urdaneta"                    "Valencia"                   
[217] "Ventanas"                    "Vinces"                     
[219] "Yacuambi"                    "Yantzaza"                   
[221] "Zamora"                      "Zapotillo"                  
[223] "Zaruma"                     

Al examinarlo, vemos que este GeoJSon incluye las areas administrativas y territorios. (Por supuesto, también podría simplificarse a otras áreas de interés simplemente adaptando el código siguiente).

Primero creamos una lista de las geografías que deseamos seleccionar y las asignamos a una variable llamada “estados”. A continuación, reasignamos nuestra variable “ecuador_ad2_sub” para incluir solo las areas seleccionadas, y finalmente trazamos los resultados.

# crear un vector de territorios que queremos en nuestro límite
estados<-
  c("Tiwintza", "Tosagua", "Tulcan" ,"Urdaneta" , "Valencia", "Riobamba" ,"Ventanas","Vinces","Yacuambi","Yantzaza","Zamora","Zapotillo")

# seleccione todos los estados y territorios que no están en la lista anterior
# en R el "!" significa NO o lo OPUESTO
ecuador_ad2_sub<-ecuador_ad2[(ecuador_ad2$shapeName %in% estados),]
# en R el "!" significa NO o lo OPUESTO
# ecuador_ad_sub<-ecuador_ad2[!(ecuador_ad2$shapeName %in% estados),]
# revisa las imágenes
plot(sf::st_geometry(ecuador_ad2_sub))

Podemos llevar esto un paso más allá y seleccionar un solo estado para el análisis. Aquí utilizamos un método ligeramente diferente creando un nuevo objeto llamado “ecu_adm1” subdividiendo por nombre.

# extraer solo un  límite
ecu_adm1 <- ecuador_ad2[ecuador_ad2$shapeName == "Zapotillo",]
# revisa las imágenes
plot(sf::st_geometry(ecu_adm1))

Desde aquí podemos recortar la dimension escapical del ráster WSIM-GLDAS utilizando el límite almacenado. Puedes llamar a la función sf::st_crop() para recortar la capa WSIM-GLDAS, pero como ves a continuación, de manera más simple, puedes usar la indexación entre corchetes [ ] para recortar un objeto stars con un objeto de sf.

# recortar el archivo wsim-gldas hasta la extensión del límite de Ecuador
wsim_gldas <- wsim_gldas[ecuador_ad2]
print(wsim_gldas)
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
              Min.   1st Qu.   Median     Mean  3rd Qu. Max.  NA's
surplus  -57.31307 -3.759839 2.188979 0.238033 3.910136   60 83028
dimension(s):
     from  to offset delta  refsys                    values x/y
x     352 420   -180  0.25  WGS 84                      NULL [x]
y     354 381     90 -0.25  WGS 84                      NULL [y]
time    1 793     NA    NA POSIXct 1948-12-01,...,2014-12-01    

Finalmente, visualizamos el último paso de tiempo en el conjunto de datos WSIM-GLDAS (1 de diciembre de 2014) y lo representamos con una superposición del límite de Ecuador para realizar una verificación visual de nuestro procesamiento.

# elimina los primeros 15 pasos de tiempo de wsim-gldas y trazalos
wsim_gldas |>
  dplyr::slice(index = 15, along = "time") |>
  # visualizar
  plot(reset = FALSE) # , breaks = c(0, 5, 10, 20, 30, 50, 60)
# agrega el límite de Ecuador en la parte superior
plot(sf::st_geometry(ecuador_ad2),
     add = TRUE,
     lwd = 1,
     fill = NA,
     border = 'purple')

¡Voilá! Recortamos con éxito el WSIM-GLDAS hasta la extensión de Ecuador y creamos un mapa superpuesto con ambos conjuntos de datos para verificar los resultados. Si estaba realizando más análisis o deseaba compartir su trabajo con colegas, es posible que desee guardar el WSIM-GLDAS procesado en el disco.

Los archivos ráster multidimensionales (déficit, tiempo, latitud, longitud) se pueden guardar con la función stars::write_mdim y los datos vectoriales se pueden guardar usando sf::st_write.

# escribir el archivo wsim-gldas procesado en el disco como netcdf
# stars::write_mdim(wsim_gldas_tex, "output/wsim_gldas_ecu.nc")
# escribir el límite de Ecuador en el disco
# sf::st_write(ecuador_ad2, "output/ecuador.geojson")

El tamaño del conjunto de datos preprocesado es de 1,6 MB en comparación con el conjunto de datos original de 1,7 GB. Esto es mucho más manejable en entornos de nube, talleres y plataformas Git. t

Visualizaciones avanzadas e integraciones de datos

Ahora que hemos introducido los conceptos básicos de manipulación y visualización de WSIM-GLDAS, podemos explorar visualizaciones e integraciones de datos más avanzadas. Limpiemos el espacio de trabajo y comencemos de nuevo con el mismo Período de retorno de doce meses de anomalía compuesta WSIM-GLDAS que usamos anteriormente. Subdividiremos espacialmente los datos para cubrir a Ecuador, lo que ayudará a minimizar nuestra huella de memoria. Podemos reducir aún más nuestra sobrecarga de memoria leyendo solo la variable que queremos analizar. En este caso, podemos leer solo el atributo déficit del archivo del período de retorno de doce meses de anomalía compuesta de WSIM-GLDAS, en lugar de leer el NetCDF completo con todos sus atributos.

Revisión de codificación

La memoria de acceso aleatorio (RAM) es donde se almacenan temporalmente los datos y programas que están actualmente en uso. Permite que la computadora acceda rápidamente a los datos, haciendo que todo lo que haga en la computadora sea más rápido y eficiente. A diferencia del disco duro, que almacena datos de forma permanente, la RAM pierde todos sus datos cuando se apaga la computadora. La RAM es como la memoria a corto plazo de una computadora y le ayuda a realizar múltiples tareas a la vez.

Para este ejercicio, podemos recorrer rápidamente pasos de preprocesamiento similares que realizamos anteriormente en esta lección y luego pasar a visualizaciones e integraciones más avanzadas. Vuelva a leer los datos de integración originales de 12 meses, filtre con una lista de fechas para cada diciembre que abarca del 1994 al 2001 y luego recorte los datos ráster con el límite de Ecuador utilizando nuestro objeto geoBoundaries.

# leer en la capa wsim-gldas del depósito SCHOOL s3
# wsim_gldas <- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_12mo.nc", proxy = FALSE, sub = 'surplus')

# generar un vector de fechas para subconjuntos
keeps<-seq(lubridate::ymd("1994-12-01"),
           lubridate::ymd("2001-12-01"), 
           by = "year")

#cambiar tipo de datos a POSIXct
keeps <- as.POSIXct(keeps)
# filtrar usando ese vector
wsim_gldas <- dplyr::filter(wsim_gldas, time %in% keeps)

Verifique nuevamente la información del objeto.

# check the basic information again
print(wsim_gldas)
stars object with 3 dimensions and 1 attribute
attribute(s):
              Min.  1st Qu.   Median     Mean  3rd Qu. Max.  NA's
surplus  -28.33879 2.153126 3.824581 7.591903 7.982027   60 12840
dimension(s):
     from to offset delta  refsys                    values x/y
x       1 69 -92.25  0.25  WGS 84                      NULL [x]
y       1 28   1.75 -0.25  WGS 84                      NULL [y]
time    1  8     NA    NA POSIXct 1994-12-01,...,2001-12-01    

Querrá revisar la copia impresa para asegurarse de que se vea bien.

  • ¿Contiene las variables que esperabas?

  • ¿Los valores de las variables parecen plausibles?

Otros análisis descriptivos básicos son útiles para verificar y comprender sus datos. Uno de ellos es producir una distribución de frecuencia (también conocida como histograma), que se analiza a continuación.

Verificación de conocimientos
  1. ¿Cuáles son las dos formas en que redujimos nuestra huella de memoria cuando cargamos el conjunto de datos WSIM-GLDAS para esta lección? (seleccione todas las que correspondan)
    1. Subconjunto espacial de los datos al Ecuador
    2. Leyendo en sólo un año de datos
    3. Leer sólo el atributo de interés (es decir, ‘surplus’)
    4. Todo lo anterior.

Serie temporal anual de Ecuador

Las propiedades básicas de los datos revisadas en el paso anterior son útiles para el análisis exploratorio de datos, pero debemos realizar una inspección más profunda. Podemos comenzar nuestra exploración visual de los excedentes anual en Ecuador creando un mapa que ilustre el período de retorno del excedente para cada uno de los años en el objeto WSIM-GLDAS.

# cargar los datos base del límite de Ecuador
ggplot2::ggplot()+
  # visualizar el objeto
  stars::geom_stars(data = wsim_gldas)+
  # establecer coordenadas iguales para los ejes
  ggplot2::coord_equal()+
  # crear múltiples paneles para cada paso de tiempo
  ggplot2::facet_wrap(~time)+
  # trazar el límite de Ecuador solo con el contorno
  # ggplot2::geom_sf(fill = NA, lwd=0)+
  # agregar la paleta para wsim-gldas
  ggplot2::scale_fill_stepsn(
    colors = c(
    # -29 a 0
    '#FFC300', 
    # 0 a 5
    '#DDFFFF',
    # 5 a 10
    '#00BFFF',
    # 10 a 30
    '#0000FF',
    # 30 a 60
    '#191970'), 
    # establecer las diviciones en la paleta
    breaks = c(-29, 0, 5, 10, 30, 60))+
  # agregar etiquetas de trama
  ggplot2::labs(
    title="Anomalías del Excedente Medio Anual del Ecuador",
    subtitle = "Utilizando datos observados de 12 meses integrados WSIM-GLDAS para 1994-2001",
    fill = "Período de retorno del excedente"
  )+
  # establecer el tema mínimo
  ggplot2::theme_minimal()+
  # desactivar algunos elementos gráficos adicionales
  ggplot2::theme(
    axis.title.x=ggplot2::element_blank(),
    axis.text.x=ggplot2::element_blank(),
    axis.ticks.x=ggplot2::element_blank(),
    axis.title.y=ggplot2::element_blank(),
    axis.text.y=ggplot2::element_blank(),
    axis.ticks.y=ggplot2::element_blank())

Esta visualización muestra que hubo varios eventos de inundacion importantes (como lo indican los valores del período de retorno del exceso en azul oscuro) durante el período 1994-2001. Los eventos de inundacion importantes incluyeron el noroeste en 1998, el sur en 2001, Las inundaciones de los años 1998 y 2001 son particularmente graves y generalizados, con períodos de retorno superiores a 30 años que abarcan varios estados. Según las normas históricas, ¡solo deberíamos esperar sequías tan fuertes cada 30 o 60 años!

Verificación de conocimientos
  1. Los mapas de salida de Anomalías del Execdente Medio Anual para Ecuador indican que…
    1. El exceso más significativo en 1998 se localiza en el oeste de Ecuador.
    2. El exceso menos significativo de 2000 se sitúa en el Sur.
    3. El exceso más significativo en 1996 se localiza en los estados vecinos.
    4. El déficit menos significativo en 2000 se sitúa en el sudeste.

Serie temporal mensual

Podemos obtener una visión más detallada de estos eventos de sequía utilizando el conjunto de datos compuesto WSIM-GLDAS de 1 mes y recortando los datos en una extensión espacial más pequeña que coincida con uno de los eventos que hemos observado en el gráfico anterior. Echemos un vistazo más de cerca a la sequía de California de 2014.

La sequía en las noticias

La sequía de California de 2012-2014 fue la peor en 1200 años (“Evidence Suggests California’s Drought Is the Worst in 1,200 Years” 2014). Esta sequía causó problemas a los propietarios e incluso conflictos entre los agricultores y el salmón salvaje. El gobernador Jerry Brown declaró una emergencia por sequía y pidió a los residentes que redujeran el consumo de agua en un 20%. El uso de agua aumentó un 8% en mayo de 2014 en comparación con 2013, en lugares como la costa de California y Los Ángeles. Debido a la escasez de agua, el estado votó a favor de multar a quienes desperdician agua hasta $500 dólares. La sequía también afectó a los residentes de manera diferente según su situación económica. Por ejemplo, en el condado de El Dorado, ubicado en una zona rural al este de Sacramento, los residentes se duchaban con baldes y los residentes rurales informaron que los pozos, de los que dependen para obtener agua dulce, se estaban secando. El gobierno federal finalmente anunció una ayuda de emergencia por sequía de $9,7 millones para esas áreas (Sanders 2014).

Además, había miles de salmones adultos que luchaban por sobrevivir en el río Klamath en el norte de California, donde el agua era escasa y cálida debido al desvío del flujo del río hacia el Valle Central, una zona agrícola donde se cultivan almendros. Las almendras son uno de los cultivos más importantes de California; el estado produce el 80% de las almendras del mundo. Sin embargo, el salmón, que migra río arriba, podría contraer una enfermedad llamada pudrición branquial, que prospera en aguas cálidas y que ya mató a decenas de miles de Chinook en 2002. Esta enfermedad se estaba propagando nuevamente entre la población de salmón debido a esta asignación de agua, afectando a los nativos locales. Tribus americanas que dependen del pescado (Bland 2014).

Para limitar la cantidad de memoria informática necesaria para la operación, primero borraremos elementos del espacio de trabajo en memoria y luego recargaremos un archivo compuesto más pequeño; comenzaremos eliminando el objeto compuesto de 12 meses.

# eliminar el objeto wsim grande del entorno
rm(wsim_gldas_1mo)
rm(wsim_gldas_ecu)
rm(wsim_gldas_2000_2014)
rm(ecuador_ad2_sub)
rm(ecu_adm1)

# borrar la memoria; R tiene un recolector de basura automatizado para borrar la memoria no utilizada
# pero puedes invocarlo usando el comando gc()
gc()
          used (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells 1414796 75.6    2697790  144.1   2697790  144.1
Vcells 4827865 36.9  634567277 4841.4 750705718 5727.5

Nuevamente nos conectaremos al depósito de nube SCHOOL y leeremos el archivo de integración de 1 mes. La estructura es la misma que la integración de 12 meses, por lo que continuaremos y cargaremos solo la capa deficitaria y verificaremos el objeto con print.

# conéctese al depósito y descargue la integración de 1 mes
# leer en la capa wsim-gldas del depósito SCHOOL s3
wsim_gldas_1mo_98 <- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_1mo.nc", proxy = FALSE, sub = 'surplus')
surplus, 
print(wsim_gldas_1mo_98)
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
         Min.   1st Qu.    Median      Mean   3rd Qu. Max.  NA's
surplus   -60 -5.997547 -4.314255 -5.601294 -3.739346   60 98724
dimension(s):
     from   to offset delta  refsys                    values x/y
x       1 1440   -180  0.25  WGS 84                      NULL [x]
y       1  600     90 -0.25  WGS 84                      NULL [y]
time    1  804     NA    NA POSIXct 1948-01-01,...,2014-12-01    

Una vez más, subdividiremos la dimensión temporal para nuestro período de interés. Sin embargo, esta vez queremos todos los meses de 1998 para poder observar más de cerca la inundacion en Ecuador.

# generar un vector de fechas para subconjuntos
keeps<-seq(lubridate::ymd("1998-01-01"),
           lubridate::ymd("1998-12-01"), 
           by = "month")
# cambiar tipo de datos a POSIXct
keeps <- as.POSIXct(keeps)
# filtrar usando ese vector
wsim_gldas_1mo_98 <- dplyr::filter(wsim_gldas_1mo_98, time %in% keeps)
# verifica la información
print(wsim_gldas_1mo_98)
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
         Min.   1st Qu.    Median      Mean  3rd Qu. Max.  NA's
surplus   -60 -4.030512 -3.502258 0.2928404 3.956679   60 98724
dimension(s):
     from   to offset delta  refsys                    values x/y
x       1 1440   -180  0.25  WGS 84                      NULL [x]
y       1  600     90 -0.25  WGS 84                      NULL [y]
time    1   12     NA    NA POSIXct 1998-01-01,...,1998-12-01    

Ahora tenemos 12 rásteres con datos mensuales para 1998. Acerquémonos a Ecuador y veamos cómo progresó esta inundacion a lo largo del año.

# subconjunto wsim-gldas hasta la extensión del límite de california
wsim_gldas_1mo_98 <- wsim_gldas_1mo_98[ecuador_ad2]
# dar etiquetas bonitas a la dimensión del tiempo
wsim_gldas_1mo_98 <-
  wsim_gldas_1mo_98 |>
    stars::st_set_dimensions("time", values = month.name)

# parcelas mensuales de Ecuador
# cargar los datos base de Ecuador
ggplot2::ggplot()+
  # cargar el objeto wsim stars
  stars::geom_stars(data = wsim_gldas_1mo_98)+
  # relación de aspecto igual para el hacha lat/long
  ggplot2::coord_equal()+
  # crear una trama secundaria para cada paso de tiempo
  ggplot2::facet_wrap(~time)+
  # configurar la paleta wsim
  ggplot2::scale_fill_stepsn(
    colors = c(
    # -60 a -30
    '#9B0039',
    #-30 a -10
    '#FF8D43',
    # -10 a -5
    '#FFC754',
    # -5 a 0
    '#FFF4C7',
    # 0 a 5
    '#DDFFFF',
    # 5 a 10
    '#00BFFF',
    # 10 a 30
    '#0000FF',
    # 30 a 60
    '#191970'), 
    # establecer las diviciones de la paleta
    breaks = c(-60, -30, -10, -5, 0, 5, 10, 30, 60))+
    
  # agregar etiquetas
  ggplot2::labs(
    title="Anomalías del Excedente Medio Mensual del Ecuador",
    subtitle = "Utilizando datos observados de 12 meses WSIM-GLDAS para 1998",
    fill = "Período de retorno del excedente"
  )+
  # comenzar con un tema mínimo base
  ggplot2::theme_minimal()+
  # eliminar elementos adicionales para hacer un mapa más limpio
  ggplot2::theme(
    axis.title.x=ggplot2::element_blank(),
    axis.text.x=ggplot2::element_blank(),
    axis.ticks.x=ggplot2::element_blank(),
    axis.title.y=ggplot2::element_blank(),
    axis.text.y=ggplot2::element_blank(),
    axis.ticks.y=ggplot2::element_blank())

Esta serie de mapas muestra una imagen sorprendente. Ecuado enfrentó enormes excedentes de agua en el Oeste el estado en Marzo, Mayo, Junio. A esto le siguieron déficits de agua en este del estado en Diciembre. Aunque el Centro y el Norte de Ecuado experimentaron cierto alivio en Septiembre, el Noroeste de Ecuador continuó experimentando excedentes hasta Noviembre.

Histogramas mensuales

Revisión de ciencia de datos

Un marco de datos (dataFrame) es una estructura de datos utilizada para almacenar datos tabulares. Organiza datos en filas y columnas. Cada columna puede tener un tipo diferente de datos (numéricos, de caracteres, de factores, etc.) y las filas representan observaciones o casos individuales. Los marcos de datos proporcionan una manera conveniente de trabajar con datos estructurados, lo que los hace esenciales para proyectos de análisis de datos y estadísticas.

Podemos explorar más los datos creando una distribución de frecuencia (también llamada histograma) de las anomalías del déficit para cualquier extensión espacial determinada; Aquí todavía estamos observando las anomalías del excedente en Ecuador. Comenzamos extrayendo los datos de la serie temporal ráster y luego creamos un marco de datos de valores que son más fáciles de manipular en un histograma. Los dataFrame R son datos que se muestran en formato de tabla, que se pueden trazar en gráficos o tablas.

# extraer los valores ráster en un marco de datos
excedente_hist <-  
  wsim_gldas_1mo_98 |>
  as.data.frame(wsim_gldas_1mo_98$surplus)

# eliminar los valores NA
excedente_hist<-stats::na.omit(excedente_hist)

# crear el histograma
ggplot2::ggplot(excedente_hist, ggplot2::aes(surplus))+
  # darle estilo a las barras
  ggplot2::geom_histogram(binwidth = 6, fill = "#325d88")+
  # subtrama para cada paso de tiempo
  ggplot2::facet_wrap(~time)+
  # usa el tema mínimo
  ggplot2::theme_minimal()

Los histogramas empiezan a cuantificar lo que vimos en los mapas de series temporales. Mientras que el mapa muestra dónde ocurren los excedentes, la distribución de frecuencia indica el número de celdas ráster para cada período de retorno del rango de excedente. El número de celdas ráster bajo un excedente de 60 años (período de retorno) es muy alto en la mayoría de los meses, superando con creces cualquier otro valor en el rango.

Verificación de conocimientos
  1. ¿Cuál es otro término para “histograma”?
    1. coropleta
    2. Distribución de frecuencias
    3. Matriz de datos
    4. Diagrama de caja

Resúmenes Zonales

Hasta este punto, hemos descrito la excedented de Ecuador de 1998 examinando el estado en su conjunto. Aunque al mirar los mapas tenemos una idea de lo que está sucediendo en diferentes ciudades o condados, estos no proporcionan resúmenes cuantitativos de las áreas locales.

Las estadísticas zonales son una forma de resumir las celdas de una capa ráster que se encuentran dentro de los límites de otra capa de datos (que puede estar en formato ráster o vectorial). Por ejemplo, agregar períodos de retorno del déficit con otro ráster que represente el tipo de cobertura terrestre o un límite vectorial (shapefile) de países, estados o condados producirá estadísticas descriptivas para esa nueva capa. Estas estadísticas podrían incluir la suma, la media, la mediana, la desviación estándar y el rango.

Revisión de codificación

El paquete R exactextractr (Daniel Baston 2023) resume los valores ráster en agrupaciones o zonas, también conocidas como estadísticas zonales. Las estadísticas zonales ayudan a evaluar las características estadísticas de una determinada región.

El paquete R terra procesa datos geoespaciales rasterizados y ofrece funcionalidades como manipulación de datos, análisis espacial, modelado y visualización, con un enfoque en la eficiencia y la escalabilidad.

Realizaremos nuestras estadísticas zonales utilizando el paquete exactextractr (Daniel Baston 2023). Es la herramienta de estadísticas zonales más rápida, precisa y flexible para el lenguaje de programación R, pero actualmente no tiene métodos predeterminados para el paquete stars y los objetos stars, por lo que cambiaremos a terra para esto. parte de la lección. Notarás que existen ligeras diferencias en la sintaxis para realizar las mismas operaciones en terra que en stars.

Vuelva a leer los datos con terra.

# eliminar el objeto del mes anterior
rm(wsim_gldas_1mo_98)

# crear una colección terra con el archivo wsim-gldas de 1 mes de SEDAC
wsim_gldas_1mo<-terra::sds("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_1mo.nc")

Verificaremos la información básica con print.

print(wsim_gldas_1mo)
class       : SpatRasterDataset 
subdatasets : 5 
dimensions  : 600, 1440 (nrow, ncol)
nlyr        : 804, 804, 804, 804, 804 
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : composite_1mo.nc 
names       : deficit, deficit_cause, surplus, surplus_cause, both 

Puedes ver que la información de terra está organizada de manera diferente, pero toda sigue ahí. Los atributos están numerados en los “subconjuntos de datos” (5) y se enumeran por nombre en la columna “nombres”. La dimensión de tiempo no aparece por nombre, pero “nlyr” nos dice que hay 804 capas para cada atributo/subconjunto de datos. Puede acceder a información más detallada para las dimensiones de tiempo usando terra::time.

Subconjunto simplemente al excedente de agua de WSIM-GLDAS y luego nuevamente a lo largo de la dimensión temporal.

# sacar solo la capa deficitaria
 wsim_gldas_1mo<-wsim_gldas_1mo["surplus"]
# crear una secuencia de fechas para mantener
keeps<-seq(lubridate::ymd("2014-01-01"), lubridate::ymd("2014-12-01"), by = "month")
# subconjunto del objeto terra solo para esos pasos de tiempo del 2014
wsim_gldas_1mo_14<-wsim_gldas_1mo[[terra::time(wsim_gldas_1mo) %in% keeps]]
# etiquetar los pasos de tiempo
names(wsim_gldas_1mo_14) <- keeps
#verifica la información
print(wsim_gldas_1mo_14)
class       : SpatRaster 
dimensions  : 600, 1440, 12  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : composite_1mo.nc:surplus 
varname     : surplus (Composite Surplus Index) 
names       : 2014-01-01, 2014-02-01, 2014-03-01, 2014-04-01, 2014-05-01, 2014-06-01, ... 
time (days) : 2014-01-01 to 2014-12-01 
# crear una secuencia de fechas para mantener
keeps<-seq(lubridate::ymd("1998-01-01"), lubridate::ymd("1998-12-01"), by = "month")
# subconjunto del objeto terra solo para esos pasos de tiempo del 1998
wsim_gldas_1mo_98<-wsim_gldas_1mo[[terra::time(wsim_gldas_1mo) %in% keeps]]
# etiquetar los pasos de tiempo
names(wsim_gldas_1mo_98) <- keeps
#verifica la información
print(wsim_gldas_1mo_98)
class       : SpatRaster 
dimensions  : 600, 1440, 12  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : composite_1mo.nc:surplus 
varname     : surplus (Composite Surplus Index) 
names       : 1998-01-01, 1998-02-01, 1998-03-01, 1998-04-01, 1998-05-01, 1998-06-01, ... 
time (days) : 1998-01-01 to 1998-12-01 
#remover el variable 
rm(wsim_gldas_1mo)

Ahora que hemos seleccionado solo la variable “surplus”, terra enumera el “varname” como “surplus” y muestra los pasos de tiempo en la columna “nombres”. Un poco confuso, pero terra cambia la visualización una vez que el objeto pasa de múltiples variables a solo una.

Hora de realizar el resumen zonal.

# ejecutar la extracción
ecu_resumen_14<-
  exactextractr::exact_extract(
    # los datos ráster
    wsim_gldas_1mo_14, 
    # los datos del límite
    ecuador_ad2, 
    # el cálculo
    'mean', 
    # no mostrar progreso
    progress = FALSE)
# darle etiquetas más bonitas al paso del tiempo
names(ecu_resumen_14)<-lubridate::month(keeps, label = TRUE, abbr = FALSE)
# Incluir 'E14_' a los nombres de columna
names(ecu_resumen_14) <- paste0('E14_', names(ecu_resumen_14))

ecu_resumen_98<-
  exactextractr::exact_extract(
    # los datos ráster
    wsim_gldas_1mo_98, 
    # los datos del límite
    ecuador_ad2, 
    # el cálculo
    'mean', 
    # no mostrar progreso
    progress = FALSE)
# darle etiquetas más bonitas al paso del tiempo
names(ecu_resumen_98)<-lubridate::month(keeps, label = TRUE, abbr = FALSE)
# Incluir 'E14_' a los nombres de columna
names(ecu_resumen_98) <- paste0('E98_', names(ecu_resumen_98))

exactextractr devolverá estadísticas resumidas en el mismo orden que el archivo de límites de entrada, por lo tanto, podemos unir los nombres de las areas admnistrativas de Ecuador a la salida de estadísticas resumidas exactextract para su visualización. También haremos una versión para ver como una tabla para inspeccionar los datos sin procesar. Podemos echar un vistazo rápido a los primeros 10 areas para ver su período medio de retorno del déficit de enero a junio.

# vincular los medios extraídos con el límite de Ecuador
ecuador_ad2<-cbind(ecuador_ad2, ecu_resumen_98)
ecuador_ad2<-cbind(ecuador_ad2, ecu_resumen_14)
# Combinar las dos tables 98 y 14
combined_data <- cbind(ecu_resumen_98, ecu_resumen_14)
# preparar una versión para crear una tabla
ecu_tabla<-cbind(Area=ecuador_ad2$shapeName,
                         round(combined_data))
# crear la tabla con solo las primeras filas y 25 columnas
kableExtra::kbl(ecu_tabla[,c(1:25)]) |>
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed"))
Area E98_January E98_February E98_March E98_April E98_May E98_June E98_July E98_August E98_September E98_October E98_November E98_December E14_January E14_February E14_March E14_April E14_May E14_June E14_July E14_August E14_September E14_October E14_November E14_December
24 de Mayo 60 60 58 15 60 60 60 54 59 58 60 -3 2 -4 -4 -2 3 3 3 3 3 2 2 -4
Aguarico -6 -3 8 -4 37 -2 5 -2 -5 -8 -2 -11 0 3 3 1 3 -3 2 -3 2 2 3 -3
Alausi 4 21 10 3 60 19 13 23 4 -2 -2 -4 -3 2 -2 -2 6 -2 -3 -3 -3 -3 2 2
Alfredo Baquerizo Moreno 55 49 59 13 60 60 60 25 6 3 9 2 1 -3 -3 -2 6 3 2 2 -2 -3 2 -14
Ambato 14 25 26 7 60 60 60 14 3 -3 4 2 -13 2 -3 -3 4 1 -3 -3 -4 -3 3 -3
Antonio Ante 4 2 51 9 60 60 53 6 2 -4 52 -2 -3 2 2 -7 5 -2 4 -2 -4 -2 2 -2
Arajuno -7 -1 5 -5 58 -2 7 2 -6 -6 -2 -16 -2 3 2 1 4 -3 2 -1 -2 3 3 -3
Archidona -1 0 13 0 60 20 16 4 -2 -5 3 -5 -5 2 2 -3 4 -1 3 -2 -4 0 3 -3
Arenillas 6 7 5 18 8 7 8 7 7 8 7 6 4 1 -3 -3 -3 -2 -2 -2 -2 -2 -2 -2
Atacames 31 5 60 60 60 60 60 25 4 -2 25 3 -3 -3 -2 -8 7 3 5 2 -2 2 3 -2
Atahualpa 3 10 6 12 6 7 10 6 4 4 0 -1 3 1 1 0 2 3 3 2 2 4 3 3
Azogues 3 15 6 3 58 10 8 20 5 2 0 -4 0 2 2 -2 6 1 -3 -3 -3 -2 2 3
Baba 53 57 60 12 60 60 58 49 18 9 12 5 -3 -3 -3 -2 5 3 3 2 -2 -3 3 -8
Babahoyo 49 49 57 10 60 59 59 29 7 4 10 -1 -1 -3 -3 -2 5 3 2 2 -3 -3 1 -23
Balao 8 44 45 14 48 27 39 8 3 3 2 -7 4 -2 -3 -3 2 2 4 -2 -3 4 4 -3
Balsas 4 8 4 11 5 7 10 8 6 6 4 3 4 1 1 -2 2 2 2 1 1 3 3 3
Balzar 55 55 60 20 60 60 60 54 26 16 26 5 -4 -3 -4 -2 4 3 3 3 2 1 2 -6
Baños de Agua Santa 7 15 18 4 60 54 46 13 3 -4 3 -2 -5 2 0 -3 4 2 0 -2 -3 -1 4 -3
Biblian 4 21 9 4 60 14 9 27 6 2 2 -4 -2 2 1 -2 6 -1 -3 -3 -3 -3 2 3
Bolivar 3 2 47 6 60 60 37 6 -2 -4 44 -3 -3 2 2 -5 5 -2 3 -2 -4 -2 1 2
Bolivar 57 39 60 19 60 60 60 23 10 5 50 1 -3 -3 -4 -2 4 3 4 3 1 2 2 -26
Buena Fe 40 21 52 26 60 60 60 17 8 4 14 9 -6 -2 -5 -2 4 3 3 3 2 2 4 -2
Caluma 33 56 49 9 60 60 60 43 9 4 8 3 -7 1 -3 -3 5 2 -1 -1 -3 -3 2 -3
Calvas 3 2 -1 3 3 4 4 4 4 6 4 3 6 3 3 2 4 4 4 4 4 5 3 3
Camilo Ponce Enriquez 7 37 29 11 35 20 29 15 9 7 4 -3 3 0 0 -1 3 3 4 0 0 3 4 -2
Cañar 7 32 13 5 60 22 15 34 8 3 2 -3 -3 2 -1 -2 5 1 -1 -2 -3 -3 1 1
Carlos Julio Arosemena Tola 3 3 13 2 60 22 17 6 -1 -5 3 -3 -4 2 2 -2 4 -1 2 1 -4 1 4 -3
Cascales -3 -3 23 0 60 28 13 4 -3 -4 6 -3 -4 3 3 -3 4 -3 3 -2 -4 -2 2 0
Catamayo 4 6 3 7 4 6 10 9 7 8 6 3 7 3 4 2 5 5 4 3 3 4 3 3
Cayambe 4 1 37 5 60 58 32 4 -2 -5 18 -2 -4 2 2 -4 5 -2 3 0 -4 -2 2 -3
Celica 5 4 2 1 3 4 5 5 5 7 5 4 4 0 0 -1 0 0 -1 -1 -2 -1 -2 -1
Centinela del Condor 4 5 3 4 7 9 7 11 5 5 5 2 11 4 6 3 7 5 3 -2 -3 0 0 2
Cevallos 6 21 20 4 60 60 60 15 3 -3 2 -2 -9 2 -3 -3 4 2 -3 -3 -3 -3 3 -3
Chaguarpamba 5 9 4 9 5 7 12 11 7 6 5 3 5 3 3 -2 4 5 4 4 4 6 3 3
Chambo 4 21 12 2 60 37 24 18 3 -3 0 -4 -5 2 -2 -3 5 0 -3 -3 -3 -3 3 -1
Chilla 4 19 10 11 10 10 17 14 8 6 3 -2 3 2 3 2 5 5 5 4 3 4 3 3
Chillanes 21 60 31 7 60 53 44 54 13 4 4 0 -5 2 -3 -2 5 2 0 -2 -2 -3 2 -3
Chimbo 27 59 44 9 60 60 60 49 10 4 5 3 -8 2 -3 -3 5 2 -2 -2 -2 -3 2 -3
Chinchipe 2 -3 -3 -2 0 2 0 -2 9 8 5 2 6 2 6 3 4 4 2 2 3 4 2 -3
Chone 50 26 60 29 60 60 60 13 5 3 47 1 -3 -3 -4 -3 4 3 4 3 0 3 3 -12
Chordeleg 3 10 5 3 43 8 7 16 4 1 0 -4 2 2 3 -2 6 1 -2 -3 -3 -2 2 3
Chunchi 5 29 12 4 60 26 15 31 6 -2 -2 -4 -3 2 -2 -2 6 -2 -3 -3 -3 -3 2 2
Colimes 60 54 60 17 60 60 60 55 52 46 54 3 -2 -3 -3 -2 4 3 3 3 2 0 1 -8
Colta 8 41 20 5 60 52 40 32 5 -2 2 -2 -6 2 -3 -3 5 2 -3 -3 -3 -3 3 -2
Cotacachi 12 3 49 28 60 60 60 12 4 -3 59 -1 -3 2 2 -7 6 1 4 2 -3 2 2 -3
Crnel. Marcelino Maridueña 31 51 38 8 60 47 57 33 6 2 3 -6 2 -2 -3 -1 5 3 3 -2 -3 -1 2 -6
Cuenca 7 27 12 6 52 16 14 28 9 4 3 -3 1 2 2 -2 5 3 2 -1 -3 -2 2 2
Cumanda 17 58 24 6 60 45 31 52 14 5 4 -2 -4 2 -2 -2 5 2 2 -2 -2 -3 2 -2
Cuyabeno -8 -4 12 -5 14 -2 4 -2 -5 -6 -2 -11 -1 3 3 -2 3 -3 -2 -4 1 0 3 -5
Daule 58 57 55 12 60 60 60 28 14 12 15 -5 2 -4 -3 -2 5 2 3 2 -2 -2 2 -27
Deleg 4 17 7 4 56 12 8 23 6 3 2 -4 0 2 2 -2 6 1 -2 -3 -3 -3 2 3
Duran 40 54 54 12 60 58 60 27 7 2 6 1 0 -2 -3 -2 5 3 3 1 -2 -1 3 -4
Echeandia 39 49 56 11 60 60 60 31 7 4 15 3 -6 0 -3 -3 5 3 2 1 -3 -3 3 -3
El Carmen 42 14 55 33 60 60 60 13 6 2 30 6 -4 -2 -4 -3 5 3 3 3 0 3 4 -5
El Chaco -1 -1 18 1 60 24 16 3 -3 -4 5 -5 -4 2 2 -3 4 -2 3 -2 -4 -1 2 -3
El Guabo 5 30 22 16 15 12 22 18 13 10 5 3 3 2 3 -1 4 4 3 3 3 3 3 2
El Pan 3 7 4 2 43 7 5 14 4 -2 -2 -5 2 2 3 -2 7 -2 -3 -3 -3 -2 2 4
El Pangui 4 4 3 3 14 9 6 12 4 4 4 0 7 3 5 2 7 4 0 -3 -3 -3 0 3
El Piedrero 18 50 27 7 60 36 41 29 6 3 3 -8 1 -1 -2 1 5 3 3 -2 -3 1 1 -9
El Tambo 4 24 10 4 60 20 12 28 5 0 0 -4 -3 2 0 -2 6 -2 -3 -3 -3 -3 2 3
El Triunfo 21 54 30 7 60 40 46 28 6 2 3 -9 1 -2 -3 0 5 3 3 -2 -3 1 2 -9
Eloy Alfaro 11 3 37 52 60 60 60 13 4 -2 60 -1 -3 0 2 -9 6 3 5 1 -3 3 3 -1
Empalme 43 33 59 24 60 60 60 25 10 5 21 8 -5 -2 -5 -2 5 3 3 3 2 2 3 -4
Esmeraldas 29 5 56 60 60 60 60 20 4 -3 37 3 -3 -2 1 -5 7 3 5 2 -2 2 3 1
Espejo 3 2 48 10 60 60 51 7 0 -3 59 -3 -3 2 2 -6 5 -1 4 -2 -4 1 -1 3
Espindola -1 -1 -2 2 0 3 3 0 3 5 4 3 7 3 3 2 4 5 4 4 4 5 3 2
Flavio Alfaro 43 15 60 33 60 60 60 12 5 2 38 3 -4 -3 -4 -3 5 3 4 3 -1 3 3 -12
Giron 4 15 6 5 34 12 10 20 6 4 3 -3 2 2 3 -2 5 3 2 -2 -3 -2 3 3
Gnral. Antonio Elizalde 22 58 28 7 60 48 39 51 13 4 4 -2 -3 1 -2 -2 5 2 2 -2 -2 -3 2 -2
Gonzalo Pizarro -1 -2 26 2 60 40 17 4 -3 -4 8 -3 -4 2 2 -3 4 -2 3 -2 -4 -2 2 0
Gonzanama 4 3 2 5 4 5 7 6 5 8 6 4 8 3 4 3 5 6 5 4 3 4 3 3
Guachapala 3 8 4 3 45 7 6 14 4 -1 -2 -5 2 2 3 -2 7 -2 -3 -3 -3 -2 2 3
Gualaceo 3 12 5 3 49 9 7 18 5 1 0 -4 2 2 2 -2 6 1 -2 -3 -3 -2 2 3
Gualaquiza 4 7 4 3 23 8 6 13 4 3 2 -3 5 3 4 0 7 3 -1 -3 -3 -2 2 3
Guamote 4 20 10 3 60 23 15 19 3 -3 -2 -4 -4 2 -2 -3 5 -2 -4 -3 -3 -3 3 2
Guano 7 28 17 4 60 56 43 20 3 -3 2 -2 -7 2 -3 -3 5 2 -3 -3 -3 -3 3 -2
Guaranda 23 48 41 9 60 60 60 33 7 1 5 2 -11 2 -4 -3 4 2 0 -2 -3 -3 3 -3
Guayaquil 29 51 51 21 46 41 51 32 26 10 8 -1 3 -3 -3 -3 0 0 3 1 0 1 3 -4
Huamboya 0 6 7 0 59 14 11 10 1 -4 -1 -6 -3 3 2 -2 5 -2 0 2 -3 3 3 -2
Huaquillas 5 11 10 28 8 8 8 6 5 5 5 4 4 -2 -3 -4 -3 -3 -2 -2 -2 -2 -2 -2
Ibarra 4 2 50 12 60 60 54 7 2 -3 53 -3 -3 2 2 -7 6 -1 4 -2 -4 0 1 -1
Isabela -1 11 22 11 7 3 -1 -2 -4 -5 -8 -13 4 0 -1 -4 -3 -2 -1 -4 -6 -12 -20 -40
Isidro Ayora 54 60 55 17 59 56 60 59 60 43 26 0 3 -3 -3 -3 4 3 3 3 2 1 3 -5
Jama 57 27 60 39 60 60 60 21 7 4 29 3 -3 -4 -4 -3 4 3 4 3 2 3 3 -4
Jaramijo 60 60 60 18 60 60 60 55 45 47 60 3 2 -3 -3 -3 2 2 3 3 3 3 2 -3
Jipijapa 60 60 60 19 60 60 60 58 57 56 59 2 3 -3 -3 -3 2 2 2 2 2 2 2 -3
Junin 60 43 60 20 60 60 60 27 12 7 47 2 -3 -3 -5 -2 3 3 4 3 3 3 2 -12
La Concordia 42 8 47 44 60 60 60 20 7 3 22 7 -4 -2 -4 -2 5 3 4 3 2 2 5 -3
La Joya de los Sachas -5 -2 13 -3 60 5 10 3 -4 -4 1 -7 -3 3 3 -3 4 -3 1 -3 -3 2 3 -4
La Libertad NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
La Mana 36 20 46 20 60 60 60 16 7 3 8 7 -9 0 -5 -2 4 3 3 2 0 -2 4 -3
La Troncal 21 53 31 7 60 38 49 23 4 3 2 -10 2 -3 -2 2 5 3 3 -2 -3 2 2 -11
Lago Agrio -5 -3 21 -3 54 11 9 3 -3 -4 4 -5 -3 3 3 -3 4 -3 0 -3 -3 -2 3 0
Las Golondrinas 24 3 43 55 60 60 60 22 5 -2 55 3 -4 2 -1 -5 6 3 5 2 -2 2 3 -3
Las Lajas 6 8 4 9 5 7 10 9 8 8 7 6 5 2 2 -2 3 3 3 2 2 2 2 2
Las Naves 42 46 58 13 60 60 60 28 7 4 17 4 -5 -1 -3 -2 5 3 2 2 -2 -3 3 -3
Latacunga 18 8 34 8 60 60 54 9 3 -3 8 3 -11 2 -3 -2 4 3 0 -2 -3 -3 4 -3
Limon Indanza 3 6 4 3 36 7 6 13 4 0 -1 -4 4 3 4 -1 7 2 -3 -3 -3 -2 3 3
Logroño -2 4 3 -3 41 4 5 10 1 -4 -3 -8 2 3 3 -1 6 -2 -3 -2 -3 1 3 2
Loja 3 3 2 4 4 5 8 8 6 8 6 2 9 4 5 2 6 6 4 3 1 3 2 2
Lomas de Sargentillo 60 60 60 18 60 60 60 59 60 60 39 5 2 -3 -3 -2 4 3 3 3 2 2 3 -4
Loreto -4 -1 10 -3 60 5 10 3 -3 -4 2 -8 -3 3 2 -2 4 -2 3 -2 -4 2 3 -3
Macara 4 3 -1 3 3 3 4 3 4 4 3 3 6 2 3 2 3 4 3 3 2 3 -2 -1
Machala 5 14 15 24 10 10 16 15 8 7 5 3 4 2 2 -2 4 3 2 2 2 2 3 3
Manta 60 60 60 18 60 60 60 55 45 47 60 3 2 -3 -3 -3 2 2 3 3 3 3 2 -3
Marcabeli 4 8 4 10 5 7 10 9 6 6 4 3 5 2 2 -2 3 3 3 2 2 3 3 3
Mejia 21 6 38 12 60 60 55 10 4 -3 13 4 -11 2 -3 -2 4 3 2 0 -3 -3 3 -3
Mera 6 13 16 3 60 51 36 12 2 -4 2 -3 -4 3 2 -2 4 1 2 -1 -3 2 4 -3
Milagro 41 52 45 9 60 55 60 39 7 2 4 -2 3 -2 -3 -2 5 2 3 -2 -3 -2 2 -2
Mira 4 2 50 16 60 60 60 7 2 -3 60 -3 -2 2 2 -7 6 0 4 -2 -3 2 0 2
Mocache 43 50 60 17 60 60 60 40 13 6 12 7 -7 -2 -4 -2 5 3 3 3 2 -2 3 -3
Mocha 7 22 20 4 60 60 60 15 3 -3 2 -2 -9 2 -3 -3 4 2 -3 -3 -3 -3 3 -3
Montalvo 48 48 59 9 60 60 60 26 5 4 9 -7 -2 -3 -3 -2 6 2 2 2 -3 -3 -1 -48
Montecristi 60 60 60 18 60 60 60 55 45 47 60 3 2 -3 -3 -3 2 2 3 3 3 3 2 -3
Montufar 2 2 47 6 60 60 34 6 -2 -4 50 -3 -3 2 2 -5 5 -2 3 -2 -4 -2 -1 3
Morona -1 4 4 -3 53 5 5 10 0 -4 -3 -9 -1 3 2 -1 5 -2 -2 0 -3 1 3 0
Muisne 33 6 60 60 60 60 60 26 4 -2 26 4 -3 -3 -2 -7 6 3 5 3 -2 2 3 -4
Nabon 4 11 5 5 18 10 10 16 6 4 3 -2 3 2 3 -2 6 3 2 -2 -3 1 3 3
Nangaritza -1 -1 0 2 3 5 5 5 4 8 6 2 18 6 9 3 7 6 3 -1 -2 -1 -1 -2
Naranjal 21 58 42 11 59 39 53 26 7 2 3 -6 3 0 -2 -2 4 3 4 -2 -3 2 3 -7
Naranjito 36 50 41 8 60 51 60 40 7 3 3 -2 2 -2 -3 -2 5 2 3 -2 -3 -2 2 -2
Nobol 56 60 57 17 60 57 60 59 59 47 30 1 3 -3 -3 -3 4 3 3 3 2 1 3 -5
Olmedo 5 9 4 9 5 7 12 11 7 6 5 2 5 3 3 -2 4 5 5 4 4 6 4 3
Olmedo 60 58 59 15 60 60 60 59 60 46 60 -2 -1 -3 -3 -2 4 3 3 3 2 -1 -1 -13
Oña 4 10 4 5 13 10 9 14 5 4 3 -2 4 3 3 -2 6 4 2 -2 -3 2 3 3
Orellana -7 -2 9 -4 60 1 8 1 -4 -5 -1 -13 -3 3 3 -2 3 -3 2 -3 -3 2 3 -3
Otavalo 8 3 50 13 60 60 52 7 1 -4 46 -1 -3 2 2 -5 5 0 4 1 -4 -1 2 -3
Pablo Sexto 1 5 6 -2 60 10 9 11 2 -4 -2 -5 -4 2 1 -3 5 -2 -2 -1 -3 -1 3 -2
Pajan 60 60 56 17 59 59 59 52 58 57 51 -1 2 -3 -3 -3 3 2 3 2 2 1 2 -4
Palanda -2 -3 -3 -2 -1 2 2 1 5 11 4 1 11 4 7 4 7 7 4 3 3 4 2 -2
Palenque 49 58 60 19 60 60 60 57 19 10 15 7 -6 -2 -4 -2 5 3 3 3 2 -1 3 -3
Palestina 60 48 60 19 60 60 60 24 9 5 24 7 2 -2 -3 -2 6 3 3 3 2 -2 4 -3
Pallatanga 7 40 18 5 60 44 29 36 6 -1 0 -3 -5 2 -3 -3 5 0 -3 -3 -3 -3 2 0
Palora 4 10 12 2 60 34 23 11 2 -4 1 -3 -3 3 1 -2 4 0 1 -1 -3 2 4 -3
Paltas 4 6 3 7 4 6 8 8 6 7 5 3 6 3 3 1 4 5 4 4 3 4 3 3
Pangua 32 41 51 15 60 60 60 29 9 3 7 5 -10 1 -4 -2 5 3 2 2 -1 -2 4 -3
Paquisha 5 4 3 4 8 9 7 10 5 5 5 2 10 4 6 2 6 4 2 -3 -3 -2 -2 1
Pasaje 4 18 12 17 11 10 17 16 9 7 4 1 3 2 3 0 5 4 4 3 3 3 3 3
Pastaza -5 2 5 -4 58 1 7 4 -5 -6 -3 -14 0 3 2 2 4 -3 2 -1 -2 3 3 -1
Patate 6 15 18 4 60 54 48 13 3 -4 3 -2 -7 2 -2 -3 4 2 -2 -2 -3 -3 4 -3
Paute 3 14 6 4 55 10 8 20 5 2 0 -4 2 2 2 -2 6 1 -2 -3 -3 -2 2 3
Pedernales 41 11 59 50 60 60 60 23 6 2 32 4 -4 -3 -3 -4 5 3 4 3 0 3 3 -8
Pedro Carbo 58 59 55 17 57 57 60 56 60 59 27 0 2 -3 -3 -3 3 3 3 3 2 0 2 -3
Pedro Moncayo 8 3 48 10 60 60 52 7 1 -4 38 0 -4 2 2 -4 5 1 3 1 -4 -2 2 -3
Pedro Vicente Maldonado 34 4 41 49 60 60 60 23 6 -2 38 4 -5 -1 -3 -3 5 3 4 3 -2 1 4 -2
Penipe 3 11 10 0 60 23 21 12 3 -4 -1 -4 -5 2 -2 -3 5 -1 -3 -2 -3 -3 4 -2
Pichincha 50 36 60 22 60 60 60 23 9 5 43 2 -4 -3 -4 -2 5 3 3 3 0 2 3 -24
Pimampiro 3 1 43 5 60 60 32 5 -2 -4 30 -3 -3 2 2 -5 5 -2 3 -2 -4 -2 2 1
Piñas 3 8 5 12 6 7 9 7 5 5 2 1 4 0 0 -1 1 1 1 1 1 4 3 3
Pindal 6 4 2 -2 2 3 4 4 5 6 5 4 2 -3 -3 -3 -3 -4 -4 -5 -5 -4 -5 -4
Playas 17 56 52 29 23 26 39 27 25 26 21 7 4 -3 -3 -3 -2 -2 2 2 2 2 2 -3
Portovelo 4 9 5 8 5 7 12 10 7 6 3 1 5 3 3 0 5 5 4 3 2 4 3 3
Portoviejo 60 54 60 17 60 60 60 47 35 28 57 -1 -2 -4 -4 -2 3 3 3 3 3 3 2 -16
Pucara 5 25 11 8 18 12 18 20 13 8 5 0 3 2 3 1 5 4 4 3 1 2 3 1
Puebloviejo 49 59 60 14 60 60 57 59 21 11 13 5 -5 -3 -3 -2 5 3 3 2 -1 -3 3 -8
Puerto Lopez 60 60 60 19 60 60 60 58 57 55 60 3 3 -3 -3 -3 2 2 3 3 2 2 2 -3
Puerto Quito 37 5 42 55 60 60 60 26 7 0 32 4 -4 -1 -3 -3 5 3 4 3 0 2 4 -2
Pujili 27 24 41 12 60 60 60 15 5 -1 6 4 -11 2 -4 -2 4 3 0 -1 -3 -3 4 -3
Putumayo -8 -4 17 -6 22 3 5 1 -4 -4 2 -7 -2 3 3 -3 3 -3 -3 -3 0 -2 3 -2
Puyango 5 7 4 5 4 5 8 8 6 7 6 5 4 1 1 -2 2 2 2 1 0 1 1 1
Quero 6 21 20 4 60 60 60 15 3 -3 2 -2 -9 2 -3 -3 4 2 -3 -3 -3 -3 3 -3
Quevedo 39 42 57 21 60 60 60 33 11 5 10 8 -7 -2 -5 -2 5 3 3 3 2 -2 4 -3
Quijos 1 1 18 1 60 30 18 4 -3 -5 5 -4 -5 2 2 -3 4 -2 3 -2 -4 -1 3 -3
Quilanga 3 2 -2 3 3 3 4 4 3 6 4 3 9 3 4 3 5 6 5 4 4 5 3 2
Quininde 32 5 50 58 60 60 60 25 6 0 40 4 -3 -2 -2 -4 6 3 5 3 -1 3 4 -3
Quinsaloma 39 45 58 20 60 60 60 33 10 5 11 7 -7 -2 -4 -2 5 3 3 3 1 -2 4 -3
Quito 16 3 44 16 60 60 58 10 3 -3 33 2 -6 2 -1 -3 5 2 3 1 -3 -2 3 -3
Riobamba 6 27 15 3 60 41 32 20 4 -3 0 -3 -6 2 -2 -3 5 0 -3 -3 -3 -3 3 -1
Rioverde 16 3 45 60 60 60 60 16 4 -3 55 2 -3 -2 2 -10 7 3 5 2 1 3 3 1
Rocafuerte 60 50 60 20 60 60 60 28 14 8 60 -2 -3 -4 -4 -2 3 3 4 4 3 3 3 -4
Rumiñahui 14 4 38 9 60 60 60 8 3 -3 15 3 -6 2 -2 -3 4 3 3 -2 -3 -3 3 -3
Salcedo 14 17 28 7 60 57 53 12 3 -3 5 2 -11 2 -3 -2 4 3 -2 -2 -3 -3 4 -3
Salinas NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
Salitre 59 54 57 13 60 60 60 22 6 3 14 -3 3 -3 -3 -2 6 3 2 3 -2 -3 3 -22
Samborondon 53 53 56 12 60 60 60 25 6 3 9 -2 1 -3 -3 -2 5 3 3 2 -3 -3 3 -17
San Cristobal NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
San Fernando 5 18 7 5 30 12 13 22 9 5 4 -2 3 2 3 -2 5 3 3 1 -3 -2 3 3
San Jacinto de Yaguachi 40 53 47 10 60 54 60 32 6 2 4 -3 2 -2 -3 -2 5 2 3 -1 -3 -2 2 -4
San Juan Bosco 3 5 3 2 24 8 6 12 4 2 1 -4 5 3 4 0 7 3 -1 -3 -3 -2 2 3
San Lorenzo 5 2 39 38 60 60 60 8 3 -2 60 -3 -2 2 3 -10 6 3 4 -1 -3 3 1 0
San Miguel 22 59 38 8 60 60 57 52 10 3 4 2 -7 2 -3 -3 5 2 -2 -2 -3 -3 2 -3
San Miguel de Los Bancos 32 5 45 31 60 60 60 17 6 -2 34 4 -6 2 -3 -3 5 3 3 2 -1 -1 3 -3
San Miguel de Urcuqui 5 2 51 17 60 60 59 8 3 -3 58 -3 -3 2 2 -8 6 0 4 -2 -3 1 0 -2
San Pedro de Huaca -2 -2 37 3 60 60 20 5 -3 -4 24 -3 -3 2 2 -4 4 -3 3 -2 -4 -2 2 4
San Pedro de Pelileo 7 21 20 4 60 60 60 15 3 -3 2 -2 -9 2 -3 -3 4 2 -3 -3 -3 -3 3 -3
San Vicente 60 40 60 27 60 60 60 21 9 5 37 2 -3 -4 -5 -2 3 3 4 4 3 3 3 -3
Santa Ana 60 52 60 16 60 60 60 51 38 29 53 0 -2 -3 -4 -2 4 3 3 3 3 2 2 -33
Santa Clara -1 2 8 -1 60 9 11 6 -3 -5 0 -7 -3 3 2 -2 4 -2 2 2 -4 3 3 -3
Santa Cruz 2 14 30 15 8 5 3 2 -2 -3 -3 -3 5 3 2 -3 -3 -3 -3 -3 -4 -5 -5 -7
Santa Elena 42 60 59 29 48 53 59 46 42 34 30 4 4 -3 -3 -3 0 1 1 1 1 1 2 -3
Santa Isabel 5 18 7 6 21 12 15 21 11 6 4 0 3 3 3 -2 5 4 3 2 -2 0 3 3
Santa Lucia 60 54 59 18 60 60 60 48 43 42 45 4 0 -3 -3 -2 5 3 3 3 2 1 3 -7
Santa Rosa 4 10 9 22 8 8 10 9 6 6 4 4 4 -1 -1 -3 -1 0 -1 -1 -1 1 1 1
Santiago 2 9 5 2 48 8 6 14 4 -1 -2 -5 2 3 3 -2 7 0 -3 -3 -3 -2 3 3
Santiago de Pillaro 7 9 19 4 60 48 36 10 3 -4 3 -1 -9 2 -3 -3 4 2 -2 -2 -3 -3 4 -3
Santo Domingo 40 8 44 33 60 60 60 16 7 2 16 7 -6 -1 -4 -2 5 3 3 3 1 1 4 -2
Saquisili 25 14 39 10 60 60 60 12 4 -3 7 4 -11 2 -4 -2 4 3 0 -2 -3 -3 4 -3
Saraguro 5 11 5 6 8 9 12 14 8 6 5 1 5 3 3 -2 5 4 3 1 -2 2 3 3
Sevilla de Oro 3 12 5 3 54 9 7 17 4 1 -2 -4 2 2 2 -2 6 1 -3 -3 -3 -2 2 3
Shushufindi -5 -3 15 -3 57 5 9 2 -4 -5 2 -7 -3 3 3 -2 3 -3 2 -3 -2 1 3 -3
Sigchos 34 15 43 18 60 60 60 14 6 -2 9 5 -10 2 -4 -2 4 3 2 1 -2 -2 4 -3
Sigsig 3 9 4 3 40 8 7 15 4 3 2 -4 3 2 3 -2 6 2 -2 -3 -3 -2 2 3
Simon Bolivar 40 49 46 8 60 53 60 37 6 3 4 -4 1 -2 -3 -2 5 2 3 -1 -3 -3 1 -15
Sozoranga 3 3 0 4 3 4 5 4 4 6 4 3 6 3 3 2 3 4 4 3 3 3 1 1
Sucre 59 40 60 25 60 60 60 20 9 5 57 -1 -3 -4 -4 -3 4 3 4 4 1 3 3 -6
Sucua 0 6 4 -2 47 6 5 12 3 -3 -3 -7 0 3 2 -2 6 -1 -3 -2 -3 1 3 3
Sucumbios -2 -2 32 3 60 54 19 5 -3 -4 14 -3 -4 2 2 -4 4 -3 3 -2 -4 -2 2 3
Suscal 12 45 18 6 60 29 21 43 12 4 3 -2 -3 2 -2 -2 5 3 2 -2 -3 -3 -2 -3
Taisha -2 4 5 -2 60 5 8 7 -4 -5 -3 -12 2 3 2 2 5 -3 -2 0 -3 3 3 -2
Tena 0 2 12 0 60 19 16 5 -2 -4 2 -5 -4 3 2 -2 4 -1 2 -1 -4 1 3 -3
Tisaleo 6 21 20 4 60 60 60 15 3 -3 2 -2 -9 2 -3 -3 4 2 -3 -3 -3 -3 3 -3
Tiwintza -1 4 3 -1 39 5 6 9 0 -3 -3 -8 3 3 3 1 6 -1 -3 -3 -3 0 4 2
Tosagua 60 45 60 21 60 60 60 25 11 7 56 0 -3 -4 -4 -2 3 3 4 4 2 3 3 -7
Tulcan 3 2 42 16 60 60 52 6 1 -3 56 -3 -2 2 2 -7 6 1 4 -2 -3 2 -1 4
Urdaneta 50 46 60 11 60 60 60 24 6 4 16 3 -3 -2 -3 -3 6 3 2 3 -3 -3 3 -4
Valencia 38 21 48 24 60 60 60 17 8 3 9 8 -7 -2 -5 -2 4 3 3 2 2 -1 4 -2
Ventanas 45 45 60 13 60 60 60 28 7 5 20 4 -4 -2 -3 -2 5 3 2 3 -2 -3 3 -3
Vinces 54 55 60 18 60 60 59 45 16 8 18 6 -2 -2 -3 -2 5 3 3 3 1 -2 3 -4
Yacuambi 5 8 4 4 12 9 8 13 5 4 4 -2 6 3 4 0 6 4 2 -2 -3 1 2 3
Yantzaza 5 5 3 4 10 9 7 11 5 5 4 1 8 4 5 2 7 4 2 -3 -3 -1 0 2
Zamora 3 4 3 5 5 7 8 10 6 6 6 3 12 5 6 3 6 5 3 1 -1 2 2 2
Zapotillo 4 3 4 2 3 3 3 1 1 3 1 0 1 0 0 2 0 -1 -2 -2 -2 0 -3 3
Zaruma 4 13 7 9 7 8 13 11 7 5 2 -1 4 2 3 1 5 5 4 3 1 4 3 3

Esto confirma la distribución generalizada de valores de excedente elevados en nuestros mapas exploratorios. Actualmente, los datos están en formato ancho, lo que facilita la visualización de una serie temporal, pero la visualización programática más avanzada generalmente requiere que los datos estén en un formato largo (más sobre esto más adelante).

Verificación de conocimientos
  1. ¿Para qué se pueden utilizar las estadísticas zonales?
    1. Subcomponer datos a un período de tiempo de interés.
    2. Creando una serie de tiempo.
    3. Resumir los valores de las celdas que se encuentran dentro de un límite.
    4. Todo lo anterior.

Coropletas de Areas Administrativas

Ahora que hemos inspeccionado los datos, podemos hacer una coropleta a partir de los datos del período de retorno del excedente medio. Anteriormente demostramos mapas más complejos usando ggplot2, sf y stars, pero también puedes hacer trazados rápidos de objetos sf con la función de trazado base. De forma predeterminada, sf creará un mapa para cada columna enumerada en el conjunto de datos.

colnames(ecuador_ad2)
 [1] "shapeName"     "shapeISO"      "shapeID"       "shapeGroup"   
 [5] "shapeType"     "E98_January"   "E98_February"  "E98_March"    
 [9] "E98_April"     "E98_May"       "E98_June"      "E98_July"     
[13] "E98_August"    "E98_September" "E98_October"   "E98_November" 
[17] "E98_December"  "E14_January"   "E14_February"  "E14_March"    
[21] "E14_April"     "E14_May"       "E14_June"      "E14_July"     
[25] "E14_August"    "E14_September" "E14_October"   "E14_November" 
[29] "E14_December"  "geometry"     

En este caso solo queremos ver las medias mensuales, por lo que solo trazaremos las columnas 6 a 18. Puede realizar modificaciones simples en la paleta de colores y la posición de la leyenda, pero los títulos de mapas personalizados y los títulos de leyendas no se logran fácilmente con múltiples -Mapas en panel (varios mapas en uno) como el que se muestra a continuación.

# trazar los datos para un cheque usando solo los 12 resúmenes mensuales en las columnas 
plot(ecuador_ad2[c(6, 18)], # 15:17
     # la paleta que hemos estado usando
     pal = c(
    # -60 a -30
    '#9B0039',
    #-30 a -10
    '#FF8D43',
    # -10 a -5
    '#FFC754',
    # -5 a 0
    '#FFF4C7',
    # 0 a 5
    '#DDFFFF',
    # 5 a 10
    '#00BFFF',
    # 10 a 30
    '#0000FF',
    # 30 a 60
    '#191970'), 
    # establecer las diviciones de la paleta
    breaks = c(-60, -30, -10, -5, 0, 5, 10, 30, 60), 
    # donde colocaremos la leyenda
    key.pos = 1, 
    # qué tan gruesa queremos que sea la leyenda
    key.width = 0.3,
    # Ajustar el grueso del borde
    border=NA)

Debido a los excedentes hídricos generalizados en los datos sin procesar, los valores medios no parecen muy diferentes de la capa ráster de excedente sin procesar; sin embargo, los mapas coropléticos, también llamados mapas temáticos, pueden facilitar a los usuarios el estudio del paisaje mediante la visualización de lugares familiares. (como areas administrativas o ciudades) que se colocan a sí mismos y a sus experiencias vividas junto a los datos.

Si bien esto muestra un panorama sorprendente de excedentes hídricos generalizados, ¿cuántas personas se ven afectadas por estas inundaciones? Aunque la superficie terrestre parece bastante grande, si uno no está familiarizado con la distribución de la población y los centros urbanos de Ecuador, puede resultar difícil tener una idea del impacto humano directo. (Esto se debe en parte a que los lugares más poblados suelen estar representados por superficies de tierra más pequeñas y los lugares menos poblados suelen estar representados por grandes fronteras administrativas que contienen mucha más superficie de tierra). Normalizar un tema determinado por superficie terrestre puede ser algo que un analista quiera hacer, pero a continuación cubrimos otro enfoque.

Verificación de conocimientos
  1. Los mapas coropléticos, también conocidos como mapas temáticos, son útiles para
    1. Visualizar datos en geografías familiares como condados o estados.
    2. Encontrar direcciones de un lugar a otro
    3. Visualización de datos en unidades geográficas uniformes, como celdas de cuadrícula ráster.
    4. Cálculo de periodos de retorno.

Integración de datos de población

Gridded Population of the World (GPW) es una recopilación de datos de SEDAC que modela la distribución de la población humana global como recuentos y densidades en formato ráster (Center For International Earth Science Information Network-CIESIN-Columbia University 2018). Aprovecharemos al máximo exactextractr para integrarse a través de WSIM-GLDAS, geoBoundaries y GPW. Para comenzar, necesitamos descargar los datos de densidad de población para el año 2015 con una resolución de 15 minutos (aproximadamente 30 kilómetros cuadrados en el ecuador) de GPW. Esta versión de GPW se acerca más a nuestro período de tiempo (2014) y a la resolución de WSIM-GLDAS (0,25 grados). Aunque en muchas aplicaciones se puede optar por utilizar las capas de datos de recuento de población de GPW, como utilizamos exactextractr, podemos lograr resultados más precisos (especialmente a lo largo de las costas) utilizando la densidad de población junto con las estimaciones de superficie terrestre delexactextractr paquete.

Revisión de datos

La población cuadriculada del mundo versión 4 está disponible en múltiples métricas objetivo (por ejemplo, recuentos, densidad), períodos (2000, 2005, 2010, 2015, 2020) y resoluciones espaciales (30 segundos, 2,5 minutos, 15 minutos, 30 minutos, 60 minutos). Lea más sobre GPW en la página de inicio de la colección en SEDAC. GPW es uno de los cuatro conjuntos de datos de población cuadriculados globales disponibles en formato ráster. Estos conjuntos de datos varían en el grado en que utilizan información adicional como variables auxiliares para modelar la distribución espacial de la población a partir de las unidades administrativas (polígonos vectoriales) en las que se originan. En un artículo de Leyk y sus colegas (Leyk et al. 2019) se encuentra una revisión de estos conjuntos de datos y sus modelos subyacentes. Puede obtener más información sobre los conjuntos de datos de población cuadriculados en POPGRID.org. La idoneidad para el uso es un principio importante para determinar el mejor conjunto de datos a utilizar para un análisis específico. La pregunta que nos hacemos aquí es: ¿cuál es la exposición de la población a diferentes niveles de exceso de agua en Ecuador? — utiliza entradas espacialmente aproximadas y es para un lugar con entradas de datos de alta calidad, GPW es una buena opción para este análisis. Los usuarios con datos censales en formato vectorial (a nivel de condado o subcondado) también podrían adoptar este enfoque para esos datos. En el caso de California, los datos del censo de EE. UU. y el GPW producirán estimaciones casi idénticas porque el PGB se basa en los datos del censo.

Cargue el GPW tif desde el almacenamiento en la nube de SCHOOL s3.

# leer el archivo con terra desde nuestro depósito s3
pob_dens<- 
terra::rast("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/gpw_v4_population_density_rev11_2015_15_min.tif")

# comprobar la información básica
print(pob_dens)
class       : SpatRaster 
dimensions  : 720, 1440, 1  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : gpw_v4_population_density_rev11_2015_15_min.tif 
name        : gpw_v4_population_density_rev11_2015_15_min 

Para este ejemplo, clasificaremos la capa ráster del período de retorno del excedented WSIM-GLDAS en ocho categorías. Combinar los datos facilitará la gestión de la salida y la interpretación de los resultados.

# establece los saltos de clase por filas (desde, hasta, nueva etiqueta)
# Por ejemplo, la primera fila indica: "todos los períodos de retorno del 0 al 5 ahora se etiquetarán como 0
m <- c(-60, -30, -30,
       -30,-10, -10,
       -10,-5,-5,
      -5, 0, 0,
       0, 5, 5,
       5, 10, 10,
       10, 30, 30,
       30, 60, 60)
# convertir el vector en una matriz
rclmat <- matrix(m, ncol=3, byrow=TRUE)
# clasificar los datos
wsim_gldas_1mo_class_98 <-
  terra::classify(wsim_gldas_1mo_98, rclmat, include.lowest = TRUE)
wsim_gldas_1mo_class_14 <-
  terra::classify(wsim_gldas_1mo_14, rclmat, include.lowest = TRUE)

En nuestro ejemplo anterior, utilizamos la función “media” incorporada de exactextractr, pero también podemos pasar funciones personalizadas a exactextractr que también llevarán a cabo varias operaciones a la vez. El siguiente código podría combinarse en una única función pasada a exactextractr, pero aquí se presenta como funciones múltiples para poder seguirlo más fácilmente. Puede leer más sobre los argumentos de exactextractr en el paquete guía de ayuda. Los argumentos clave a tener en cuenta son los llamados a:

  1. weights = pob_dens: resume el período de retorno del excedente de cada celda WSIM-GLDAS con el valor de densidad de población correspondiente.
  2. coverage_area = TRUE: calcula el área correspondiente de la celda ráster WSIM-GLDAS que está cubierta por el límite de Ecuador
# ejecutar la extracción
pop_by_rp_98 <-
  exactextractr::exact_extract(wsim_gldas_1mo_class_98, ecuador_ad2, function(df) {
    df <- data.table::setDT(df)}, 
  # convertir la salida a un marco de datos
  summarize_df = TRUE, 
  # especificamos los pesos que estamos usando
  weights = pob_dens, 
  # devolver el área de cobertura (m^2)
  coverage_area = TRUE,
  # devolver el nombre del area con los datos de salida
  include_cols = 'shapeISO', 
  # no mostrar progreso
  progress = FALSE)

# ejecutar la extracción
pop_by_rp_14 <-
  exactextractr::exact_extract(wsim_gldas_1mo_class_14, ecuador_ad2, function(df) {
    df <- data.table::setDT(df)}, 
  # convertir la salida a un marco de datos
  summarize_df = TRUE, 
  # especificamos los pesos que estamos usando
  weights = pob_dens, 
  # devolver el área de cobertura (m^2)
  coverage_area = TRUE,
  # devolver el nombre del area con los datos de salida
  include_cols = 'shapeISO', 
  # no mostrar progreso
  progress = FALSE)

Esto devuelve un data.frame con una fila para cada celda ráster en la capa WSIM-GLDAS que se superpone con el límite de Ecuador. Echemos un vistazo a las primeras 6 filas.

head(pop_by_rp_98)
   shapeISO 1998-01-01 1998-02-01 1998-03-01 1998-04-01 1998-05-01 1998-06-01
     <char>      <num>      <num>      <num>      <num>      <num>      <num>
1:                  60         60         60         30         60         60
2:                  60         60         60         30         60         60
3:                  60         60         60         30         60         60
4:                  60         60         60         30         60         60
5:                  -5          0         30         -5         10          0
6:                  -5          0         30          0         10          0
   1998-07-01 1998-08-01 1998-09-01 1998-10-01 1998-11-01 1998-12-01
        <num>      <num>      <num>      <num>      <num>      <num>
1:         60         60         60         60         60          5
2:         60         60         60         60         60          0
3:         60         60         60         60         60          0
4:         60         60         60         60         60          0
5:          5          0          0         -5          0        -10
6:          5          0          0         -5          0        -10
        weight coverage_area
         <num>         <num>
1: 165.4402771  7.000643e+03
2: 257.8148193  3.607563e+07
3:  62.8166161  4.526563e+08
4:  43.7354660  4.049716e+07
5:   0.9948775  1.421737e+08
6:   0.2794557  7.428570e+07
  • shapeISO: La etiqueta del límite del polígono donde se encuentra la celda. Esto se transmitió desde el límite geojson de Ecuador como se especifica en el argumento include_cols = 'shapeISO'. En este caso, no es muy útil porque utilizamos el límite estatal de Ecuador, pero si pasamos el límite ADM2 con los condados, proporcionaría el nombre del condado donde se encuentra la celda.

  • 2014-01-01 a 2014-12-01: Las siguientes 12 columnas enumeran el valor de clasificación del período de retorno del déficit para la celda en cada uno de los 12 meses correspondientes a la dimensión temporal de la capa ráster wsim_gldas_1mo.

  • weights: la columna weight enumera el valor de densidad de población correspondiente (personas por km^2) para esa celda WSIM-GLDAS. Las capas ráster WSIM-GLDAS y GPW tienen la misma proyección y resolución. Por lo tanto, debido a que están perfectamente alineadas, cada celda del período de retorno de WSIM tiene un peso de población de GPW correspondiente que se puede superponer exactamente sobre ella.

  • coverage_area: el área total (m^2) de la celda WSIM-GLDAS que está cubierta por la capa límite de California. Dada la superficie total de la celda del WSIM que está cubierta y el peso de personas del GPW por unidad de área, podemos calcular el número de personas que se estima que vivirán dentro de esta celda durante este período de retorno del déficit del WSIM.

Necesitaremos realizar algunos pasos de procesamiento más para preparar este “marco de datos” para una visualización de series de tiempo que integre todos los datos. Usaremos la función melt para transformar los datos de formato ancho a formato largo para producir una visualización en ggplot2. Específicamente, necesitamos usar melt para convertir las columnas de 12 meses (2014-01-01 a 2014-12-01) en 2 columnas nuevas: 1) especificando el valor del período de retorno del excedente de WSIM-GLDAS y 2 ) el mes de donde vino.

Coding Review

Revisión de codificación

La conversión de datos de formatos ancho a largo o de largo a ancho es un componente clave para el procesamiento de datos; sin embargo, puede resultar confuso. Para leer más sobre cómo fundir/pivotar más tiempo (de ancho a largo) y cómo fundir/pivotar más ancho (de largo a ancho), consulte la viñeta data.table Remodelación eficiente usando data.tables y dplyr pivot_longer y pivot_wider páginas de referencia.

# convertir el conjunto de datos de ancho a largo (melt)
pop_by_rp_98 <-
  data.table::melt(
    # datos que estamos derritiendo
    pop_by_rp_98,
    # las columnas que deben permanecer como columnas/ids
    id.vars = c("shapeISO", "coverage_area", "weight"),
    # nombre de las columnas del mes que vamos a fusionar en una sola columna
    variable.name = "month",
    # nombre de los valores que se convierten en una sola columna
    value.name = "return_period")

# comprobar las primeras 5 filas de datos derretidos
head(pop_by_rp_98)
   shapeISO coverage_area      weight      month return_period
     <char>         <num>       <num>     <fctr>         <num>
1:           7.000643e+03 165.4402771 1998-01-01            60
2:           3.607563e+07 257.8148193 1998-01-01            60
3:           4.526563e+08  62.8166161 1998-01-01            60
4:           4.049716e+07  43.7354660 1998-01-01            60
5:           1.421737e+08   0.9948775 1998-01-01            -5
6:           7.428570e+07   0.2794557 1998-01-01            -5
# convertir el conjunto de datos de ancho a largo (melt)
pop_by_rp_14 <-
  data.table::melt(
    # datos que estamos derritiendo
    pop_by_rp_14,
    # las columnas que deben permanecer como columnas/ids
    id.vars = c("shapeISO", "coverage_area", "weight"),
    # nombre de las columnas del mes que vamos a fusionar en una sola columna
    variable.name = "month",
    # nombre de los valores que se convierten en una sola columna
    value.name = "return_period")

# comprobar las primeras 5 filas de datos derretidos
head(pop_by_rp_14)
   shapeISO coverage_area      weight      month return_period
     <char>         <num>       <num>     <fctr>         <num>
1:           7.000643e+03 165.4402771 2014-01-01             5
2:           3.607563e+07 257.8148193 2014-01-01             0
3:           4.526563e+08  62.8166161 2014-01-01             5
4:           4.049716e+07  43.7354660 2014-01-01             5
5:           1.421737e+08   0.9948775 2014-01-01             0
6:           7.428570e+07   0.2794557 2014-01-01             5

Cada fila enumera la superficie terrestre (área de cobertura) cubierta por la zona, el valor de densidad de población (peso) de la zona, el mes al que corresponde el período de retorno del déficit y la clase del período de retorno del déficit real para la zona.

A continuación, resumiremos los datos por clase de período de retorno.

# crear una nueva columna con totales de población para cada mes y combinación de período de retorno
# dividir por 1000000
pop_by_rp_98 <-
  pop_by_rp_98[, .(pop_rp = round(sum(coverage_area * weight) / 1e6)), by = .(month, return_period)]
# algunas celdas no tienen valores de período de retorno deficitario y dan como resultado NaN; simplemente elimínelas
pop_by_rp_98 <- na.omit(pop_by_rp_98)
# crear una nueva columna con la población total para ese mes (será la misma para cada mes pero será necesaria para crear una fracción de clase wsim)
pop_by_rp_98[, total_pop := sum(pop_rp), by = month]
# revisa las primeras filas
head(pop_by_rp_98)
        month return_period  pop_rp total_pop
       <fctr>         <num>   <num>     <num>
1: 1998-01-01            60 5583311  15279575
2: 1998-01-01            -5  151649  15279575
3: 1998-01-01             0  580917  15279575
4: 1998-01-01            30 4861930  15279575
5: 1998-01-01             5 2748874  15279575
6: 1998-01-01            10 1350026  15279575
# crear una nueva columna con totales de población para cada mes y combinación de período de retorno
# dividir por 1000000
pop_by_rp_14 <-
  pop_by_rp_14[, .(pop_rp = round(sum(coverage_area * weight) / 1e6)), by = .(month, return_period)]
# algunas celdas no tienen valores de período de retorno deficitario y dan como resultado NaN; simplemente elimínelas
pop_by_rp_14 <- na.omit(pop_by_rp_14)
# crear una nueva columna con la población total para ese mes (será la misma para cada mes pero será necesaria para crear una fracción de clase wsim)
pop_by_rp_14[, total_pop := sum(pop_rp), by = month]
# revisa las primeras filas
head(pop_by_rp_14)
        month return_period  pop_rp total_pop
       <fctr>         <num>   <num>     <num>
1: 2014-01-01             5 5284826  15279575
2: 2014-01-01             0 4486888  15279575
3: 2014-01-01           -10 1335779  15279575
4: 2014-01-01            -5 3626857  15279575
5: 2014-01-01            10  397448  15279575
6: 2014-01-01            30  147777  15279575

Ahora tenemos una fila para cada combinación única de mes, clase de período de retorno y población total. Podemos calcular el porcentaje de la población total representada por esta clase del período de retorno con 2 líneas más.

# calcular la fracción de esa combinación de clase mes-wsim en relación con la población total
pop_by_rp_98[, pop_frac := pop_rp / total_pop][, total_pop := NULL]
# calcular la fracción de esa combinación de clase mes-wsim en relación con la población total
pop_by_rp_14[, pop_frac := pop_rp / total_pop][, total_pop := NULL]


# revisa las primeras filas
head(pop_by_rp_14)
        month return_period  pop_rp    pop_frac
       <fctr>         <num>   <num>       <num>
1: 2014-01-01             5 5284826 0.345875196
2: 2014-01-01             0 4486888 0.293652670
3: 2014-01-01           -10 1335779 0.087422523
4: 2014-01-01            -5 3626857 0.237366353
5: 2014-01-01            10  397448 0.026011718
6: 2014-01-01            30  147777 0.009671539

Antes de trazar, haremos que las etiquetas de los meses sean más legibles para el trazado, convertiremos la clase del período de retorno WSIM-GLDAS en un factor y configuraremos la paleta de clases WSIM-GLDAS.

Revisión de codificación

Los factores son la forma más común de manejar datos categóricos en R. Aunque convertir sus variables categóricas en factores no siempre es la mejor opción, en muchos casos (especialmente al trazar con ggplot2) los beneficios superarán cualquier molestia. Para obtener más información sobre los factores y R, consulte el capítulo de Hadley Wickham sobre factores en R para ciencia de datos, segunda edición.

# ggplot es más fácil con factores
pop_by_rp_98$return_period<-
  factor(pop_by_rp_98$return_period, 
         levels = c("-60", "-30", "-10", "-5","0", "5", "10", "30", "60"))

# agregar los datos base
pop_by_rp_98[,month:=lubridate::month(pop_by_rp_98$month, label = TRUE)]

# ggplot es más fácil con factores
pop_by_rp_14$return_period<-
  factor(pop_by_rp_14$return_period, 
         levels = c("-60", "-30", "-10", "-5","0", "5", "10", "30", "60"))

# agregar los datos base
pop_by_rp_14[,month:=lubridate::month(pop_by_rp_14$month, label = TRUE)]


#renombrar las columnas
colnames(pop_by_rp_98) <- c("mes", "periodo_retorno", "pop_rp", "pop_frac")
colnames(pop_by_rp_14) <- c("mes", "periodo_retorno", "pop_rp", "pop_frac")

head(pop_by_rp_14)
     mes periodo_retorno  pop_rp    pop_frac
   <ord>          <fctr>   <num>       <num>
1:   Jan               5 5284826 0.345875196
2:   Jan               0 4486888 0.293652670
3:   Jan             -10 1335779 0.087422523
4:   Jan              -5 3626857 0.237366353
5:   Jan              10  397448 0.026011718
6:   Jan              30  147777 0.009671539

Ahora podemos ponerlo todo junto en una visualización.

# Anadir una columna del año. 
pop_by_rp_98$year <- "1998"
pop_by_rp_14$year <- "2014"

# Combinar las tablas
combined_data <- rbind(pop_by_rp_98, pop_by_rp_14)
# crear la paleta para pasar a ggplot
leg_colors<-c(
    # -60 a -30
    '#9B0039',
    #-30 a -10
    '#FF8D43',
    # -10 a -5
    '#FFC754',
    # -5 a 0
    '#FFF4C7',
    # 0 a 5
    '#DDFFFF',
    # 5 a 10
    '#00BFFF',
    # 10 a 30
    '#0000FF',
    # 30 a 60
    '#191970') 


# Crear una combinacion the los resulados
ggplot2::ggplot(combined_data, 
       ggplot2::aes(x = mes, 
           y = pop_frac,
           group = periodo_retorno, 
           fill = periodo_retorno)) +
  ggplot2::geom_bar(stat = "identity", 
           position = "stack", 
           color = "darkgrey") +
  ggplot2::scale_fill_manual(values = leg_colors) +
  ggplot2::ylim(0, 1) +
  ggplot2::labs(title = "Fracción mensual de la población con excedente de agua en Ecuador",
       subtitle = "1998 y 2014, Categorizado por intensidad del período de retorno del excedente",
       x = "",
       y = "Fracción de la población*",
       caption = "*Población derivada de la población mundial cuadriculada (2015)",
       color = "Periodo Retorno", fill = "Periodo Retorno", group = "Periodo Retorno", alpha = "Periodo Retorno") +
  ggplot2::theme_minimal() +
  ggplot2::facet_wrap(~year)+
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust =0.25))

Esta cifra realmente ilustra el impacto humano de la inundacion de 1998 y 2014.

Revisión de conocimientos
  1. Hay varias opciones para subconjuntos (o recortes) espaciales de un objeto ráster a una región de interés. ¿Qué método se utilizó en esta lección?
    1. Usando un vector de fechas.
    2. Usando otro objeto ráster.
    3. Especificación de un cuadro delimitador.
    4. Usando un conjunto de datos de límites vectoriales.
  2. Cuando tenga problemas de memoria, ¿qué puede hacer para reducir la carga computacional?
    1. Trabaje con un período de tiempo o región a la vez.
    2. Guárdelo como un archivo nuevo.
    3. Subconjunto de datos en una región de interés/marco de tiempo.
    4. Encuentre otros datos con los que trabajar.
  3. ¿Cuál es la importancia de subconjuntos de datos?
    1. Liberando espacio.
    2. Analizar un determinado momento o área de interés.
    3. Hacer que el código se ejecute más rápido.
    4. Todo lo anterior.
  4. Conjuntos de datos de población cuadriculados… (seleccione todos los que correspondan)
    1. Muestre dónde están los afluentes de los ríos.
    2. Modele la distribución de la población humana global en celdas de cuadrícula ráster.
    3. Permitir análisis del número de personas impactadas por una amenaza como la sequía.
    4. Variar en la medida en que se utilizan datos auxiliares en su producción.

En esta lección, aprendiste…

¡Felicidades! Ahora deberías poder:

  • Navegue por el sitio web de SEDAC para buscar y descargar conjuntos de datos.
  • Acceda a límites administrativos desde datos de geoBoundaries utilizando API.
  • Subconjunto temporal de un ráster NetCDF utilizando paquetes R como dplyr y lubridate.
  • Recortar una ráster NetCDF con un límite espacial.
  • Escriba un conjunto de datos subconjunto en el disco y cree una imagen para compartir resultados.
  • Identificar áreas de sequía severa y seleccionarlas para su posterior análisis.
  • Resumir datos por condado utilizando la herramienta exactaextractr.
  • Integre el excedente de WSIM-GLDAS, la población de GPW y los datos de límites administrativos de geoBoundaries para crear visualizaciones complejas de series de tiempo.

Lección 2

En esta lección exploramos la sequía de California de 2014. En nuestra próxima lección, examinaremos datos de inundaciones casi en tiempo real en California utilizando el producto de datos MODIS.

Producto de inundaciones globales en tiempo casi real (NRT) de LANCE MODIS

Paso 2, verifique qué días están disponibles para el producto MCDWD L3 F3 accediendo a este enlace: Servidor primario LANCE NRT o Secundario. Nota: debe tener una cuenta de NASA Earthdata.

Según la disponibilidad, edite la variable día_año AAAA-DD. Ejemplo: ‘2022-01’

#agregue el año y la fecha que desea buscar (AAAA-DD, 2022-01)
dia <- '2024-269'

Determinar mosaicos de interés:

Mapa de mosaicos MODIS NRT

Un mapa global con mosaicos etiquetados

Según la disponibilidad, edite la variable Tile_code:

#añadir código de mosaico del mapa de arriba (escrito como h00v00)
codigo_azulejo <- 'h10v09'

Esta es la URL de la API NRT Flood F3 (MCDWD_L3_F3):

# Servidor primario
API_URL <-paste0('https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')

# si el servidor primario no funciona, el servidor secundario puede estar disponible:
#API_URL <-paste0('https://nrt4.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')

Podemos combinar la URL de API anterior con el día_año proporcionado e imprimir los conjuntos de datos disponibles:

#pegando juntos URL y dia
url <- paste0(API_URL, dia)
print(url)
[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2024-269"

Cargando los datos

NASA Earthdata requiere un nombre de usuario y un token de la NASA. Vaya aquí para aprender [Cómo generar un token de usuario] (https://wiki.earthdata.nasa.gov/display/EL/How+to+Generate+a+User+Token).

Una vez que haya generado un token de usuario, podrá acceder a NASA Earthdata utilizando la función GET (wilsjame?). Reemplace USER_TOKEN con el token que generó:

token <- "TOKEN_DEL_USUARIO"

Combine el token con el protador "Bearer" en una sola cadena de texto.

api_key <- paste("Bearer", token)

Este código es para usar en una computadora local, no en 2i2c. Úselo para proteger su token localmente.

# Definir encabezados (incluye tu token)
encabezados <- c (Authorization = api_key)

# Realizar solicitud httr::GET
contenido_respuesta <- httr::GET(url, httr::add_headers(.headers=encabezados))

Verifique el estado de la respuesta desde la función GET:

contenido_respuesta
Response [https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2024-269]
  Date: 2024-09-25 21:13
  Status: 200
  Content-Type: application/json;charset=UTF-8
  Size: 172 kB

A partir de la respuesta del servidor, comprobaremos si la respuesta fue exitosa con if (http_status(contenido_respuesta)$category == "Success"). Si esta afirmación es verdadera, entonces el contenido se asignará a la variable datos en formato JSON, que luego se analiza en un marco de datos usando data_parsed <- jsonlite::fromJSON(data). El marco de datos contiene data_parsed$content, una columna con contenido. Filtramos el contenido por código de mosaico usando el comando content_items <- data_parsed$content[grepl(tile_code, data_parsed$content$name, ignore.case = TRUE),] y agregamos los resultados a un marco de datos.

if (httr::http_status(contenido_respuesta)$category == "Success") {
  # Leer el contenido JSON en un marco de datos
  contenido_respuesta <- httr::content(contenido_respuesta, as = "parsed", simplifyVector = TRUE)
  contenido_respuesta <- contenido_respuesta$content
} else {
  print("La solicitud falló con el código de estado", httr::http_status(contenido_respuesta)$status_code)
}
names(contenido_respuesta)
 [1] "archiveSets"   "cksum"         "collections"   "dataDay"      
 [5] "downloadsLink" "fileId"        "md5sum"        "mtime"        
 [9] "name"          "products"      "resourceType"  "self"         
[13] "size"          "start"         "status"       

Busque el marco de datos y subconjunto de las filas que contienen “codigo_azulejo” en la columna “downloadsLink”.

# Subconjunto de filas donde la columna downloadsLink contiene el código del mosaico
contenido_respuesta <- contenido_respuesta[grepl(codigo_azulejo, contenido_respuesta$downloadsLink),]

contenido_respuesta
    archiveSets      cksum    collections               dataDay
217          61 1681188629 modis-nrt-c6.1 2024-268 = 2024-09-24
                                                                                          downloadsLink
217 https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2024268.h10v09.061.tif
        fileId                           md5sum      mtime
217 2388804642 513b281c967029a292b519ecfede824e 1727231010
                                       name        products resourceType
217 MCDWD_L3_F3_NRT.A2024268.h10v09.061.tif MCDWD_L3_F3_NRT         File
                                                                                                  self
217 https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details/MCDWD_L3_F3_NRT.A2024268.h10v09.061.tif
      size                     start status
217 526436 2024-09-25 02:23:30.99025 Online

Si solo hay 1 fila, seleccione la columna de cadena ‘enlace de descargas’:

contenido_respuesta <- contenido_respuesta$downloadsLink
print(contenido_respuesta)
[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2024268.h10v09.061.tif"

Leer los datos

Utilice la función read_stars() de la biblioteca R stars para leer el ráster geoTiff. El ráster se asigna a la variable raster_df. Esta es una capa derivada de la web que se encuentra en el Sistema de coordenadas geográficas (GCS) WGS84 - Sistema geodésico mundial 1984 EPSG:4326 (número EPSG 4326). Podemos asignar el Sistema de Referencia de Coordenadas (CRS) al objeto estrella:

#read the raster data from the download_link
raster_df <- stars::read_stars(contenido_respuesta)
#Set a CRS for the stars object based on the data documentation
sf::st_crs(raster_df) <- sf::st_crs(4326)

Visualización de datos de inundaciones NRT

Traza el ráster para verlo rápidamente. downsample permite una representación más rápida para la visualización:

plot(raster_df, downsample=30)

Create NRT Flood Plot with Classification

Crear gráfico de inundación NRT con clasificación

Consulte la [Guía del usuario del producto MODIS NRT Global Flood] (https://www.earthdata.nasa.gov/s3fs-public/2023-01/MCDWD_UserGuide_RevC.pdf) para obtener más información.

Los datos de inundaciones de NRT tienen 5 clasificaciones:

Codigo Definicion
0 Sin agua
1 Agua Superficial
2 Inundación recurrente2
3 Inundación (inusual)
255 Datos insuficientes

Para ver los datos de esta clasificación, necesitaremos crear una leyenda clasificada; sin embargo, los datos de NRT Flood se almacenan en números decimales (conocidos como punto flotante). Los saltos de clase nos permiten agrupar datos en función de los valores. Cree saltos de clasificación que divida los datos utilizando los saltos que se enumeran a continuación y sus colores y etiquetas correspondientes. Finalmente, agregue un título para el gráfico que incluya el año y el día, “dia” y “codigo_azulejo”:

#quiebramentos para la categorización del diagrama de tabla
quiebramentos<- c( 0, 1, 2, 3, Inf)
#colores para las pausas argumentales
colores <- c( "gray", "blue", "yellow", "red")
#etiquetas para las quiebras argumentales
etiquetas = c("0: Sin agua", "1: Agua superficial", "2: Inundación recurrente", "3: Inundación (inusual)")
titulo = paste("Inundación casi en tiempo real F3", dia, codigo_azulejo)

Generar un mapa base a partir de Esri Imagery

La función basemap_stars() del paquete stars nos permite acceder a las capas de imágenes de Esri. Para generar un mapa base que muestre la ubicación de nuestro ráster, primero debemos definir el cuadro delimitador con raster_df. Elija “world_imagery” como nuestro fondo y asígnelo al objeto bm_m. La función st_rgb combina la imagen RGB del mapa base en un solo archivo. Finalmente, st_transform transforma el CRS en WGS84.

#definir mapa base usando el ráster como cuadro delimitador
bm_m <- basemaps::basemap_stars(raster_df, map_service = "esri", map_type = "world_imagery")
Loading basemap 'world_imagery' from map service 'esri'...
#La función `st_rgb` nos permite convertir el elemento de estrellas RGB en una sola imagen
bm_m <- stars::st_rgb(bm_m)
#transformar el CRS del mapa base para que coincida con el ráster de inundación
bm_m <- sf::st_transform(bm_m, 4326)

Trazar mapa base y datos de inundaciones NRT

Genere un gráfico desde la biblioteca tmap usando la función tm_shape(). Trazaremos el mapa base y los elementos raster_df en un solo gráfico.

## modo tmap configurado en "trazar"
tmap::tmap_mode("plot")

## el modo tmap también se puede configurar en "ver"
#tmap_mode("view")

#crear un objeto que trace el mapa base y el ráster de inundación NRT
#con la biblioteca tmap, llame a la función tm_shape() para el mapa base
tm_plot <-  tmap::tmap_options(max.raster = c(plot = 50000, view = 50000))+
  tmap::tm_shape(bm_m)+
  #el mapa base bm_m se trata como un ráster
  tmap::tm_raster()+
  #cree una nueva forma de tmap para el ráster de inundación NRT con un estilo como "cat", que significa categórico.
  tmap::tm_shape(raster_df, style="cat")+
  #agregar el estilo de clasificación al ráster, incluidos colores, títulos, saltos de clase y etiquetas
  tmap::tm_raster( palette = c(colores),
  title = titulo, 
  breaks = quiebramentos,
  labels = etiquetas )+
  #estilizar la trama
  tmap::tm_layout(legend.outside = TRUE) +
  tmap::tm_graticules(lines=FALSE)
 
 ##ver la trama
 tm_plot

Para observar una ubicación más de cerca, podemos crear un nuevo cuadro delimitador usando valores de latitud y longitud, recortar los datos y volver a trazarlos en un mapa. Primero, seleccione del mapa cuatro esquinas en grados que incluyan información Norte/Sur y Este/Oeste. Usamos estas coordenadas para crear una matriz de puntos que representan cuatro esquinas; Luego, las cuatro esquinas se utilizan para crear un cuadro delimitador. Las esquinas se colocan primero en el GCS WGS84 para que coincidan con nuestro ráster.

# Definir las coordenadas NWES
norte <- -2 #se utilizan valores negativos para el sur
sur <- -7 #se utilizan valores negativos para el sur
oeste <- -71 #se utilizan números negativos para Oeste
este <- -76 #se utilizan números negativos para Oeste

# Crea un cuadro delimitador
bbox_subset <- sf::st_bbox(c(xmin = oeste, ymin = sur, xmax = este, ymax = norte), crs = 4326)

Rehacemos el mismo proceso anterior para trazar los datos con un mapa base pero con un cuadro delimitador diferente:

#generar un nuevo mapa base
bm_m <- basemaps::basemap_stars(bbox_subset, map_service = "esri", map_type = "world_imagery")
Loading basemap 'world_imagery' from map service 'esri'...
#combinar bandas RGB del mapa base
bm_m <- stars::st_rgb(bm_m)
#transformarse a CRS 4326
bm_m <- sf::st_transform(bm_m, 4326)
## modo tmap configurado en "trazar"
tmap::tmap_mode("plot")

## el modo tmap también se puede configurar en "ver"
#tmap_mode("view")

#crear un objeto que trace el mapa base y el ráster de inundación NRT
#con la biblioteca tmap, llame a la función tm_shape() para el mapa base
tm_plot <- tmap::tmap_options(max.raster = c(plot = 50000, view = 50000))+ 
  tmap::tm_shape(bm_m, raster.downsample=TRUE)+ 
  tmap::tm_raster()+ 
  #cree una nueva forma de tmap para el ráster de inundación NRT con el estilo "gato", que significa categórico.
  tmap::tm_shape(raster_df, style="cat", raster.downsample=TRUE)+ 
  #agregar el estilo de clasificación al ráster 
  tmap::tm_raster( palette = c(colores), title = titulo, breaks = quiebramentos, labels = etiquetas )+ 
  #estilizar la trama
  tmap::tm_layout(legend.outside = TRUE) + 
  tmap::tm_graticules(lines=FALSE)

#Ver la trama:
tm_plot

¡Felicitaciones! Pudimos descargar y explorar los datos de inundaciones del NRT de la NASA para elegir un área de interés en Ecuador y visualizar los resultados. Ahora debería poder:

  • Navegar por el sitio web de datos de LANCE y determinar qué datos están disponibles.

  • Seleccionar un mosaico y una fecha para descargar los datos del NRT.

  • Crear una solicitud HTTP GET para descargar datos casi en tiempo real.

  • Trazar en un mapa y clasificar los datos ráster para determinar las áreas con inundaciones inusuales.

References

Bland, Alastair. 2014. “California Drought Has Wild Salmon Competing with Almonds for Water.” National Public Radio. https://www.npr.org/sections/thesalt/2014/08/21/342167846/california-drought-has-wild-salmon-competing-with-almonds-for-water.
Center For International Earth Science Information Network-CIESIN-Columbia University. 2018. “Gridded Population of the World, Version 4 (GPWv4): Population Count, Revision 11.” Palisades, NY: NASA Socioeconomic Data; Applications Center (SEDAC). https://doi.org/10.7927/H4JW8BX5.
Daniel Baston. 2023. Exactextractr: Fast Extraction from Raster Datasets Using Polygons. https://isciences.gitlab.io/exactextractr/.
“Evidence Suggests California’s Drought Is the Worst in 1,200 Years.” 2014. Woods Hole Oceanographic Institution. https://phys.org/news/2014-12-evidence-california-drought-worst-years.html.
Guillen, George. 2002. “Klamath River Fish Die-Off.” Mortality Report AFWO-01-03. Arcata, CA: U.S. Fish & Wildlife Service. https://www.trrp.net/DataPort/doc.php?id=301.
IPCC. 2023. “Climate Change 2021 the Physical Science Basis,” June. https://doi.org/10.1017/9781009157896.
ISciences, and Center For International Earth Science Information Network-CIESIN-Columbia University. 2022a. “Documentation for the Water Security Indicator Model - Global Land Data Assimilation System (WSIM-GLDAS) Monthly Grids, Version 1.” Palisades, NY: NASA Socioeconomic Data; Applications Center (SEDAC). https://doi.org/10.7927/X7FJ-JJ41.
———. 2022b. “Water Security Indicator Model - Global Land Data Assimilation System (WSIM-GLDAS) Monthly Grids, Version 1.” Palisades, NY: NASA Socioeconomic Data; Applications Center (SEDAC). https://doi.org/10.7927/Z1FN-KF73.
Leyk, Stefan, Andrea E. Gaughan, Susana B. Adamo, Alex de Sherbinin, Deborah Balk, Sergio Freire, Amy Rose, et al. 2019. “The Spatial Allocation of Population: A Review of Large-Scale Gridded Population Data Products and Their Fitness for Use.” Earth System Science Data 11 (3): 1385–1409. https://doi.org/10.5194/essd-11-1385-2019.
NOAA. 2019. “Water Cycle.” National Oceanic; Atmospheric Administration. https://www.noaa.gov/education/resource-collections/freshwater/water-cycle.
Rodgers, Alison Ince. 2023. “Understanding Droughts.” National Geographic Society. https://education.nationalgeographic.org/resource/understanding-droughts/.
Sanders, Sam. 2014. “Despite California’s Drought, Taps Still Flowing in LA County.” National Public Radio. https://www.npr.org/2014/07/20/333019977/despite-californias-drought-taps-still-flowing-in-la-county.
Zhou, Haddeland, T. 2016. “Human-Induced Changes in the Global Water Cycle.” Geophysical Monograph Series. https://doi.org/10.1002/9781118971772.ch4.

Footnotes

  1. Crédito de la foto, NASA/JPL.↩︎

  2. El valor 2 (inundación recurrente) no se completa en la versión beta.↩︎