<- c("stars", "terra", "sf", "cubelyr", "lubridate", "httr", "data.table", "exactextractr", "ggplot2", "kableExtra", "jsonlite", "tmap", "basemaps", "sp")
packages_to_check
# 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)
}
WSIM-GLDAS: Adquisición, Exploración e Integración
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).
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.
- ¿Cómo describirías mejor el ciclo del agua?
- Un período prolongado de poca o ninguna lluvia.
- Baja precipitación combinada con temperaturas récord.
- La circulación del agua en y sobre la superficie de la Tierra.
- Un ciclo que ocurre debido a la sequía.
- ¿Qué intervenciones humanas afectan el ciclo del agua (selecciona todas las que apliquen)?
- Emisiones de gases de efecto invernadero.
- Cambios en el uso del suelo.
- Desarrollo de presas y embalses.
- Sobreexplotación de aguas subterráneas.
- ¿Qué es un déficit de precipitación?
- Un período de precipitación por debajo del promedio.
- Un período prolongado de poca o ninguna lluvia.
- Un período de reacciones en cadena.
- Un período de precipitación por encima del promedio.
Adquisición de los 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
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.
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
<- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_12mo.nc", proxy = TRUE) wsim_gldas
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')
<- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_12mo.nc", proxy = FALSE, sub = 'surplus') wsim_gldas
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
<-seq(lubridate::ymd("2000-12-01"),
keeps::ymd("2014-12-01"),
lubridateby = "year")
#cambiar tipo de datos a POSIXct
<- as.POSIXct(keeps)
keeps # filtrar el conjunto de datos usando el vector de fechas
<- dplyr::filter(wsim_gldas, time %in% keeps)
wsim_gldas_2000_2014 # 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
::slice(index = 1:6, along = "time") |>
dplyr# 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.
- ¿Cuál de estos describe mejor un conjunto de datos ráster?
Un tipo de datos geográficos formados por celdas cuadradas.
Una tabla o lista de números geográficos.
Una región geográfica de interés.
Un tipo de datos geográficos formados por líneas y puntos
- ¿Cuál de las siguientes afirmaciones es cierta sobre la información que pueden almacenar los rásteres? (seleccione todas las que correspondan)
- Atributos (contenido temático)
- Dimensiones (información que expresa información de extensión espacial o temporal)
- Coordenadas geográficas
- Una lista de números
Selección espacial
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/
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.
<- httr::GET("https://www.geoboundaries.org/api/current/gbOpen/ECU/ADM2/") ecuador_ad2
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
<- httr::content(ecuador_ad2)
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
<- sf::st_read(ecuador_ad2$gjDownloadURL) ecuador_ad2
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
$shapeName ecuador_ad2
[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
<-
estadosc("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[(ecuador_ad2$shapeName %in% estados),]
ecuador_ad2_sub# 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
<- ecuador_ad2[ecuador_ad2$shapeName == "Zapotillo",]
ecu_adm1 # 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[ecuador_ad2]
wsim_gldas 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 ::slice(index = 15, along = "time") |>
dplyr# 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.
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
<-seq(lubridate::ymd("1994-12-01"),
keeps::ymd("2001-12-01"),
lubridateby = "year")
#cambiar tipo de datos a POSIXct
<- as.POSIXct(keeps)
keeps # filtrar usando ese vector
<- dplyr::filter(wsim_gldas, time %in% keeps) wsim_gldas
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.
- ¿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)
- Subconjunto espacial de los datos al Ecuador
- Leyendo en sólo un año de datos
- Leer sólo el atributo de interés (es decir, ‘surplus’)
- 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
::ggplot()+
ggplot2# visualizar el objeto
::geom_stars(data = wsim_gldas)+
stars# establecer coordenadas iguales para los ejes
::coord_equal()+
ggplot2# crear múltiples paneles para cada paso de tiempo
::facet_wrap(~time)+
ggplot2# trazar el límite de Ecuador solo con el contorno
# ggplot2::geom_sf(fill = NA, lwd=0)+
# agregar la paleta para wsim-gldas
::scale_fill_stepsn(
ggplot2colors = 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
::labs(
ggplot2title="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
::theme_minimal()+
ggplot2# desactivar algunos elementos gráficos adicionales
::theme(
ggplot2axis.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!
- Los mapas de salida de Anomalías del Execdente Medio Anual para Ecuador indican que…
- El exceso más significativo en 1998 se localiza en el oeste de Ecuador.
- El exceso menos significativo de 2000 se sitúa en el Sur.
- El exceso más significativo en 1996 se localiza en los estados vecinos.
- 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 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
<- stars::read_stars("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_1mo.nc", proxy = FALSE, sub = 'surplus') wsim_gldas_1mo_98
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
<-seq(lubridate::ymd("1998-01-01"),
keeps::ymd("1998-12-01"),
lubridateby = "month")
# cambiar tipo de datos a POSIXct
<- as.POSIXct(keeps)
keeps # filtrar usando ese vector
<- dplyr::filter(wsim_gldas_1mo_98, time %in% keeps)
wsim_gldas_1mo_98 # 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[ecuador_ad2]
wsim_gldas_1mo_98 # dar etiquetas bonitas a la dimensión del tiempo
<-
wsim_gldas_1mo_98 |>
wsim_gldas_1mo_98 ::st_set_dimensions("time", values = month.name)
stars
# parcelas mensuales de Ecuador
# cargar los datos base de Ecuador
::ggplot()+
ggplot2# cargar el objeto wsim stars
::geom_stars(data = wsim_gldas_1mo_98)+
stars# relación de aspecto igual para el hacha lat/long
::coord_equal()+
ggplot2# crear una trama secundaria para cada paso de tiempo
::facet_wrap(~time)+
ggplot2# configurar la paleta wsim
::scale_fill_stepsn(
ggplot2colors = 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
::labs(
ggplot2title="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
::theme_minimal()+
ggplot2# eliminar elementos adicionales para hacer un mapa más limpio
::theme(
ggplot2axis.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
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
<-stats::na.omit(excedente_hist)
excedente_hist
# crear el histograma
::ggplot(excedente_hist, ggplot2::aes(surplus))+
ggplot2# darle estilo a las barras
::geom_histogram(binwidth = 6, fill = "#325d88")+
ggplot2# subtrama para cada paso de tiempo
::facet_wrap(~time)+
ggplot2# usa el tema mínimo
::theme_minimal() ggplot2
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.
- ¿Cuál es otro término para “histograma”?
- coropleta
- Distribución de frecuencias
- Matriz de datos
- 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.
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
<-terra::sds("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/composite_1mo.nc") wsim_gldas_1mo
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["surplus"]
wsim_gldas_1mo# crear una secuencia de fechas para mantener
<-seq(lubridate::ymd("2014-01-01"), lubridate::ymd("2014-12-01"), by = "month")
keeps# subconjunto del objeto terra solo para esos pasos de tiempo del 2014
<-wsim_gldas_1mo[[terra::time(wsim_gldas_1mo) %in% keeps]]
wsim_gldas_1mo_14# 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
<-seq(lubridate::ymd("1998-01-01"), lubridate::ymd("1998-12-01"), by = "month")
keeps# subconjunto del objeto terra solo para esos pasos de tiempo del 1998
<-wsim_gldas_1mo[[terra::time(wsim_gldas_1mo) %in% keeps]]
wsim_gldas_1mo_98# 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::exact_extract(
exactextractr# 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::exact_extract(
exactextractr# 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
<-cbind(ecuador_ad2, ecu_resumen_98)
ecuador_ad2<-cbind(ecuador_ad2, ecu_resumen_14) ecuador_ad2
# Combinar las dos tables 98 y 14
<- cbind(ecu_resumen_98, ecu_resumen_14)
combined_data # preparar una versión para crear una tabla
<-cbind(Area=ecuador_ad2$shapeName,
ecu_tablaround(combined_data))
# crear la tabla con solo las primeras filas y 25 columnas
::kbl(ecu_tabla[,c(1:25)]) |>
kableExtra::kable_styling(
kableExtrabootstrap_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).
- ¿Para qué se pueden utilizar las estadísticas zonales?
- Subcomponer datos a un período de tiempo de interés.
- Creando una serie de tiempo.
- Resumir los valores de las celdas que se encuentran dentro de un límite.
- 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.
- Los mapas coropléticos, también conocidos como mapas temáticos, son útiles para
- Visualizar datos en geografías familiares como condados o estados.
- Encontrar direcciones de un lugar a otro
- Visualización de datos en unidades geográficas uniformes, como celdas de cuadrícula ráster.
- 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.
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::rast("/vsicurl/https://s3-us-west-2.amazonaws.com/tops-school/water-module/gpw_v4_population_density_rev11_2015_15_min.tif")
terra
# 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
<- c(-60, -30, -30,
m -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
<- matrix(m, ncol=3, byrow=TRUE)
rclmat # clasificar los datos
<-
wsim_gldas_1mo_class_98 ::classify(wsim_gldas_1mo_98, rclmat, include.lowest = TRUE)
terra<-
wsim_gldas_1mo_class_14 ::classify(wsim_gldas_1mo_14, rclmat, include.lowest = TRUE) terra
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:
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.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 ::exact_extract(wsim_gldas_1mo_class_98, ecuador_ad2, function(df) {
exactextractr<- data.table::setDT(df)},
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 ::exact_extract(wsim_gldas_1mo_class_14, ecuador_ad2, function(df) {
exactextractr<- data.table::setDT(df)},
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 argumentoinclude_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
a2014-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ásterwsim_gldas_1mo
.weights
: la columnaweight
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.
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 ::melt(
data.table# 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 ::melt(
data.table# 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_rp = round(sum(coverage_area * weight) / 1e6)), by = .(month, return_period)]
pop_by_rp_98[, .(# algunas celdas no tienen valores de período de retorno deficitario y dan como resultado NaN; simplemente elimínelas
<- na.omit(pop_by_rp_98)
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)
:= sum(pop_rp), by = month]
pop_by_rp_98[, total_pop # 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_rp = round(sum(coverage_area * weight) / 1e6)), by = .(month, return_period)]
pop_by_rp_14[, .(# algunas celdas no tienen valores de período de retorno deficitario y dan como resultado NaN; simplemente elimínelas
<- na.omit(pop_by_rp_14)
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)
:= sum(pop_rp), by = month]
pop_by_rp_14[, total_pop # 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_rp / total_pop][, total_pop := NULL]
pop_by_rp_98[, pop_frac # calcular la fracción de esa combinación de clase mes-wsim en relación con la población total
:= pop_rp / total_pop][, total_pop := NULL]
pop_by_rp_14[, pop_frac
# 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.
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
$return_period<-
pop_by_rp_98factor(pop_by_rp_98$return_period,
levels = c("-60", "-30", "-10", "-5","0", "5", "10", "30", "60"))
# agregar los datos base
:=lubridate::month(pop_by_rp_98$month, label = TRUE)]
pop_by_rp_98[,month
# ggplot es más fácil con factores
$return_period<-
pop_by_rp_14factor(pop_by_rp_14$return_period,
levels = c("-60", "-30", "-10", "-5","0", "5", "10", "30", "60"))
# agregar los datos base
:=lubridate::month(pop_by_rp_14$month, label = TRUE)]
pop_by_rp_14[,month
#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.
$year <- "1998"
pop_by_rp_98$year <- "2014"
pop_by_rp_14
# Combinar las tablas
<- rbind(pop_by_rp_98, pop_by_rp_14) combined_data
# crear la paleta para pasar a ggplot
<-c(
leg_colors# -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
::ggplot(combined_data,
ggplot2::aes(x = mes,
ggplot2y = pop_frac,
group = periodo_retorno,
fill = periodo_retorno)) +
::geom_bar(stat = "identity",
ggplot2position = "stack",
color = "darkgrey") +
::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",
ggplot2subtitle = "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") +
::theme_minimal() +
ggplot2::facet_wrap(~year)+
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust =0.25)) ggplot2
Esta cifra realmente ilustra el impacto humano de la inundacion de 1998 y 2014.
- 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?
- Usando un vector de fechas.
- Usando otro objeto ráster.
- Especificación de un cuadro delimitador.
- Usando un conjunto de datos de límites vectoriales.
- Cuando tenga problemas de memoria, ¿qué puede hacer para reducir la carga computacional?
- Trabaje con un período de tiempo o región a la vez.
- Guárdelo como un archivo nuevo.
- Subconjunto de datos en una región de interés/marco de tiempo.
- Encuentre otros datos con los que trabajar.
- ¿Cuál es la importancia de subconjuntos de datos?
- Liberando espacio.
- Analizar un determinado momento o área de interés.
- Hacer que el código se ejecute más rápido.
- Todo lo anterior.
- Conjuntos de datos de población cuadriculados… (seleccione todos los que correspondan)
- Muestre dónde están los afluentes de los ríos.
- Modele la distribución de la población humana global en celdas de cuadrícula ráster.
- Permitir análisis del número de personas impactadas por una amenaza como la sequía.
- 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)
<- '2024-269' dia
Determinar mosaicos de interés:
Según la disponibilidad, edite la variable Tile_code:
#añadir código de mosaico del mapa de arriba (escrito como h00v00)
<- 'h10v09' codigo_azulejo
Esta es la URL de la API NRT Flood F3 (MCDWD_L3_F3):
# Servidor primario
<-paste0('https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')
API_URL
# 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
<- paste0(API_URL, dia)
url 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_DEL_USUARIO" token
Combine el token con el protador "Bearer"
en una sola cadena de texto.
<- paste("Bearer", token) api_key
Este código es para usar en una computadora local, no en 2i2c. Úselo para proteger su token localmente.
# Definir encabezados (incluye tu token)
<- c (Authorization = api_key)
encabezados
# Realizar solicitud httr::GET
<- httr::GET(url, httr::add_headers(.headers=encabezados)) contenido_respuesta
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
<- httr::content(contenido_respuesta, as = "parsed", simplifyVector = TRUE)
contenido_respuesta <- contenido_respuesta$content
contenido_respuesta 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[grepl(codigo_azulejo, contenido_respuesta$downloadsLink),]
contenido_respuesta
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$downloadsLink
contenido_respuesta 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
<- stars::read_stars(contenido_respuesta)
raster_df #Set a CRS for the stars object based on the data documentation
::st_crs(raster_df) <- sf::st_crs(4326) sf
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
<- c( 0, 1, 2, 3, Inf)
quiebramentos#colores para las pausas argumentales
<- c( "gray", "blue", "yellow", "red")
colores #etiquetas para las quiebras argumentales
= c("0: Sin agua", "1: Agua superficial", "2: Inundación recurrente", "3: Inundación (inusual)")
etiquetas = paste("Inundación casi en tiempo real F3", dia, codigo_azulejo) titulo
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
<- basemaps::basemap_stars(raster_df, map_service = "esri", map_type = "world_imagery") bm_m
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
<- stars::st_rgb(bm_m)
bm_m #transformar el CRS del mapa base para que coincida con el ráster de inundación
<- sf::st_transform(bm_m, 4326) bm_m
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_mode("plot")
tmap
## 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
<- tmap::tmap_options(max.raster = c(plot = 50000, view = 50000))+
tm_plot ::tm_shape(bm_m)+
tmap#el mapa base bm_m se trata como un ráster
::tm_raster()+
tmap#cree una nueva forma de tmap para el ráster de inundación NRT con un estilo como "cat", que significa categórico.
::tm_shape(raster_df, style="cat")+
tmap#agregar el estilo de clasificación al ráster, incluidos colores, títulos, saltos de clase y etiquetas
::tm_raster( palette = c(colores),
tmaptitle = titulo,
breaks = quiebramentos,
labels = etiquetas )+
#estilizar la trama
::tm_layout(legend.outside = TRUE) +
tmap::tm_graticules(lines=FALSE)
tmap
##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
<- -2 #se utilizan valores negativos para el sur
norte <- -7 #se utilizan valores negativos para el sur
sur <- -71 #se utilizan números negativos para Oeste
oeste <- -76 #se utilizan números negativos para Oeste
este
# Crea un cuadro delimitador
<- sf::st_bbox(c(xmin = oeste, ymin = sur, xmax = este, ymax = norte), crs = 4326) bbox_subset
Rehacemos el mismo proceso anterior para trazar los datos con un mapa base pero con un cuadro delimitador diferente:
#generar un nuevo mapa base
<- basemaps::basemap_stars(bbox_subset, map_service = "esri", map_type = "world_imagery") bm_m
Loading basemap 'world_imagery' from map service 'esri'...
#combinar bandas RGB del mapa base
<- stars::st_rgb(bm_m)
bm_m #transformarse a CRS 4326
<- sf::st_transform(bm_m, 4326) bm_m
## modo tmap configurado en "trazar"
::tmap_mode("plot")
tmap
## 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
<- tmap::tmap_options(max.raster = c(plot = 50000, view = 50000))+
tm_plot ::tm_shape(bm_m, raster.downsample=TRUE)+
tmap::tm_raster()+
tmap#cree una nueva forma de tmap para el ráster de inundación NRT con el estilo "gato", que significa categórico.
::tm_shape(raster_df, style="cat", raster.downsample=TRUE)+
tmap#agregar el estilo de clasificación al ráster
::tm_raster( palette = c(colores), title = titulo, breaks = quiebramentos, labels = etiquetas )+
tmap#estilizar la trama
::tm_layout(legend.outside = TRUE) +
tmap::tm_graticules(lines=FALSE)
tmap
#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.