Taller 1: Preparación de datos Región de Coquimbo

Preparación de datos para predicción de temperatura en la Región de Coquimbo

Objetivo general

El objetivo de este taller es probar varios metodos de predicción espacial para espacializar la temperatura en la Región de Coquimbo, CHILE, para el año 2023 con enfoque en los meses de enero, junio, octubre y diciembre. En esta sección se desea preparar los datos para realizar la predicción de temperatura basada en diversos predictores, principalmente de la temperatura in situ obtenida desde estaciones meteorologicas, como variables predictorias obtenidas a partir de datos raster tales como la elevación, indices topográficos (pendiente y orientación de laderas), NDVI, distancia de la costa y temperatura superficial mediante sensores satelitales (MODIS).

Objetivo Especifico

1.- Preparación de datos tabulados y raster de los meses indicados y por las estaciones ubicadas en la Región de Coquimbo. 2.- Realizar un analisis exploratorio de datos básico.

1. Librerias

Se llaman las librerias a utilizar para este proyecto

library(remotes) #paquete que trae función para instala un paquete desde github
library(tidyverse) #colección de paquetes para análisis de datos

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf) #manejo de datos vectoriales

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(terra) #manejo de datos raster y vectoriales

terra 1.7.78

Adjuntando el paquete: 'terra'

The following object is masked from 'package:tidyr':

    extract

library(tmap) # creación de mapas temáticos

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')

library(geodata) # acceso a descarga de diferentes datos raster y vectoriales
library(rstac) # acceder a spatial temporal assets catalogs (STAC)
library(agrometR) #acceder a datos de estaciones red Agromet
library(earthdatalogin) #para acceder a productos NASA 'EarthData' 
library(psych) # para estadistica

Adjuntando el paquete: 'psych'

The following objects are masked from 'package:terra':

    describe, distance, rescale

The following objects are masked from 'package:ggplot2':

    %+%, alpha

library(gt) #paquete para la creación de tablas
library(dplyr) #para manejo de dataframe

Es necesario establecer un directorio de trabajo

setwd("~/00_tallegeo/taller01_geost")

2. Descargar datos vectoriales y de elevación de los límites de Coquimbo

2.1 Datos vectoriales

regiones <- gadm('chile', path = 'data/raw')# descarga limites administrativos Chile

regiones_sf <- regiones |> 
  st_as_sf()
class(regiones_sf) #se identifica el tipo de dato

[1] "sf"         "data.frame"

#se plotean las regiones del pais
plot(sf::st_geometry(regiones_sf)) 

El signo $ es para traer indexación y un elemento especifico

regiones_sf$NAME_1

 [1] "Antofagasta"                      "Araucanía"                       
 [3] "Arica y Parinacota"               "Atacama"                         
 [5] "Aysén del General Ibañez del Cam" "Bío-Bío"                         
 [7] "Coquimbo"                         "Libertador General Bernardo O'Hi"
 [9] "Los Lagos"                        "Los Ríos"                        
[11] "Magallanes y Antártica Chilena"   "Maule"                           
[13] "Ñuble"                            "Santiago Metropolitan"           
[15] "Tarapacá"                         "Valparaíso"                      

Selecciona la geometria (fila 7) corresponde a la regiona de Coquimbo

reg_coq <- regiones_sf[7, ] 
plot(sf::st_geometry(reg_coq), main = "Región de Coquimbo")

Se exportan los datos vectoriales a utilizar

write_sf(reg_coq,'data/raw/region_coquimbo.gpkg')
write_rds(reg_coq,'data/raw/region_coquimbo.rds')

2.2. Descargar datos de elevación

Se preparan los datos raster

coords <- reg_coq |> 
  st_centroid() |> #funcion de centroide devuelve un objeto sf
  st_coordinates() #convierte el objeto sf en matrix

Warning: st_centroid assumes attributes are constant over geometries

#descargar datos de elevacion a 30 m/pix 
elev <- geodata::elevation_30s('chile', path='data/raw')

#reproyectar (raster) a 32719 utm zona 19s datum wgs84 
elev_proj <- project(elev,'EPSG:32719')

#reproyecta objeto sf (vectorial) a utm zona 19s datum wgs84
reg_coq_utm <- reg_coq |> 
  st_transform(32719)

Cortar por área de interes

elev_proj_coq <- crop(elev_proj,reg_coq_utm)
plot(elev_proj_coq)
plot(reg_coq_utm,add= TRUE)

Warning in plot.sf(reg_coq_utm, add = TRUE): ignoring all but the first
attribute

Aplicacion de máscara

elev_proj_coq <- mask(elev_proj_coq, reg_coq_utm)

#Plotear finalmente el dem con la mascara aplicada
plot(elev_proj_coq, main="Modelo Digital de Elevación")

2.3. Obtener datos de Slope y Aspect

Generación de datos

# obtencion de mapa de laderas
elev_asp <- terrain(elev_proj_coq,'aspect')
# obtencion capa raster de pendientes
elev_slope <- terrain(elev_proj_coq,'slope')

Ploteo rápido de mapa de Pendientes (Slppe)

plot(elev_slope, main="Modelo de Pendientes")

Ploteo rápido de mapa de Laderas (Aspect)

plot(elev_asp, main="Modelo de Orientación de Laderas")

Se pueden unir las tres capas y visualizar

elev_pred <- c(elev_proj_coq,elev_asp,elev_slope)
names(elev_pred) <- c('dem','aspect','slope')
plot(elev_pred)

Los datos son exportados

writeRaster(elev_pred,'data/predictorias/elev_predict_Coq.tif',overwrite=TRUE)

3. Preparacion de datos para determinar la distancia a la costa

Preparacion de limites de chile

lim_ch <- gadm('chile',level=0, path= 'data/raw')

lim_chl <- lim_ch |> 
  st_as_sf() |> 
  st_geometry() |> 
  st_transform(32719)

#Ver mapa general de CHile y la ubicación de la Región
plot(st_geometry(lim_chl))
plot(st_geometry(reg_coq_utm),col='brown', add= TRUE)

Obtención de datos para el cálculo de extracción de costa

lim_chl_linea <- lim_chl |> 
  st_simplify(dTolerance = 1000) |> 
  st_cast('MULTILINESTRING')
reg_coq_linea <- reg_coq_utm |> 
  st_simplify(dTolerance = 1000) |> 
  st_cast('MULTILINESTRING')

costa <- st_intersection(reg_coq_linea,lim_chl_linea)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

m <- matrix(c(-70.511, -72.581, -74.759, -72.549,-70.511, -28.834, -39.706, -39.294, -28.575,-28.834),ncol=2)

pol <- st_polygon(list(m)) |> 
  st_sfc(crs=4326) |> 
  st_transform(32719)
costa <- st_intersection(costa,pol)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

4. Cálculo de línea de Costa

Creacion raster de 500 metros de resolución

bb_num <- reg_coq_utm |> 
  sf::st_bbox()
r<- rast(xmin=bb_num[1],
         xmax=bb_num[3],
         ymin=bb_num[2],
         ymax=bb_num[4],
         res=500, crs='EPSG:32719')
xy <-terra::xyFromCell(r,1:ncell(r))
r_cent <- xy |> 
  as_tibble() |> 
  st_as_sf(coords=c('x','y'),crs=32719)

Calculo de linea de costa y ploeto de mapa

dist_cost <- st_distance(costa,r_cent)
values(r) <- as.numeric(dist_cost)
dcosta <- mask(r,vect(reg_coq_utm))
names(dcosta) <- 'dist_costa'
plot(dcosta, main="Distancia de Costa en metros")

Se exporta el resultado

write_rds(dcosta,'data/predictorias/distancia_costa_Coq.rds')

5. Descarga de datos STAC-NASA (cmr.earthdata.nasa.gov)

Logeo con pagina de descarga de datos satelital y preparacion de descarga de datos. Se deja programado los meses de interés y los directorios para almacenar las imágenes satelitales.

earthdatalogin::edl_netrc() #login al servicio cmr.earthdata.nasa.gov

#productos a descargar almacenados como carácter en un vector
products <- c('MOD13A3.v061','MOD11A2.v061')

# meses de interes para extraer observaciones
mes <- c('enero'="2023-01-01T00:00:00Z/2023-01-31T00:00:00Z",
         'junio'="2023-06-01T00:00:00Z/2023-06-30T00:00:00Z",
         'octubre'="2023-10-01T00:00:00Z/2023-10-31T00:00:00Z",
         'diciembre'="2023-12-01T00:00:00Z/2023-12-31T00:00:00Z")# meses de interes para extraer observaciones
         
region <- read_rds('data/raw/region_coquimbo.rds')

bb <- region |> st_transform(4326) |> 
  st_bbox() |> 
  as.numeric()

o <- lapply(products, \(prod) {
  lapply(1:4, \(i) {
    m <- mes[i]
    print(paste(m, prod))
    
    items <- stac("https://cmr.earthdata.nasa.gov/stac/LPCLOUD") |> 
      stac_search(collections = prod,
                  bbox = bb,
                  datetime = m) |> 
      post_request() |> 
      items_fetch()
    
    items$features |> 
      purrr::map(list("assets", 6, 'href')) |> 
      purrr::map_chr(\(x) edl_download(x, dest = file.path('data/raw/modis/', prod, names(m), basename(x))))
  })
})

[1] "2023-01-01T00:00:00Z/2023-01-31T00:00:00Z MOD13A3.v061"
[1] "2023-06-01T00:00:00Z/2023-06-30T00:00:00Z MOD13A3.v061"
[1] "2023-10-01T00:00:00Z/2023-10-31T00:00:00Z MOD13A3.v061"
[1] "2023-12-01T00:00:00Z/2023-12-31T00:00:00Z MOD13A3.v061"
[1] "2023-01-01T00:00:00Z/2023-01-31T00:00:00Z MOD11A2.v061"
[1] "2023-06-01T00:00:00Z/2023-06-30T00:00:00Z MOD11A2.v061"
[1] "2023-10-01T00:00:00Z/2023-10-31T00:00:00Z MOD11A2.v061"
[1] "2023-12-01T00:00:00Z/2023-12-31T00:00:00Z MOD11A2.v061"

5.1. Obtención de NDVI

Preparación de rutas

# Rutas
files_ene <- list.files('data/raw/modis/MOD13A3.v061/enero',full.names = TRUE)
files_jun <- list.files('data/raw/modis/MOD13A3.v061/junio',full.names = TRUE)
files_oct <- list.files('data/raw/modis/MOD13A3.v061/octubre',full.names = TRUE)
files_dic <- list.files('data/raw/modis/MOD13A3.v061/diciembre',full.names = TRUE)

NDVI Enero

# NDVI Enero
ndvi_ene <- files_ene |> 
  purrr::map(rast,lyrs=1) |> #itera y carga cada archivo y los almacena en una lista
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |>  #crea el mosaico
  terra::project('EPSG:32719') |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |> #aplica mascara para el área de estudio
  terra::trim() |> #elimina los NA sobrantes alrededor de la extensión
  terra::app(\(x) x*1e-8) # transforma los valorea al rango 0-1

ndvi_ene <- terra::resample(ndvi_ene,r) #resamplea para igual resolución
names(ndvi_ene) <- 'ndvi_ene' #cambia nombre
plot(ndvi_ene, main="NDVI enero")

NDVI Junio

ndvi_jun <- files_jun |> 
  purrr::map(rast,lyrs=1) |> #itera y carga cada archivo y los almacena en una lista
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |>  #crea el mosaico
  terra::project('EPSG:32719') |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |> #aplica mascara para el área de estudio
  terra::trim() |> #elimina los NA sobrantes alrededor de la extensión
  terra::app(\(x) x*1e-8) # transforma los valorea al rango 0-1

ndvi_jun <- terra::resample(ndvi_jun,r) #resamplea para igual resolución
names(ndvi_jun) <- 'ndvi_jun' #cambia nombre
plot(ndvi_jun, main="NDVI junio")

NDVI Octubre

ndvi_oct <- files_oct |> 
  purrr::map(rast,lyrs=1) |> #itera y carga cada archivo y los almacena en una lista
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |>  #crea el mosaico
  terra::project('EPSG:32719') |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |> #aplica mascara para el área de estudio
  terra::trim() |> #elimina los NA sobrantes alrededor de la extensión
  terra::app(\(x) x*1e-8) # transforma los valorea al rango 0-1

ndvi_oct <- terra::resample(ndvi_oct,r) #resamplea para igual resolución
names(ndvi_oct) <- 'ndvi_oct' #cambia nombre
plot(ndvi_oct, main="NDVI octubre")

NDVI Diciembre

ndvi_dic <- files_dic |> 
  purrr::map(rast,lyrs=1) |> #itera y carga cada archivo y los almacena en una lista
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |>  #crea el mosaico
  terra::project('EPSG:32719') |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |> #aplica mascara para el área de estudio
  terra::trim() |> #elimina los NA sobrantes alrededor de la extensión
  terra::app(\(x) x*1e-8) # transforma los valorea al rango 0-1

ndvi_dic <- terra::resample(ndvi_dic,r) #resamplea para igual resolución
names(ndvi_dic) <- 'ndvi_dic' #cambia nombre
plot(ndvi_dic, main="NDVI diciembre")

Ploteo general de los cuatro meses

ndvi_coq <- c(ndvi_ene,ndvi_jun,ndvi_oct,ndvi_dic)
plot(ndvi_coq)

5.2. Obtención de LST (Land Surface Temperature)

Tile de la área de estudio son:

tiles <- c('h11v11','h11v12','h12v12') #tiles que corresponden al área de estudio

LST Enero

lst_out_ene <- purrr::map(tiles,\(tile){
  files_lst_ene <- list.files('data/raw/modis/MOD11A2.v061/enero',full.names = TRUE,pattern = tile)
  
  lst_regs_ene <- files_lst_ene |> 
    purrr::map(rast,lyrs = 'LST_Day_1km') |> #extrae la banda por el nombre, también se podría haber usado el índice (1)
    terra::rast() |> 
    terra::app(mean)
})

lst_ene <- lst_out_ene |> 
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |> # aplica el mosaico para unir los tres tiles
  terra::project("EPSG:32719") |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |>  #aplica mascara para el área de estudio
  terra::trim() |>  #elimina NA sobrantes en los bordes
  terra::app(\(x) x-273.15) #convierte la T° de Kelvin a °C

# resampleo

lst_ene <- resample(lst_ene,r)
names(lst_ene) <- 'lst_ene'
plot(lst_ene, main="LST enero")

LST Junio

lst_out_jun <- purrr::map(tiles,\(tile){
  files_lst_jun <- list.files('data/raw/modis/MOD11A2.v061/junio',full.names = TRUE,pattern = tile)
  
  lst_regs_jun <- files_lst_jun |> 
    purrr::map(rast,lyrs = 'LST_Day_1km') |> #extrae la banda por el nombre, también se podría haber usado el índice (1)
    terra::rast() |> 
    terra::app(mean)
})

lst_jun <- lst_out_jun |> 
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |> # aplica el mosaico para unir los tres tiles
  terra::project("EPSG:32719") |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |>  #aplica mascara para el área de estudio
  terra::trim() |>  #elimina NA sobrantes en los bordes
  terra::app(\(x) x-273.15) #convierte la T° de Kelvin a °C

# resampleo

lst_jun <- resample(lst_jun,r)
names(lst_jun) <- 'lst_jun'
plot(lst_jun, main="LST junio")

LST Octubre

lst_out_oct <- purrr::map(tiles,\(tile){
  files_lst_oct <- list.files('data/raw/modis/MOD11A2.v061/octubre',full.names = TRUE,pattern = tile)
  
  lst_regs_oct <- files_lst_oct |> 
    purrr::map(rast,lyrs = 'LST_Day_1km') |> #extrae la banda por el nombre, también se podría haber usado el índice (1)
    terra::rast() |> 
    terra::app(mean)
})

lst_oct <- lst_out_oct |> 
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |> # aplica el mosaico para unir los tres tiles
  terra::project("EPSG:32719") |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |>  #aplica mascara para el área de estudio
  terra::trim() |>  #elimina NA sobrantes en los bordes
  terra::app(\(x) x-273.15) #convierte la T° de Kelvin a °C

# resampleo

lst_oct <- resample(lst_oct,r)
names(lst_oct) <- 'lst_oct'
plot(lst_oct, main="LST octubre")

LST Diciembre

lst_out_dic <- purrr::map(tiles,\(tile){
  files_lst_dic <- list.files('data/raw/modis/MOD11A2.v061/diciembre',full.names = TRUE,pattern = tile)
  
  lst_regs_dic <- files_lst_dic |> 
    purrr::map(rast,lyrs = 'LST_Day_1km') |> #extrae la banda por el nombre, también se podría haber usado el índice (1)
    terra::rast() |> 
    terra::app(mean)
})

lst_dic <- lst_out_dic |> 
  terra::sprc() |>  # crea una spatRasterCollection para poder hacer el mosaico
  terra::mosaic() |> # aplica el mosaico para unir los tres tiles
  terra::project("EPSG:32719") |> #reproyecta a UTM
  terra::mask(vect(reg_coq_utm)) |>  #aplica mascara para el área de estudio
  terra::trim() |>  #elimina NA sobrantes en los bordes
  terra::app(\(x) x-273.15) #convierte la T° de Kelvin a °C

# resampleo

lst_dic <- resample(lst_dic,r)
names(lst_dic) <- 'lst_dic'
plot(lst_dic, main="LST diciembre")

LST de los cuatros meses

lst_coq <- c(lst_ene,lst_jun,lst_oct,lst_dic)
plot(lst_coq)

Resamplear a la misma resolución los raster

#resamplea para ajustar a resolución de 500 m

elev_pred_res <- terra::resample(elev_pred,r)
predictores <- c(elev_pred_res,dcosta,ndvi_ene,ndvi_jun,ndvi_oct,ndvi_dic,lst_ene,lst_jun,lst_oct,lst_dic)

Exportar y visualizar

writeRaster(predictores,'data/raw/predictorias/predictores_Coq.tif',overwrite = TRUE) #guarda los predictores
plot(predictores)

6. Datos de estaciones in situ

obtención de datos agromet

estaciones_agromet |> glimpse()

Rows: 417
Columns: 8
$ ema           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
$ institucion   <chr> "FDF", "FDF", "FDF", "FDF", "FDF", "FDF", "FDF", "FDF", …
$ nombre_ema    <chr> "Azapa1", "Azapa2", "Tranque Lautaro", "Jotabeche", "Hor…
$ comuna        <chr> "Arica", "Arica", "Tierra Amarilla", "Tierra Amarilla", …
$ region        <chr> "Arica y Parinacota", "Arica y Parinacota", "Atacama", "…
$ latitud       <dbl> -18.50964, -18.52044, -27.97558, -27.58861, -27.72889, -…
$ longitud      <dbl> -70.24806, -70.23267, -70.00000, -70.24472, -70.19667, -…
$ fecha_de_alta <dttm> 2013-03-08 06:49:10, 2013-03-08 06:49:10, 2013-03-08 06…

estas_sf <- estaciones_agromet |> 
  st_as_sf(coords = c('longitud','latitud'),crs = 4326) |> #crea objeto sf
  st_transform(32719) #reproyecta a UTM WGS84 huso 19 Sur

estas_regs <- st_intersection(estas_sf,reg_coq_utm) #intersección estaciones con área de estudio

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

ids <- estas_regs |> pull(ema) #extrae la columna ema de las estaciones

# leer los datos de agromet del 2023
data_agro <- read_rds('data/raw/agromet/datos_agromet_2023.rds')

data_agro |> glimpse()

Rows: 3,651,647
Columns: 13
$ station_id            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ fecha_hora            <dttm> 2023-01-01 00:00:00, 2023-01-01 01:00:00, 2023-…
$ temp_promedio_aire    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ precipitacion_horaria <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ humed_rel_promedio    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ presion_atmosferica   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ radiacion_solar_max   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ veloc_max_viento      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ temp_minima           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ temp_maxima           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ direccion_del_viento  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ grados_dia            <dbl> 425.75, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ horas_frio            <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

Tratamiento de datos anómalos y faltantes

# filtrado estadistico
filt_out2 <- function(x,c,...){
  M <- .6745*(x - median(x,...))/mad(x,...)
  
  x[M > c] <- NA 
  x[M < -c] <- NA
  
  return(x)
}

data_agro_estaciones <- data_agro |> 
  filter(station_id %in% ids) |> #filtra estaciones que corresponden a ids
  select(station_id,fecha_hora,temp_promedio_aire) |> #selecciona las columnas indicadas  
  mutate(temp_promedio_aire = filt_out2(temp_promedio_aire,c=1.7)) |> #aplica filtro simple para outliers
  mutate(mes = month(fecha_hora)) |> # crea una variable mes
  filter(mes %in% c(1,6,10,12)) |>  #filtra los meses 
  group_by(station_id,mes) |> # agrupa por codigo de estación y mes
  summarize(temp = mean(temp_promedio_aire,na.rm = TRUE)) |> #calcula el promedio mensual
  drop_na() #elimina las filas que tienen valores NA

`summarise()` has grouped output by 'station_id'. You can override using the
`.groups` argument.

data_estas_sf <- estas_regs |> 
  select(ema) |> 
  dplyr::rename(station_id = ema) |> 
  left_join(data_agro_estaciones) 

Joining with `by = join_by(station_id)`

#Asegurarse de eleminar valores NA en temperatura
data_estas_sf <- data_estas_sf |> filter(!is.na(temp))

Rápida revisión de datos con un diagrama de densidad de temperatura

data_estas_sf |> 
  st_drop_geometry() |> 
  ggplot(aes(temp)) + 
  geom_density() +
  theme_bw()

Visualización de la temperatura segun latitud

data_estas_sf |> 
  cbind(st_coordinates(data_estas_sf |> st_transform(4326))) |> 
  st_drop_geometry() |> 
  ggplot(aes(temp,Y)) + 
  geom_point() +
  theme_bw()

Mapa dinámico de los predictores de elevación y costa. El usuario puede ajustar el mes que desea ver. Es importante considerar que los meses a filtrar son 1, 6, 10 y 12

tmap_mode('view')

tmap mode set to interactive viewing

#install.packages("viridis")
library(viridis)

Cargando paquete requerido: viridisLite

preds <- rast('data/raw/predictorias/predictores_Coq.tif')
tm_shape(preds) + 
  tm_raster(style = 'cont') +
  tm_shape(data_estas_sf |> filter(mes == 1)) + # mes de enero, cambiar el numero para otro mes
  tm_dots(col='temp',style = 'jenks',palette = viridis::viridis(10)) +
  tm_facets(free.scales.fill = TRUE,as.layers = FALSE) +
  tm_layout(legend.show = FALSE)

Tip: rasters can be shown as layers instead of facets by setting tm_facets(as.layers = TRUE).

7. Extracción de predictores

Es necesario extraer los valores de los datos puntuales de las estaciones que midieron temperatura con los datos ráster que se obtuvieron recientemente. Esto con el fin de crear una base de datos para luego trabajar con ella en el siguiente taller.

#extrae los datos de los predictores en las ubicaciones de las estaciones
data_temp_pred <- terra::extract(predictores,data_estas_sf)

#une con los datos de temperatura in-situ
data_unida <- cbind(data_estas_sf,data_temp_pred) |> 
  select(-ID) |> 
  drop_na()

#almacena en el disco el resultado
write_rds(data_unida,'data/procesada/data_temp_con_predictores_Coq.rds')

Visualizar tabla dinamica

#crea una tabla dinámica con {gt}

data_unida |> 
  sf::st_drop_geometry() |> 
  gt::gt() |> 
  gt::fmt_number(decimals = 1, drop_trailing_zeros = TRUE) |>
  gt::opt_interactive()

8. Ánalisis estadistico rápido

Se desea visualizar la estadistica de las mediciones de temperatura y de los predictores. Esto se puede hacer de la siguiente manera. Se quita el campo geometria para efecto del funcionamiento de la herramienta estadística

# Crear una copia del dataframe sin la columna de geometría para la estadistica
data_sin_geom <- st_drop_geometry(data_unida)

# Resumen usando summary
summary(data_sin_geom)

   station_id         mes             temp             dem         
 Min.   : 10.0   Min.   : 1.00   Min.   :-3.341   Min.   :  29.37  
 1st Qu.:323.0   1st Qu.: 2.25   1st Qu.:13.279   1st Qu.: 231.94  
 Median :337.0   Median : 8.00   Median :15.723   Median : 482.77  
 Mean   :305.4   Mean   : 7.23   Mean   :15.269   Mean   : 774.28  
 3rd Qu.:353.0   3rd Qu.:10.00   3rd Qu.:18.065   3rd Qu.: 909.05  
 Max.   :664.0   Max.   :12.00   Max.   :22.646   Max.   :4470.75  
     aspect            slope           dist_costa          ndvi_ene      
 Min.   :  9.197   Min.   : 0.3427   Min.   :   600.9   Min.   :0.09522  
 1st Qu.:188.014   1st Qu.: 2.2136   1st Qu.: 19451.6   1st Qu.:0.22780  
 Median :238.823   Median : 3.4919   Median : 40340.6   Median :0.27647  
 Mean   :233.551   Mean   : 3.9659   Mean   : 46821.1   Mean   :0.27617  
 3rd Qu.:304.381   3rd Qu.: 4.6179   3rd Qu.: 70821.4   3rd Qu.:0.34396  
 Max.   :351.015   Max.   :12.1583   Max.   :134478.1   Max.   :0.51612  
    ndvi_jun          ndvi_oct          ndvi_dic          lst_ene     
 Min.   :0.07932   Min.   :0.05976   Min.   :0.08497   Min.   :25.01  
 1st Qu.:0.21925   1st Qu.:0.18249   1st Qu.:0.17173   1st Qu.:33.83  
 Median :0.26012   Median :0.23851   Median :0.23930   Median :36.22  
 Mean   :0.25526   Mean   :0.23897   Mean   :0.23429   Mean   :35.58  
 3rd Qu.:0.31767   3rd Qu.:0.29393   3rd Qu.:0.29898   3rd Qu.:37.91  
 Max.   :0.40071   Max.   :0.49559   Max.   :0.45717   Max.   :42.71  
    lst_jun           lst_oct          lst_dic     
 Min.   :-0.9634   Min.   : 9.574   Min.   :24.05  
 1st Qu.:18.2015   1st Qu.:29.861   1st Qu.:34.21  
 Median :19.5787   Median :32.181   Median :37.65  
 Mean   :18.7973   Mean   :31.053   Mean   :36.31  
 3rd Qu.:21.5643   3rd Qu.:33.910   3rd Qu.:39.03  
 Max.   :23.4239   Max.   :38.357   Max.   :44.41  

Otra forma de ver un resumen con otros estadistico es con describe():

# Revisar resumen detallado con describe
describe(data_sin_geom)

           vars   n     mean       sd   median  trimmed      mad    min
station_id    1 178   305.42   173.69   337.00   300.69    23.72  10.00
mes           2 178     7.23     4.22     8.00     7.40     4.45   1.00
temp          3 178    15.27     4.38    15.72    15.71     3.61  -3.34
dem           4 178   774.28   951.49   482.77   556.45   443.91  29.37
aspect        5 178   233.55    74.87   238.82   239.42    96.00   9.20
slope         6 178     3.97     2.71     3.49     3.55     1.78   0.34
dist_costa    7 178 46821.11 33581.79 40340.59 43874.81 37454.22 600.93
ndvi_ene      8 178     0.28     0.09     0.28     0.28     0.09   0.10
ndvi_jun      9 178     0.26     0.08     0.26     0.26     0.08   0.08
ndvi_oct     10 178     0.24     0.09     0.24     0.24     0.08   0.06
ndvi_dic     11 178     0.23     0.08     0.24     0.23     0.09   0.08
lst_ene      12 178    35.58     3.69    36.22    35.89     3.07  25.01
lst_jun      13 178    18.80     4.35    19.58    19.62     2.83  -0.96
lst_oct      14 178    31.05     5.01    32.18    31.84     3.06   9.57
lst_dic      15 178    36.31     4.25    37.65    36.69     3.08  24.05
                 max     range  skew kurtosis      se
station_id    664.00    654.00 -0.11    -0.10   13.02
mes            12.00     11.00 -0.39    -1.35    0.32
temp           22.65     25.99 -1.54     4.11    0.33
dem          4470.75   4441.39  2.17     4.18   71.32
aspect        351.02    341.82 -0.65     0.12    5.61
slope          12.16     11.82  1.44     1.80    0.20
dist_costa 134478.10 133877.17  0.64    -0.34 2517.06
ndvi_ene        0.52      0.42  0.12    -0.10    0.01
ndvi_jun        0.40      0.32 -0.39    -0.41    0.01
ndvi_oct        0.50      0.44  0.28     0.28    0.01
ndvi_dic        0.46      0.37  0.13    -0.45    0.01
lst_ene        42.71     17.71 -0.90     0.90    0.28
lst_jun        23.42     24.39 -2.30     6.27    0.33
lst_oct        38.36     28.78 -1.89     4.56    0.38
lst_dic        44.41     20.36 -0.86     0.76    0.32

Como los datos tabulado estan todos juntos, se debe separar por mes

data_ene <- data_sin_geom[data_sin_geom$mes == 1, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_ene", "lst_ene")]
data_jun <- data_sin_geom[data_sin_geom$mes == 6, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_jun", "lst_jun")]
data_oct <- data_sin_geom[data_sin_geom$mes == 10, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_oct", "lst_oct")]
data_dic <- data_sin_geom[data_sin_geom$mes == 12, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_dic", "lst_dic")]

Estadística de enero

describe(data_ene)

           vars  n     mean       sd   median  trimmed      mad    min
station_id    1 45   306.04   174.32   337.00   301.38    23.72  10.00
temp          2 45    17.97     3.92    18.59    18.66     2.02   0.00
dem           3 45   791.08   995.19   482.77   577.70   443.91  29.37
aspect        4 45   232.67    75.64   238.82   238.31    96.00   9.20
slope         5 45     4.02     2.78     3.49     3.61     1.78   0.34
dist_costa    6 45 47081.08 34456.06 40340.59 44087.58 37454.22 600.93
ndvi_ene      7 45     0.28     0.09     0.28     0.27     0.09   0.10
lst_ene       8 45    35.50     3.80    36.22    35.81     3.07  25.01
                 max     range  skew kurtosis      se
station_id    664.00    654.00 -0.12    -0.17   25.99
temp           22.65     22.65 -2.54     8.17    0.58
dem          4470.75   4441.39  2.11     3.85  148.35
aspect        351.02    341.82 -0.61    -0.03   11.28
slope          12.16     11.82  1.39     1.55    0.41
dist_costa 134478.10 133877.17  0.64    -0.41 5136.41
ndvi_ene        0.52      0.42  0.09    -0.21    0.01
lst_ene        42.71     17.71 -0.88     0.68    0.57

Histograma temperatura de enero

hist(data_ene$temp, 
     main = "Histograma de Temperaturas en Enero", 
     xlab = "Temperatura (°C)", 
     ylab = "Frecuencia", 
     col = "lightblue", 
     border = "black")

Estadística de junio

describe(data_jun)

           vars  n     mean       sd   median  trimmed      mad    min
station_id    1 44   304.07   175.83   336.00   298.83    24.46  10.00
temp          2 44    12.41     3.62    13.15    12.99     1.66  -3.13
dem           3 44   806.76  1001.05   485.60   590.94   429.22  29.37
aspect        4 44   233.64    76.24   241.51   239.65    95.56   9.20
slope         5 44     4.00     2.81     3.38     3.57     1.79   0.34
dist_costa    6 44 48015.68 34272.57 41218.37 45146.09 36827.50 600.93
ndvi_jun      7 44     0.25     0.08     0.25     0.26     0.08   0.08
lst_jun       8 44    18.71     4.66    19.85    19.56     2.50  -0.96
                 max     range  skew kurtosis      se
station_id    664.00    654.00 -0.09    -0.21   26.51
temp           17.04     20.17 -2.45     7.15    0.55
dem          4470.75   4441.39  2.08     3.70  150.91
aspect        351.02    341.82 -0.64    -0.03   11.49
slope          12.16     11.82  1.40     1.50    0.42
dist_costa 134478.10 133877.17  0.63    -0.41 5166.78
ndvi_jun        0.40      0.32 -0.36    -0.53    0.01
lst_jun        23.42     24.39 -2.27     5.79    0.70

Histograma temperatura de junio

hist(data_jun$temp, 
     main = "Histograma de Temperaturas en Junio", 
     xlab = "Temperatura (°C)", 
     ylab = "Frecuencia", 
     col = "lightblue", 
     border = "black")

Estadística de octubre

describe(data_oct)

           vars  n     mean       sd   median  trimmed      mad    min
station_id    1 45   306.04   174.32   337.00   301.38    23.72  10.00
temp          2 45    13.35     4.19    14.24    14.27     2.00  -3.34
dem           3 45   791.08   995.19   482.77   577.70   443.91  29.37
aspect        4 45   232.67    75.64   238.82   238.31    96.00   9.20
slope         5 45     4.02     2.78     3.49     3.61     1.78   0.34
dist_costa    6 45 47081.08 34456.06 40340.59 44087.58 37454.22 600.93
ndvi_oct      7 45     0.24     0.09     0.24     0.24     0.08   0.06
lst_oct       8 45    30.91     5.29    32.18    31.69     3.06   9.57
                 max     range  skew kurtosis      se
station_id    664.00    654.00 -0.12    -0.17   25.99
temp           17.34     20.68 -2.52     6.24    0.62
dem          4470.75   4441.39  2.11     3.85  148.35
aspect        351.02    341.82 -0.61    -0.03   11.28
slope          12.16     11.82  1.39     1.55    0.41
dist_costa 134478.10 133877.17  0.64    -0.41 5136.41
ndvi_oct        0.50      0.44  0.31     0.20    0.01
lst_oct        38.36     28.78 -1.89     4.33    0.79

Histograma temperatura de octubre

hist(data_oct$temp, 
     main = "Histograma de Temperaturas en Octubre", 
     xlab = "Temperatura (°C)", 
     ylab = "Frecuencia", 
     col = "lightblue", 
     border = "black")

Estadística de diciembre

describe(data_dic)

           vars  n     mean       sd   median  trimmed      mad    min
station_id    1 44   305.50   176.30   337.50   300.58    23.72  10.00
temp          2 44    17.33     2.82    17.15    17.59     2.32   8.51
dem           3 44   707.45   831.50   431.71   521.79   369.11  29.37
aspect        4 44   235.27    74.47   241.51   240.96    86.26   9.20
slope         5 44     3.83     2.52     3.38     3.47     1.79   0.34
dist_costa    6 44 45094.79 32142.79 40315.82 42603.53 35489.38 600.93
ndvi_dic      7 44     0.24     0.08     0.24     0.24     0.09   0.08
lst_dic       8 44    36.50     4.00    37.65    36.79     2.83  24.97
                 max     range  skew kurtosis      se
station_id    664.00    654.00 -0.11    -0.24   26.58
temp           22.30     13.79 -0.96     1.50    0.43
dem          3447.35   3417.99  2.03     3.41  125.35
aspect        351.02    341.82 -0.67     0.19   11.23
slope          11.41     11.07  1.37     1.83    0.38
dist_costa 124083.46 123482.53  0.53    -0.63 4845.71
ndvi_dic        0.46      0.37  0.13    -0.49    0.01
lst_dic        44.41     19.44 -0.69     0.45    0.60

Histograma temperatura de diciembre

hist(data_dic$temp, 
     main = "Histograma de Temperaturas en Diciembre", 
     xlab = "Temperatura (°C)", 
     ylab = "Frecuencia", 
     col = "lightblue", 
     border = "black")

Exportar los datos en un csv

# Crear una lista con los nombres de los dataframes y los nombres de los archivos CSV
meses <- list("ene" = data_ene, "jun" = data_jun, "oct" = data_oct, "dic" = data_dic)
ruta <- "data/procesada/"

# Usar un bucle para calcular el resumen, convertir a data.frame y guardar el CSV
for (mes in names(meses)) {
  # Calcular el resumen usando describe
  resumen <- describe(meses[[mes]])
  
  # Convertir a data.frame
  resumen_df <- as.data.frame(resumen)
  
  # Escribir el resumen como archivo CSV
  write.csv(resumen_df, paste0(ruta, "resumen_describe_", mes, "_Coq.csv"), row.names = TRUE)
}

# Fin del Script