Taller 3: Análisis de autocorrelación espacial de la variable temperatura en la Región de Coquimbo

Author

Victor Contreras, Cristian Aguirre y Christian Chacon

Published

October 12, 2023

Taller 3: Análisis de autocorrelación espacial

1. Introducción

En este laboratorio se desarrollan procesos para examinar la autocorrelación espacial de los datos de temperatura en la Región de Coquimbo durante los meses de enero, junio, octubre y diciembre de 2023. Se emplean tanto el Índice de Morán Global como el Índice de Morán Local, los cuales permiten identificar y evaluar patrones espaciales de correlación en las temperaturas registradas. Según Anselin (1995), estas herramientas son fundamentales para detectar conglomerados espaciales o zonas de dispersión en los datos geográficos. Además, se complementa el análisis con el uso de variogramas, los cuales modelan la relación espacial entre puntos geográficos, permitiendo observar cómo varía una variable, en este caso la temperatura, en función de la distancia entre estaciones de medición (ESRI, s.f.). Este método es clave para evaluar la estructura espacial de los datos, ya que proporciona una representación gráfica que ayuda a identificar patrones de autocorrelación a diferentes escalas espaciales (ArcGIS Pro, s.f.). Estos enfoques de autocorrelación espacial permiten captar las dinámicas térmicas en relación con la influencia geográfica de las zonas costeras y montañosas, además de observar las diferencias estacionales, lo que resulta fundamental para comprender la distribución térmica a lo largo del año y su relación con la geografía de la región

2. Objetivos

1.- Realizar el análisis de autocorrelación espacial global calculando el índice de Moran global para los datos de temperatura de los meses de enero, junio, octubre y diciembre en la región de Coquimbo. Aplicar la prueba de Montecarlo para determinar la significancia estadística de la autocorrelación espacial y generar los gráficos correspondientes para su interpretación. 2.- Evaluar la autocorrelación espacial local mediante el cálculo del índice de Moran local para los mismos meses, generando mapas con el paquete {tmap} que muestren los valores presentes en los cuadrantes de autocorrelación espacial (high-high, high-low, low-low, low-high). Comparar estos resultados con el análisis global e interpretar las diferencias. 3.- Analizar la estructura espacial a través de variogramas experimentales para los datos de temperatura mensual, ajustando modelos de variograma y evaluando la anisotropía espacial. Almacenar los modelos ajustados en archivos .rds, describiendo la estructura espacial y los resultados obtenidos. 4.- A partir del modelo no geoestadístico con menor RMSE y mayor R² del taller 2, obtener el variograma de sus residuos.

3. Desarrollo

3.1. Preparación de datos

library(tidyverse)

── 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(terra)

terra 1.7.78

Adjuntando el paquete: 'terra'

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

    extract

library(tmap)

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

library(sf)

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

library(viridis)

Cargando paquete requerido: viridisLite

library(viridisLite)
library(spData)

To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

library(raster)

Cargando paquete requerido: sp

Adjuntando el paquete: 'raster'

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

    select

library(purrr)
library(gstat)
library(stars)

Cargando paquete requerido: abind

library(raster)
library(ggplot2)
library(patchwork)


Adjuntando el paquete: 'patchwork'

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

    area

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

    area

library(spdep)
library(units)

udunits database from C:/Users/criaguirre/AppData/Local/Programs/R/R-4.4.1/library/units/share/udunits/udunits2.xml

library(sp)
library(abind)
library(gt)

# Directorio
setwd("~/00_tallegeo/taller01_geost")

# Base de datos
data_df <- read_rds('data/procesada/data_temp_con_predictores_Coq.rds') # vectorial puntos
predictorias_ras <- rast('data/raw/predictorias/predictores_Coq.tif') # raster
reg_coq_utm <- read_rds('data/raw/region_coquimbo.rds') # vectorial poligono
reg_coq_utm <- reg_coq_utm |> 
  st_transform(32719)
reg_coq_geo <- reg_coq_utm |> 
  st_transform(4326)

# Filtrar datos vectoriales por mes

data_ene_df <- data_df[data_df$mes == 1, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_ene", "lst_ene", "geometry")]
data_jun_df <- data_df[data_df$mes == 6, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_jun", "lst_jun", "geometry")]
data_oct_df <- data_df[data_df$mes == 10, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_oct", "lst_oct", "geometry")]
data_dic_df <- data_df[data_df$mes == 12, c("station_id", "temp", "dem", "aspect", "slope", "dist_costa", "ndvi_dic", "lst_dic", "geometry")]

# filtrar datos raster por mes

pred_ene_ras <- subset(predictorias_ras, c("dem", "aspect", "slope", "dist_costa", "ndvi_ene", "lst_ene"))
pred_jun_ras <- subset(predictorias_ras, c("dem", "aspect", "slope", "dist_costa", "ndvi_jun", "lst_jun"))
pred_oct_ras <- subset(predictorias_ras, c("dem", "aspect", "slope", "dist_costa", "ndvi_oct", "lst_oct"))
pred_dic_ras <- subset(predictorias_ras, c("dem", "aspect", "slope", "dist_costa", "ndvi_dic", "lst_dic"))

3.2. Indice de Moran Global

Calculo de distancia y eliminación de vecinos cercanos, con el fin evitar problemas de auocorrelación espacial a corta distancia.

# Enero
d_ene <- st_distance(data_ene_df) #Creación de matriz de distancia euclidiana entre dos puntos
diag(d_ene) <- NA # Establece en NA los valores de la diagonal de la matriz de distancia
d_ene <- units::drop_units(d_ene) # Elimina las unidades en metro de la matriz de distancia
ind <- which(d_ene < 300, arr.ind = TRUE) # Encuentra los indices de distancia menores a 300 m
data_ene_df <- data_ene_df[-ind[,1],] # Elimina los puntos menores a 300 m
#elimina dos estaciones

#Junio
d_jun <- st_distance(data_jun_df)
diag(d_jun) <- NA
d_jun <- units::drop_units(d_jun)
ind <- which(d_jun < 300, arr.ind = TRUE)
data_jun_df <- data_jun_df[-ind[,1],]
#elimina dos estaciones

#Octubre
d_oct <- st_distance(data_oct_df)
diag(d_oct) <- NA
d_oct <- units::drop_units(d_oct)
ind <- which(d_oct < 300, arr.ind = TRUE)
data_oct_df <- data_oct_df[-ind[,1],]
#elimina dos estaciones

#Diciembre
d_dic <- st_distance(data_dic_df)
diag(d_dic) <- NA
d_dic <- units::drop_units(d_dic)
ind <- which(d_dic < 300, arr.ind = TRUE)
data_dic_df <- data_dic_df[-ind[,1],]
#elimina dos estaciones

Creación de pesos espaciales para los indices de Moran

# Enero
w_ene <- 1/as.matrix(st_distance(data_ene_df)) # Crea una matriz de pesos esapciales inversamente proporcional a la distancia
diag(w_ene) <- 0 # Establece los valores de la diagonal en cero (sin autocorrelacion entre puntos consigo mismos)
w_ene <- units::drop_units(w_ene) # Elimina las unidades d ela matriz de distancia, dejándola como una matriz numérica pura
w_ene <- mat2listw(w_ene,  style = "W") # Convierte la matriz de pesos espaciales en un objeto de lista de pesos. Para este tipo de datos se utiliza el estilo W, que es para filas estandarizadas, donde los pesos suman 1 en cada fila

# Junio
w_jun <- 1/as.matrix(st_distance(data_jun_df))
diag(w_jun) <- 0
w_jun <- units::drop_units(w_jun)
w_jun <- mat2listw(w_jun,  style = "W")

# Octubre
w_oct <- 1/as.matrix(st_distance(data_oct_df))
diag(w_oct) <- 0
w_oct <- units::drop_units(w_oct)
w_oct <- mat2listw(w_oct,  style = "W")

# Diciembre
w_dic <- 1/as.matrix(st_distance(data_dic_df))
diag(w_dic) <- 0
w_dic <- units::drop_units(w_dic)
w_dic <- mat2listw(w_dic,  style = "W")

Índice Global de Moran. Este mide la autocorralación espacial en todo el conjunto de datos. Este indice envaúa si las temperaturas de las estaciones cercanas tienden a ser similares (autocorrelación positivo o diferente) o diferentes (autocorrelación negativa)

# Enero
moran_ene <- moran(data_ene_df$temp, w_ene , n = length(w_ene$neighbours), S0 = Szero(w_ene))# estimación indice de moran
graf_moran_ene <- moran.plot(data_ene_df$temp,w_ene)
title(main = "Moran's I para Enero", col.main = "blue", font.main = 4) #Ploteo grafico auto-corr Moran

# Junio
moran_jun <- moran(data_jun_df$temp, w_jun , n = length(w_jun$neighbours), S0 = Szero(w_jun))# estimación indice de moran
graf_moran_jun <- moran.plot(data_jun_df$temp,w_jun)
title(main = "Moran's I para Junio", col.main = "blue", font.main = 4) #Ploteo grafico auto-corr Moran

# Octubre
moran_oct <- moran(data_oct_df$temp, w_oct , n = length(w_oct$neighbours), S0 = Szero(w_oct))# estimación indice de moran
graf_moran_oct <- moran.plot(data_oct_df$temp,w_oct) 
title(main = "Moran's I para Octubre", col.main = "blue", font.main = 4) #Ploteo grafico auto-corr Moran

# Diciembre
moran_oct <- moran(data_dic_df$temp, w_dic , n = length(w_dic$neighbours), S0 = Szero(w_dic))# estimación indice de moran
graf_moran_dic <- moran.plot(data_dic_df$temp,w_dic) 
title(main = "Moran's I para Diciembre", col.main = "blue", font.main = 4)  # Título principal

#Ploteo grafico auto-corr Moran

Indice de Moran Global para enero

La línea de tendencia muestra una pendiente positiva entre los valores de temperatura y sus valores espacialmente retardados (lagged). Esto sugiere que, en general, durante el mes de enero, las temperaturas altas tienden a agruparse con otras temperaturas altas en áreas cercanas. Sin embargo, la pendiente es leve, lo que podría indicar que la correlación no es particularmente fuerte.

En el cuadrante superior derecho (High-High), se encuentran valores de temperatura altos y superiores al promedio, donde se observa una notable agrupación de puntos con mediciones elevadas. Por otro lado, en el cuadrante inferior izquierdo (Low-Low), hay agrupaciones de temperaturas bajas con otras bajas, lo cual podría estar relacionado con zonas de la cordillera principal y áreas costeras. La distribución de puntos en estos dos cuadrantes contribuye a la pendiente positiva de la línea de tendencia, aunque de manera moderada.

El cuadrante superior izquierdo (Low-High) representa puntos con temperaturas bajas (por debajo del promedio) que están rodeados de vecinos con temperaturas altas, lo cual podría reflejar transiciones abruptas de temperatura, probablemente debido a cambios bruscos en la topografía.

En el cuadrante inferior derecho (High-Low), se observan puntos con temperaturas altas rodeados de vecinos con temperaturas bajas, lo que provoca una ligera disminución de la pendiente de la línea de tendencia. Esto podría estar relacionado con la influencia de factores locales como la vegetación, la exposición solar o condiciones geográficas específicas.

Además, hay algunos valores destacados, como los días 17, 18, 23 y 31, que podrían ser considerados como valores atípicos debido a su contraste con el entorno.

En general, enero muestra una ligera autocorrelación positiva, lo que indica que las temperaturas en la región tienden a ser similares a las de sus localidades vecinas, salvo algunos puntos donde la influencia de factores geográficos, la proximidad a la costa o la exposición solar podrían modificar este patrón.

Indice de Moran para Junio

A diferencia del mes anterior, la temperatura promedio entre lo observado y retardado baja considerablemente pasando de ~18 °C en enero a ~13 °C. La distribución de los puntos y la linea de tendencia sugiere una autocorrelacion positiva, es decir, que las áresa con temperaturas sumilares tienden a agruparse. La presencia de valores en los cuadrantes II (Low-High) y (High-Low), al igual que en enero, presentan cambios de temperatura abruptos probablemente ligado a cambios abruptos en zonas muy localizadas. Hay presencia de valores atípicos en los cuadrantes I (High-High), II y en el III (Low-Low), siendo este último probablemente relacionado con zonas de alta montaña.

Indice de Moran para Octubre

Se aprecia un aumento de la tempertaura promedio, relacionado con el cambio de estacionalidad, y la distribución de los puntos del cuadrante I y III sugiere una autocorrelación espacial positiva. La menor concentración de puntos en los cuadrantes II y IV indica que hay menos transiciones abruptas entre áreas de diferente temperatura, lo que podría sugerir una relativa homogeneidad en la distribución espacial de la temperatura durante octubre.

Indice de Moran para Diciembre Diciembre presenta una autocorrelación espacial positiva en la distribución de las temperaturas, lo cual se refleja en la pendiente positiva de la línea de ajuste.

La mayoría de los puntos se concentran en los cuadrantes I y III, sugiriendo la agrupación de temperaturas similares (altas y bajas) en determinadas regiones.

La dispersión de algunos puntos en los cuadrantes II y IV indica áreas de transición entre diferentes temperaturas, aunque estas son menos comunes. Esto puede reflejar condiciones locales específicas que afectan las temperaturas de manera desigual.

A grandes rasgos, los gráficos de dispersión complementan el análisis del Índice de Morán para las cuatro estaciones, mostrando una autocorrelación espacial positiva. Sin embargo, también se observan algunos cambios abruptos, posiblemente relacionados con variaciones en las condiciones geográficas, lo que provoca que la pendiente no sea lo suficientemente pronunciada como para evidenciar una fuerte autocorrelación.

A continuación se muestran los cuadro gráficos ploteados en una misma figura

# Crear el layout para 4 gráficos en una figura
par(mfrow = c(2, 2))  # 2 filas, 2 columnas

# Graficar los 4 gráficos de Moran
moran.plot(data_ene_df$temp, w_ene, main = "Moran's I - Enero")
moran.plot(data_jun_df$temp, w_jun, main = "Moran's I - Junio")
moran.plot(data_oct_df$temp, w_oct, main = "Moran's I - Octubre")
moran.plot(data_dic_df$temp, w_dic, main = "Moran's I - Diciembre")

# Restablecer la configuración original de gráficos (opcional)
par(mfrow = c(1, 1))

3.3. Indice de Moran Local

Los resultados de las simulaciones de Monte-Carlo proprociona una evaluación de autocorrelacion espacial de la temepratura para cada uno de los meses en cuestión. El valor de Moran I indica el grado de autocorrelacion espacial de los datos, por lo que si entrega un valor positiva sugiere una agrupación de valores similares, caso contrario indicaria una distribución aleatoria. El P-value evalúa la significancia estadistica del indice de Moran, un valor p bajo (<0.05) sugiere una autocorrelacion observada y que esto no es fruto del azar, por ende, es significativo.

# Enero
imc_moran_ene <- moran.mc(data_ene_df$temp,w_ene,nsim=10000) 
imc_moran_ene


    Monte-Carlo simulation of Moran I

data:  data_ene_df$temp 
weights: w_ene  
number of simulations + 1: 10001 

statistic = 0.081411, observed rank = 9847, p-value = 0.0154
alternative hypothesis: greater

# Junio
imc_moran_jun <- moran.mc(data_jun_df$temp,w_jun,nsim=10000) 
imc_moran_jun


    Monte-Carlo simulation of Moran I

data:  data_jun_df$temp 
weights: w_jun  
number of simulations + 1: 10001 

statistic = 0.098219, observed rank = 9895, p-value = 0.0106
alternative hypothesis: greater

# Octubre
imc_moran_oct <- moran.mc(data_oct_df$temp,w_oct,nsim=10000) 
imc_moran_oct


    Monte-Carlo simulation of Moran I

data:  data_oct_df$temp 
weights: w_oct  
number of simulations + 1: 10001 

statistic = 0.078678, observed rank = 9824, p-value = 0.0177
alternative hypothesis: greater

# Diciembre
imc_moran_dic <- moran.mc(data_dic_df$temp,w_dic,nsim=10000) 
imc_moran_dic


    Monte-Carlo simulation of Moran I

data:  data_dic_df$temp 
weights: w_dic  
number of simulations + 1: 10001 

statistic = 0.09922, observed rank = 9938, p-value = 0.006299
alternative hypothesis: greater

De acuerdo con los resultados se puede decir lo siguiente

Enero

Moran I = 0.081411: Muestra una autocorrelación espacial positiva, lo que significa que las temperaturas similares tienden a agruparse. P-value = 0.0131: Este valor sugiere que la autocorrelación observada es significativa

Junio

Moran I = 0.098219: Refleja una autocorrelación espacial positiva, un poco más fuerte que en enero, lo que sugiere que las temperaturas similares se agrupan con mayor tendencia. P-value = 0.009299: La autocorrelación es significativa, reforzando la idea de que hay una estructura espacial en cómo se distribuyen las temperaturas.

Octubre

Moran I = 0.078678: La autocorrelación positiva es un poco menor que en junio, lo que indica que las temperaturas similares tienden a agruparse, pero con una fuerza ligeramente inferior. P-value = 0.0166: Al ser menor a 0.05, indica que la agrupación de temperaturas no es un resultado al azar, sino que hay un patrón real

**Diciembre*:**

Moran I = 0.09922: Refleja una autocorrelación espacial positiva. P-value = 0.007199: La significancia de este valor confirma que la distribución de las temperaturas no es aleatoria.

En general, durante todos los meses analizados, los valores de Moran I son positivos, lo que sugiere una tendencia a la autocorrelación espacial positiva: las temperaturas similares tienden a agruparse en el espacio. Además, los p-valores menores a 0.05 en todos los casos indican que esta tendencia es estadísticamente significativa y no producto del azar.

Indice Local de Moran, tambien conocido como LISA, mide la autocorrelación esoacial en un nivel local. A difrencia del índice global, este mide la autocorrelacion para cada punto en el espacio

Ii: El valor local de Moran.
Z.Ii: Valor estandarizado de la autocorrelación para ese punto.
Pr(z != E(Ii)): El p-value para evaluar si la autocorrelación en ese punto es significativa.

# Enero
lmoran_ene <- localmoran(data_ene_df$temp, w_ene)

# Junio
lmoran_jun <- localmoran(data_jun_df$temp, w_jun)


# Octubre
lmoran_oct <- localmoran(data_oct_df$temp, w_oct)

# Diciembre
lmoran_dic <- localmoran(data_dic_df$temp, w_dic)

Incorporar los valores de Indice de Moran local a los dataframe de cada mes.

# Enero
data_ene_df$Ii <- lmoran_ene[,'Ii']
data_ene_df$p <- lmoran_ene[,5]
data_ene_df$quad <- attr(lmoran_ene,'quad')$median

# Junio
data_jun_df$Ii <- lmoran_jun[,'Ii']
data_jun_df$p <- lmoran_jun[,5]
data_jun_df$quad <- attr(lmoran_jun,'quad')$median

# Octubre
data_oct_df$Ii <- lmoran_oct[,'Ii']
data_oct_df$p <- lmoran_oct[,5]
data_oct_df$quad <- attr(lmoran_oct,'quad')$median

# Diciembre

data_dic_df$Ii <- lmoran_dic[,'Ii']
data_dic_df$p <- lmoran_dic[,5]
data_dic_df$quad <- attr(lmoran_dic,'quad')$median

A continuación de generan los mapas de autocorrelación local.

Mapas de autocorrelación

Preparación de los valores

# Enero

data_ene_df <- data_ene_df |> 
  mutate(quad_2 = case_when(
    p > 0.05 ~ 'No significativo',
    .default = quad))

# Junio
data_jun_df <- data_jun_df |> 
  mutate(quad_2 = case_when(
    p > 0.05 ~ 'No significativo',
    .default = quad))

# Octubre
data_oct_df <- data_oct_df |> 
  mutate(quad_2 = case_when(
    p > 0.05 ~ 'No significativo',
    .default = quad))

# Diciembre
data_dic_df <- data_dic_df |> 
  mutate(quad_2 = case_when(
    p > 0.05 ~ 'No significativo',
    .default = quad))

Confección de mapa para Enero

# Mapa 0: Temperatura Enero
m0_ene <- tm_shape(data_ene_df) +
  tm_dots(col = "temp", 
          title = "Temperatura Enero", 
          style = "quantile",
          palette = viridis::magma(4),
          midpoint = NA, 
          size = 0.1) +  # Ajusta el tamaño si es necesario
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Temperatura Estaciones - Enero",
            legend.outside = TRUE)

# Mapa 1: Índice de Moran Local
m1_ene <- tm_shape(data_ene_df) +
  tm_dots(col = 'Ii', 
          style = 'quantile', 
          midpoint = NA, 
          size = 0.1, 
          title = 'Índice de Moran Local') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Índice de Moran Local - Enero",
            legend.outside = TRUE) 

# Mapa 2: p-values de Moran Local - Enero
m2_ene <- tm_shape(data_ene_df) +
  tm_dots(col = 'p', 
          midpoint = NA,
          breaks = c(-Inf, 0.05, 1), 
          size = 0.1, 
          title = 'P-score') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "P-valores Moran Local - Enero",
            legend.outside = TRUE) 

# Mapa 3: Cuadrantes de Moran Local - Enero
m3_ene <- tm_shape(data_ene_df) +
  tm_dots(col = 'quad_2', 
          palette = viridis::viridis(5), 
          size = 0.1, 
          title = 'Cuadrantes') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Cuadrantes Moran Local - Enero",
            legend.outside = TRUE)

NOTA: En el renderizado se tuvo problemas con el modo de mapa view (interactivo), por lo que se deja por defecto como plot, aunque si lo desea, puede cambiar a view

tmap_mode('plot')

tmap mode set to plotting

tmap_arrange(m0_ene, m1_ene, m2_ene, m3_ene)

Legend labels were too wide. The labels have been resized to 0.27, 0.25, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.23, 0.24, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.24, 0.30. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.38, 0.40, 0.24. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Mapas Enero formato plot

tmap_mode('plot')

tmap mode set to plotting

m0_ene2 <- tm_shape(data_ene_df,bbox = reg_coq_geo ) +
  tm_dots(col = 'temp', 
          style = 'quantile', 
          palette = viridis::magma(4), size = 0.4) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Temp. Estaciones - Enero",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m1_ene2 <- tm_shape(data_ene_df, bbox = reg_coq_geo) +
  tm_dots(col = 'Ii', style = 'quantile', midpoint = NA, size = 0.4, title = 'Índice de Moran Local' ) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. de Moran Local (Ii) - Enero",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m2_ene2 <- tm_shape(data_ene_df, bbox = reg_coq_geo) +
  tm_dots(col = 'p', midpoint = NA,
          breaks = c(-Inf, 0.05, 1), size = 0.4, title = 'P-valores Moran Local') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. Moran Local p<0.05 - Enero",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

m3_ene2 <- tm_shape(data_ene_df, bbox = reg_coq_geo) +
  tm_dots(col = 'quad_2', palette = viridis::viridis(5), size = 0.4, title = 'Cuadrantes') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Cuadrantes I. Moran Local - Enero",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

tmap_arrange(m0_ene2, m1_ene2, m2_ene2, m3_ene2, nrow=2, ncol=2)

Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
deprecated. It might return a CRS with a non-EPSG compliant axis order.

Legend labels were too wide. The labels have been resized to 0.61, 0.56, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.51, 0.53, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Exportar mapas

tmap_save(
  tmap_arrange(m0_ene2, m1_ene2, m2_ene2, m3_ene2,nrow=2, ncol=2), 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapas_enero_ac_completo.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapas_enero_ac_completo.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m0_ene2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_temp_enero.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_temp_enero.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m1_ene2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_IMoranloc_enero.png", 
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_IMoranloc_enero.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m2_ene2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_p-values_enero.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_p-values_enero.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m3_ene2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_quad_enero.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_quad_enero.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

Mapas de Junio

# Mapa 0: Temperatura Junio
m0_jun <- tm_shape(data_jun_df) +
  tm_dots(col = "temp", 
          title = "Temperatura Junio", 
          style = "quantile", 
          midpoint = NA,
          palette = viridis::magma(4),
          size = 0.1) +  
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Temperatura Estaciones - Junio",
            legend.outside = TRUE)

# Mapa 1: Índice de Moran Local
m1_jun <- tm_shape(data_jun_df) +
  tm_dots(col = 'Ii', 
          style = 'quantile', 
          midpoint = NA, 
          size = 0.1, 
          title = 'Índice de Moran Local') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Índice de Moran Local - Junio",
            legend.outside = TRUE) 

# Mapa 2: p-values de Moran Local - Enero
m2_jun <- tm_shape(data_jun_df) +
  tm_dots(col = 'p', 
          midpoint = NA,
          breaks = c(-Inf, 0.05, 1), 
          size = 0.1, 
          title = 'P-score') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "P-valores Moran Local - Junio",
            legend.outside = TRUE) 

# Mapa 3: Cuadrantes de Moran Local - Enero
m3_jun<- tm_shape(data_jun_df) +
  tm_dots(col = 'quad_2', 
          palette = viridis::viridis(5), 
          size = 0.1, 
          title = 'Cuadrantes') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Cuadrantes Moran Local - Junio",
            legend.outside = TRUE)

Visualización modo interactivo mes de Junio

tmap_mode('plot')

tmap mode set to plotting

tmap_arrange(m0_jun, m1_jun, m2_jun, m3_jun)

Legend labels were too wide. The labels have been resized to 0.26, 0.25, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.23, 0.24, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.24, 0.30. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.40, 0.24. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Mapas Junio formato plot

tmap_mode('plot')

tmap mode set to plotting

m0_jun2 <- tm_shape(data_jun_df,bbox = reg_coq_geo ) +
  tm_dots(col = 'temp', 
          style = 'quantile', 
          palette = viridis::magma(4), size = 0.4) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Temp. Estaciones - Junio",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m1_jun2 <- tm_shape(data_jun_df, bbox = reg_coq_geo) +
  tm_dots(col = 'Ii', style = 'quantile', midpoint = NA, size = 0.4, title = 'Índice de Moran Local' ) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. de Moran Local (Ii) - Junio",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m2_jun2 <- tm_shape(data_jun_df, bbox = reg_coq_geo) +
  tm_dots(col = 'p', midpoint = NA,
          breaks = c(-Inf, 0.05, 1), size = 0.4, title = 'P-valores Moran Local') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. Moran Local p<0.05 - Junio",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

m3_jun2 <- tm_shape(data_jun_df, bbox = reg_coq_geo) +
  tm_dots(col = 'quad_2', palette = viridis::viridis(5), size = 0.4, title = 'Cuadrantes') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Cuadrantes I. Moran Local - Junio",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

tmap_arrange(m0_jun2, m1_jun2, m2_jun2, m3_jun2, nrow=2, ncol=2)

Legend labels were too wide. The labels have been resized to 0.58, 0.56, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.51, 0.53, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Exportar mapas

tmap_save(
  tmap_arrange(m0_jun2, m1_jun2, m2_jun2, m3_jun2,nrow=2, ncol=2), 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapas_junio_ac_completo.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapas_junio_ac_completo.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m0_jun2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_temp_junio.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_temp_junio.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m1_jun2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_IMoranloc_junio.png", 
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_IMoranloc_junio.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m2_jun2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_p-values_junio.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_p-values_junio.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m3_jun2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_quad_junio.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_quad_junio.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

Mapas de Octubre

# Mapa 0: Temperatura Octubre
m0_oct <- tm_shape(data_oct_df) +
  tm_dots(col = "temp", 
          title = "Temperatura Octubre", 
          style = "quantile", 
          midpoint = NA,
          palette = viridis::magma(4),
          size = 0.1) +  
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Temperatura Estaciones - Octubre",
            legend.outside = TRUE)

# Mapa 1: Índice de Moran Local
m1_oct <- tm_shape(data_oct_df) +
  tm_dots(col = 'Ii', 
          style = 'quantile', 
          midpoint = NA, 
          size = 0.1, 
          title = 'Índice de Moran Local') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Índice de Moran Local - Octubre",
            legend.outside = TRUE) 

# Mapa 2: p-values de Moran Local - Octubre
m2_oct <- tm_shape(data_oct_df) +
  tm_dots(col = 'p', 
          midpoint = NA,
          breaks = c(-Inf, 0.05, 1), 
          size = 0.1, 
          title = 'P-score') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "P-valores Moran Local - Octubre",
            legend.outside = TRUE) 

# Mapa 3: Cuadrantes de Moran Local 
m3_oct<- tm_shape(data_oct_df) +
  tm_dots(col = 'quad_2', 
          palette = viridis::viridis(5), 
          size = 0.1, 
          title = 'Cuadrantes') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Cuadrantes Moran Local - Octubre",
            legend.outside = TRUE)

Visualización interactiva

tmap_mode('plot')

tmap mode set to plotting

tmap_arrange(m0_oct, m1_oct, m2_oct, m3_oct)

Legend labels were too wide. The labels have been resized to 0.26, 0.25, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.23, 0.24, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.24, 0.30. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.40, 0.24. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Formato plot

tmap_mode('plot')

tmap mode set to plotting

m0_oct2 <- tm_shape(data_oct_df,bbox = reg_coq_geo ) +
  tm_dots(col = 'temp', 
          style = 'quantile', 
          palette = viridis::magma(4), size = 0.4) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Temp. Estaciones - Octubre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m1_oct2 <- tm_shape(data_oct_df, bbox = reg_coq_geo) +
  tm_dots(col = 'Ii', style = 'quantile', midpoint = NA, size = 0.4, title = 'Índice de Moran Local' ) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. de Moran Local (Ii) - Octubre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m2_oct2 <- tm_shape(data_oct_df, bbox = reg_coq_geo) +
  tm_dots(col = 'p', midpoint = NA,
          breaks = c(-Inf, 0.05, 1), size = 0.4, title = 'P-valores Moran Local') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. Moran Local p<0.05 - Octubre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

m3_oct2 <- tm_shape(data_oct_df, bbox = reg_coq_geo) +
  tm_dots(col = 'quad_2', palette = viridis::viridis(5), size = 0.4, title = 'Cuadrantes') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Cuadrantes I. Moran Local - Octubre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

tmap_arrange(m0_oct2, m1_oct2, m2_oct2, m3_oct2, nrow=2, ncol=2)

Legend labels were too wide. The labels have been resized to 0.58, 0.56, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.51, 0.53, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Guardar mapas de octubre

tmap_save(
  tmap_arrange(m0_oct2, m1_oct2, m2_oct2, m3_oct2,nrow=2, ncol=2), 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapas_octubre_ac_completo.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapas_octubre_ac_completo.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m0_oct2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_temp_octubre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_temp_octubre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m1_oct2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_IMoranloc_octubre.png", 
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_IMoranloc_octubre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m2_oct2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_p-values_octubre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_p-values_octubre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m3_oct2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_quad_octubre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_quad_octubre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

Mapas de diciembre

# Mapa 0: Temperatura Diciembre
m0_dic <- tm_shape(data_dic_df) +
  tm_dots(col = "temp", 
          title = "Temperatura Diciembre", 
          style = "quantile", 
          midpoint = NA,
          palette = viridis::magma(4),
          size = 0.1) +  
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Temperatura Estaciones - Diciembre",
            legend.outside = TRUE)

# Mapa 1: Índice de Moran Local
m1_dic <- tm_shape(data_dic_df) +
  tm_dots(col = 'Ii', 
          style = 'quantile', 
          midpoint = NA, 
          size = 0.1, 
          title = 'Índice de Moran Local') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Índice de Moran Local - Diciembre",
            legend.outside = TRUE) 

# Mapa 2: p-values de Moran Local 
m2_dic <- tm_shape(data_dic_df) +
  tm_dots(col = 'p', 
          midpoint = NA,
          breaks = c(-Inf, 0.05, 1), 
          size = 0.1, 
          title = 'P-score') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "P-valores Moran Local - Diciembre",
            legend.outside = TRUE) 

# Mapa 3: Cuadrantes de Moran Local 
m3_dic<- tm_shape(data_dic_df) +
  tm_dots(col = 'quad_2', 
          palette = viridis::viridis(5), 
          size = 0.1, 
          title = 'Cuadrantes') +
  tm_shape(reg_coq_utm) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_layout(main.title = "Cuadrantes Moran Local - Diciembre",
            legend.outside = TRUE)

Visualización interactiva

tmap_mode('plot')

tmap mode set to plotting

tmap_arrange(m0_dic, m1_dic, m2_dic, m3_dic)

Legend labels were too wide. The labels have been resized to 0.27, 0.25, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.23, 0.24, 0.25, 0.25, 0.25. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.24, 0.30. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.36, 0.38, 0.24. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Modo plot

tmap_mode('plot')

tmap mode set to plotting

m0_dic2 <- tm_shape(data_dic_df,bbox = reg_coq_geo ) +
  tm_dots(col = 'temp', 
          style = 'quantile', 
          palette = viridis::magma(4), size = 0.4) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Temp. Estaciones - Diciembre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m1_dic2 <- tm_shape(data_dic_df, bbox = reg_coq_geo) +
  tm_dots(col = 'Ii', style = 'quantile', midpoint = NA, size = 0.4, title = 'Índice de Moran Local' ) +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. de Moran Local (Ii) - Diciembre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)
  
m2_dic2 <- tm_shape(data_dic_df, bbox = reg_coq_geo) +
  tm_dots(col = 'p', midpoint = NA,
          breaks = c(-Inf, 0.05, 1), size = 0.4, title = 'P-valores Moran Local') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "I. Moran Local p<0.05 - Diciembre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

m3_dic2 <- tm_shape(data_dic_df, bbox = reg_coq_geo) +
  tm_dots(col = 'quad_2', palette = viridis::viridis(5), size = 0.4, title = 'Cuadrantes') +
  tm_shape(reg_coq_geo) +
  tm_borders(col = "black", alpha = 0.65, lwd = 1.5) +
  tm_fill(col= 'grey', alpha=0.3 )+
  tm_layout(main.title = "Cuadrantes I. Moran Local - Diciembre",
            legend.show = TRUE, 
            legend.outside = TRUE,
            legend.bg.color = "white",           
            legend.frame = TRUE) +
  tm_compass(size = 5, type = "arrow", position = c("right", "top"), text.size = 0.4) +
  tm_grid(projection = "+init=epsg:4326",
          lines = TRUE,
          col = "gray",
          alpha = 0.5)

tmap_arrange(m0_dic2, m1_dic2, m2_dic2, m3_dic2, nrow=2, ncol=2)

Legend labels were too wide. The labels have been resized to 0.61, 0.56, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Legend labels were too wide. The labels have been resized to 0.51, 0.53, 0.56, 0.56, 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Some legend labels were too wide. These labels have been resized to 0.54. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Guardar mapas de diciembre

tmap_save(
  tmap_arrange(m0_dic2, m1_dic2, m2_dic2, m3_dic2,nrow=2, ncol=2), 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapas_diciembre_ac_completo.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapas_diciembre_ac_completo.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m0_dic2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_temp_diciembre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_temp_diciembre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m1_dic2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_IMoranloc_diciembre.png", 
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_IMoranloc_diciembre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m2_dic2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_p-values_diciembre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_p-values_diciembre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

tmap_save(
  m3_dic2, 
  filename = "C:/Users/criaguirre/Documents/00_tallegeo/taller01_geost/SalidaTaller3/Mapas/mapa_quad_diciembre.png",
  width = 10, height = 8, 
  dpi = 300                
)

Map saved to C:\Users\criaguirre\Documents\00_tallegeo\taller01_geost\SalidaTaller3\Mapas\mapa_quad_diciembre.png

Resolution: 3000 by 2400 pixels

Size: 10 by 8 inches (300 dpi)

En los mapas de enero, se observan patrones bastante uniformes en las zonas costeras, con altos valores del Índice de Moran Local (Ii). Esto indica una fuerte relación espacial entre las estaciones cercanas, influenciada por la proximidad al océano Pacífico. En cambio, en las zonas montañosas, esta relación disminuye notablemente. Las áreas de mayor altitud muestran una mayor dispersión en los valores de Ii, lo que refleja una variabilidad térmica más intensa, influida por la topografía. Esto se traduce en la presencia de cuadrantes “baja-alta” o “baja-baja”, que resaltan la complejidad térmica debido a las diferencias de altitud.

En junio, las zonas costeras siguen mostrando una relación espacial significativa, aunque menos intensa que en enero. En las montañas, la dispersión de los valores de Ii es aún mayor y la relación espacial prácticamente desaparece en muchas áreas. Los p-valores elevados en estas regiones sugieren que las correlaciones observadas no son significativas. En estas áreas predominan los cuadrantes “baja-baja”, mientras que las zonas costeras mantienen una tendencia de agrupación “alta-alta”.

En octubre, los patrones son similares a los de enero, con una fuerte relación espacial en las zonas costeras, favorecida por la estabilidad térmica que aporta la cercanía al mar. En las montañas, como en los meses anteriores, la correlación espacial es más débil, especialmente en las zonas de mayor altitud, donde los valores de Ii son bajos y los p-valores altos, lo que indica una baja correlación en estas áreas. Las costas siguen mostrando una tendencia “alta-alta”, mientras que las montañas alternan entre los cuadrantes “baja-alta” y “baja-baja”.

En diciembre, los patrones se parecen nuevamente a los de junio, con una correlación positiva notable en las zonas costeras, aunque menos marcada que en los meses de verano. En las montañas, la relación espacial sigue siendo baja o prácticamente nula, lo que refleja la complejidad térmica y la dificultad de capturar con precisión las variaciones en estas áreas. Los p-valores altos refuerzan la idea de que la correlación en estas zonas no es significativa, mientras que los cuadrantes costeros continúan siendo “alta-alta” y en las montañas predominan los cuadrantes “baja-alta” o “baja-baja”.

El análisis del Índice de Moran Local a lo largo de los cuatro meses revela áreas con patrones de temperatura consistentes, que varían según la estación y la geografía de la región. La alta cordillera y los valles centrales presentan los clústeres más marcados, subrayando la influencia de la altitud en la distribución de la temperatura. Además, la significancia estadística de algunos puntos sugiere que la distribución de la temperatura depende en gran medida de la ubicación geográfica y la topografía, lo que destaca la importancia de considerar estos factores al modelar la distribución espacial de la temperatura, especialmente en regiones con tanta variabilidad altitudinal como la Región de Coquimbo.

Sin embargo, a diferencia del análisis de Moran Global, donde los p-valores fueron menores a 0.05, en el análisis local la mayoría de los puntos mostraron p-valores superiores a 0.05. Esto significa que gran parte de los puntos no presentaron una autocorrelación local significativa. Esto podría explicarse porque el p-valor de LISA mide la significancia de la autocorrelación para cada punto específico, es decir, si la temperatura en una estación está correlacionada de forma significativa con las estaciones vecinas. Un p-valor mayor a 0.05 sugiere que no hay un patrón claro de agrupación a nivel local para ese punto en particular.

3.4. Analisis de correlación mediante Variograma

El variograma es una herramienta que muestra cómo varía una característica (como la temperatura o la concentración de minerales) entre diferentes lugares a medida que aumenta la distancia entre ellos. En geostatística, es fundamental para entender cómo se comportan fenómenos naturales en el espacio, como el análisis de suelos o la distribución de depósitos minerales, ayudando a predecir cómo cambian sus propiedades de un punto a otro (Bohling, 2005).

Un variograma experimental muestra cómo varía una característica entre puntos a diferentes distancias, representando gráficamente la semivarianza o diferencia promedio entre estos puntos. A medida que la distancia entre ellos aumenta, la semivarianza también crece, hasta que llega a un punto en el que se estabiliza.

Para entender mejor esta relación, se ajusta un variograma teórico al experimental usando modelos matemáticos como el esférico, exponencial o gaussiano. Esto permite predecir la variación de la característica en áreas donde no hay mediciones. El proceso implica ajustar tres parámetros:

Nugget: Representa la variación a distancias muy cortas, que puede incluir errores de medición. Sill: Es el valor máximo de variación alcanzado. Range: Es la distancia a partir de la cual los puntos dejan de estar correlacionados.

El variograma modelado es esencial para métodos de interpolación como el kriging, que estima valores en lugares no muestreados usando la dependencia espacial de los datos. Los variogramas ajustados serán útiles para el taller 4.

A continuación, se detalla el proceso de la confección de los variogramas para los meses objetivos.

Generación de variogramas experimentales

v_ene <- variogram(temp~1, data_ene_df)
v_jun <- variogram(temp~1, data_jun_df)
v_oct <- variogram(temp~1, data_oct_df)
v_dic <- variogram(temp~1, data_dic_df)

Variograma Experimental Enero

plot(v_ene, main = "Variograma Experimental Enero", plot.numbers = TRUE)

Variograma Empírico Junio

plot(v_jun, main = "Variograma Experimental Junio", plot.numbers = TRUE)

Variograma Empírico Octubre

plot(v_oct, main = "Variograma Experimental Octubre", plot.numbers = TRUE)

Variograma Empírico Diciembre

plot(v_dic, main = "Variograma Experimental Diciembre", plot.numbers = TRUE)

Ajuste de un modelo teórico al variograma empirico

mv_ene <- vgm(psill = 12, 'Sph', range = 80000, nugget = 0)
mv_jun <- vgm(psill = 15, 'Gau', range = 70000, nugget = 1)
mv_oct <- vgm(psill = 20, 'Gau', range = 70000, nugget = 1)
mv_dic <- vgm(psill = 12, 'Gau', range = 60000, nugget = 1)

plot(v_ene,mv_ene, main = "Variograma Teórico Enero", plot.numbers = TRUE)

plot(v_jun,mv_jun,  main = "Variograma Teórico Junio", plot.numbers = TRUE)

plot(v_oct,mv_oct,  main = "Variograma Teórico Octubre", plot.numbers = TRUE)

plot(v_dic,mv_dic,  main = "Variograma Teórico Diciembre", plot.numbers = TRUE)

Ajuste el modelo teórico al Variograma empírico

fitv_ene <- fit.variogram(v_ene, mv_ene)
fitv_jun <- fit.variogram(v_jun, mv_jun)
fitv_oct <- fit.variogram(v_oct, mv_oct)
fitv_dic <- fit.variogram(v_dic, mv_dic)

Variograma ajustado Enero

plot(v_ene, fitv_ene, main = "Variograma Ajustado Enero", plot.numbers = TRUE)

Variograma ajustado Junio

plot(v_jun, fitv_jun,main = "Variograma Ajustado Junio", plot.numbers = TRUE)

Variograma ajustado Octubre

plot(v_oct, fitv_oct, main = "Variograma Ajustado Octubre", plot.numbers = TRUE)

Variograma ajustado Diciembre

plot(v_dic, fitv_dic, main = "Variograma Ajustado Diciembre", plot.numbers = TRUE)

Exportar variogramas ajustado

# Guardar el variograma ajustado fitv_ene en un archivo .rds
saveRDS(fitv_ene, file = "fitv_ene_variogram.rds")
saveRDS(fitv_jun, file = "fitv_jun_variogram.rds")
saveRDS(fitv_oct, file = "fitv_oct_variogram.rds")
saveRDS(fitv_dic, file = "fitv_dic_variogram.rds")

Evaluación de anisotropía

v.dir_ene <- variogram(temp~1, st_as_sf(data_ene_df,coords =c('x','y')),alpha=(0:3)*45)
v.dir_jun <- variogram(temp~1, st_as_sf(data_jun_df,coords =c('x','y')),alpha=(0:3)*45)
v.dir_oct <- variogram(temp~1, st_as_sf(data_oct_df,coords =c('x','y')),alpha=(0:3)*45)
v.dir_dic <- variogram(temp~1, st_as_sf(data_dic_df,coords =c('x','y')),alpha=(0:3)*45)

plot(v.dir_ene, main = "Anisotropía Variograma Enero")

plot(v.dir_jun, main = "Anisotropía Variograma Junio")

plot(v.dir_oct, main = "Anisotropía Variograma octubre")

plot(v.dir_dic, main = "Anisotropía Variograma Diciembre")

Ajuste del modelo experimental al teórico Anisótropo

v.anis_ene <- vgm(psill = 12, model = "Sph", range = 80000, nugget = 0, anis = c(45,0.3))
v.anis_jun <- vgm(psill = 15, model = "Gau", range = 70000, nugget = 1, anis = c(45,0.3))
v.anis_oct <- vgm(psill = 20, model = "Gau", range = 70000, nugget = 1, anis = c(45,0.3))
v.anis_dic <- vgm(psill = 12, model = "Gau", range = 60000, nugget = 1, anis = c(45,0.3))

plot(v.dir_ene, v.anis_ene, main = "Anisotropia Variograma Ajustado Enero")

Junio

plot(v.dir_jun, v.anis_jun, main = "Anisotropia Variograma Ajustado Junio")

Octubre

plot(v.dir_oct, v.anis_oct, main = "Anisotropia Variograma Ajustado Octubre")

Diciembre

plot(v.dir_dic, v.anis_dic, main = "Anisotropia Variograma Ajustado Dicimbre")

Exportar variogramas ajustado

# Guardar el variograma ajustado fitv_ene en un archivo .rds
saveRDS(v.anis_ene, file = "anisv_ene_variogram.rds")
saveRDS(v.anis_jun, file = "anisv_jun_variogram.rds")
saveRDS(v.anis_oct, file = "anisv_oct_variogram.rds")
saveRDS(v.anis_dic, file = "anisv_dic_variogram.rds")

Mapa variograma

Enero

v_ene_map <- variogram(temp~1, st_as_sf(data_ene_df,coords =c('x','y')),map=TRUE, cutoff=100000, width= 10000)
plot(v_ene_map, main = "Mapa de variograma para Enero")

Junio

v_jun_map <- variogram(temp~1, st_as_sf(data_jun_df,coords =c('x','y')),map=TRUE, cutoff=100000, width= 10000)
plot(v_jun_map, main = "Mapa de variograma para Junio")

Octubre

v_oct_map <- variogram(temp~1, st_as_sf(data_oct_df,coords =c('x','y')),map=TRUE, cutoff=100000, width= 10000)
plot(v_oct_map, main = "Mapa de variograma para Octubre")

Diciembre

v_dic_map <- variogram(temp~1, st_as_sf(data_dic_df,coords =c('x','y')),map=TRUE, cutoff=100000, width= 10000)
plot(v_dic_map, main = "Mapa de variograma para Diciembre")

En general, el modelo de variograma teórico que mejor se ajustó fue el gaussiano, excepto en el mes de enero, donde se utilizó el modelo esférico. La alta variabilidad observada en algunos puntos cercanos genera ruido al modelar los variogramas, lo que podría estar relacionado con cambios bruscos en las condiciones geográficas locales, como se ha visto en la distribución de los puntos y su relación con la geografía. Esto sugiere que los variogramas reflejan una estructura de dependencia espacial asociada a un fenómeno natural como la temperatura, mientras que los valores atípicos podrían estar vinculados a variaciones abruptas en las condiciones geográficas. Estos variogramas ajustados serán útiles para futuros análisis y modelado.

En cuanto a la anisotropía, se observa una continuidad de los valores a lo largo de la dirección norte-sur (0°), especialmente en las zonas costeras, valles centrales y la cordillera principal. Por otro lado, de este a oeste (90° y 135°), se aprecia una mayor variabilidad de los datos, estabilizándose a una distancia de entre 60,000 y 70,000 metros, mientras que en la diagonal, la distancia requerida para alcanzar la meseta es mayor, debido a las condiciones mencionadas. En general, se evidencia una anisotropía norte-sur, con una mayor continuidad de valores en esa dirección, lo que también es reflejado a los mapas de los variograma para cada mes.

3.5. Analisis de variograma de los residuos Taller 2

En el taller 2, el modelamiento de regresión lineal o Trend Surface Fitting fue el que presento un menor RMSE por lo que presenta un mejor modelo y un R2 significativamente mayor, por lo que se tiene un mejor ajuste el modelo a los datos.

Se replica las superficies interpoladas en este taller y se obtienen los valores residuales.

# Enero
pred_ene_ras_st <- pred_ene_ras |> st_as_stars() |>  split('band')
temp_lm_ene <- krige(temp ~ dem + slope + dist_costa + aspect + ndvi_ene + lst_ene, 
                     data_ene_df, 
                     pred_ene_ras_st)

[ordinary or weighted least squares prediction]

temp_lm_ene <- temp_lm_ene |> rast() |> mask(vect(reg_coq_utm)) |> trim()

# Junio
pred_jun_ras_st <- pred_jun_ras |> st_as_stars() |>  split('band')
temp_lm_jun <- krige(temp ~ dem + slope + dist_costa + aspect + ndvi_jun + lst_jun, 
                     data_jun_df, 
                     pred_jun_ras_st)

[ordinary or weighted least squares prediction]

temp_lm_jun <- temp_lm_jun |> rast() |> mask(vect(reg_coq_utm)) |> trim()

# Octubre
pred_oct_ras_st <- pred_oct_ras |> st_as_stars() |>  split('band')
temp_lm_oct <- krige(temp ~ dem + slope + dist_costa + aspect + ndvi_oct + lst_oct, 
                     data_oct_df, 
                     pred_oct_ras_st)

[ordinary or weighted least squares prediction]

temp_lm_oct <- temp_lm_oct |> rast() |> mask(vect(reg_coq_utm)) |> trim()

# Diciembre
pred_dic_ras_st <- pred_dic_ras |> st_as_stars() |>  split('band')
temp_lm_dic <- krige(temp ~ dem + slope + dist_costa + aspect + ndvi_dic + lst_dic, 
                     data_dic_df, 
                     pred_dic_ras_st)

[ordinary or weighted least squares prediction]

temp_lm_dic <- temp_lm_dic |> rast() |> mask(vect(reg_coq_utm)) |> trim()

Extraer valores predichos y calculo del residual

# Extraer los valores del raster para las coordenadas de data_ene_df
valores_ras_ene <- terra::extract(temp_lm_ene, vect(data_ene_df))
valores_ras_jun <- terra::extract(temp_lm_jun, vect(data_jun_df))
valores_ras_oct <- terra::extract(temp_lm_oct, vect(data_oct_df))
valores_ras_dic <- terra::extract(temp_lm_dic, vect(data_dic_df))

# Unir los valores extraídos al dataframe original
data_ene_df$pred_lm <- valores_ras_ene[, 2]
data_jun_df$pred_lm <- valores_ras_jun[, 2]
data_oct_df$pred_lm <- valores_ras_oct[, 2]
data_dic_df$pred_lm <- valores_ras_dic[, 2]

data_ene_df$var_lm <- valores_ras_ene[, 3]
data_jun_df$var_lm <- valores_ras_jun[, 3]
data_oct_df$var_lm <- valores_ras_oct[, 3]
data_dic_df$var_lm <- valores_ras_dic[, 3]

# Calculo de residuales
data_ene_df$residual <- data_ene_df$temp - data_ene_df$pred_lm
data_jun_df$residual <- data_jun_df$temp - data_jun_df$pred_lm
data_oct_df$residual <- data_oct_df$temp - data_oct_df$pred_lm
data_dic_df$residual <- data_dic_df$temp - data_dic_df$pred_lm

Calculo de variograma empírico

v_res_ene <- variogram(residual~1, data_ene_df,cutoff = 80000) # se ajusta el cutoff a 80000 para convergere variograma
v_res_jun <- variogram(residual~1, data_jun_df,cutoff = 80000) # se ajusta el cutoff a 80000 para convergere variograma
v_res_oct <- variogram(residual~1, data_oct_df)
v_res_dic <- variogram(residual~1, data_dic_df)

Variograma enero

plot(v_res_ene, main = "Variograma de Residuales Enero", plot.numbers = TRUE)

Variograma junio

plot(v_res_jun, main = "Variograma de Residuales Junio", plot.numbers = TRUE)

Variograma octubre

plot(v_res_oct, main = "Variograma de Residuales Octubre", plot.numbers = TRUE)

Variograma diciembre

plot(v_res_dic, main = "Variograma de Residuales Diciembre", plot.numbers = TRUE)

Variograma teórico

mv_res_ene <- vgm(psill = 2, 'Pen', range = 14000, nugget = 1.5)
mv_res_jun <- vgm(psill = 1.5, 'Gau', range = 40000, nugget = 0.4)
mv_res_oct <- vgm(psill = 2.5, 'Gau', range = 60000, nugget = 0.5)
mv_res_dic <- vgm(psill = 3, 'Gau', range = 60000, nugget = 0.5)

Variograma enero

plot(v_res_ene,mv_res_ene, main = "Variograma de Residuales Enero", plot.numbers = TRUE)

Variograma junio

plot(v_res_jun,mv_res_jun, main = "Variograma de Residuales Junio", plot.numbers = TRUE)

Variograma octubre

plot(v_res_oct,mv_res_oct, main = "Variograma de Residuales Octubre", plot.numbers = TRUE)

Variograma diciembre

plot(v_res_dic,mv_res_dic, main = "Variograma de Residuales Diciembre", plot.numbers = TRUE)

Ajustar variograma

fitv_res_ene <- fit.variogram(v_res_ene, mv_res_ene, fit.method = 7)

Warning in fit.variogram(v_res_ene, mv_res_ene, fit.method = 7): No convergence
after 200 iterations: try different initial values?

fitv_res_jun <- fit.variogram(v_res_jun, mv_res_jun, fit.method = 2)
fitv_res_oct <- fit.variogram(v_res_oct, mv_res_oct)
fitv_res_dic <- fit.variogram(v_res_dic, mv_res_dic)

NOTA: No se pudo ajustar del todo el variograma de los residuos de enero, pese a que se probaron diferentes modelos.

Enero

plot(v_res_ene, fitv_res_ene, main = "Variograma Ajustado de Residuales Enero", plot.numbers = TRUE)

Junio

plot(v_res_jun, fitv_res_jun, main = "Variograma Ajustado de Residuales Junio", plot.numbers = TRUE)

Octubre

plot(v_res_oct, fitv_res_oct, main = "Variograma Ajustado de Residuales Octubre", plot.numbers = TRUE)

Diciembre

plot(v_res_dic, fitv_res_dic, main = "Variograma Ajustado de Residuales Diciembre", plot.numbers = TRUE)

Exportar variograma de los residuales

saveRDS(fitv_res_ene, file = "res_ene_variogram.rds")
saveRDS(fitv_res_jun, file = "res_jun_variogram.rds")
saveRDS(fitv_res_oct, file = "res_oct_variogram.rds")
saveRDS(fitv_res_dic, file = "res_dic_variogram.rds")

4. Discusión y conclusión

El análisis de la autocorrelación espacial, realizado con diferentes herramientas como Moran’s I global y local, así como el uso de variogramas, reveló una relación compleja entre las temperaturas y su distribución espacial a lo largo de los meses. El análisis global de Moran’s I mostró que las temperaturas en la región tienden a agruparse de manera significativa durante todos los meses, lo que indica que valores de temperatura similares se encuentran más cerca unos de otros de lo que sería esperable por azar, según lo sugerido por los bajos p-valores obtenidos. Sin embargo, al estudiar el Índice de Moran local, se observó que la significancia de la autocorrelación varía mucho dependiendo de cada punto de muestreo. Mientras algunos lugares mostraron agrupamientos claros de altas o bajas temperaturas, en otros puntos la autocorrelación no fue significativa, lo que refleja que, aunque hay una tendencia general de agrupamiento, las condiciones locales, como la geografía y la topografía, juegan un papel importante a escala más pequeña.

El uso de variogramas, tanto experimentales como ajustados a modelos teóricos, permitió entender mejor la dependencia espacial de la temperatura. El modelo gaussiano fue el que mejor se ajustó para la mayoría de los meses, lo que sugiere que la variabilidad de la temperatura crece de forma gradual con la distancia, típico de muchos fenómenos naturales. Sin embargo, en enero se utilizó un modelo esférico, lo que indica una transición más brusca de la variabilidad, posiblemente causada por las condiciones en la alta cordillera, donde la variabilidad de las temperaturas es mayor y genera más “ruido” en los datos. Además, se detectó una anisotropía marcada: las temperaturas tienden a ser más continuas de norte a sur, mientras que de este a oeste presentan mayores fluctuaciones. Esto se aprecia especialmente en las zonas costeras, los valles centrales y la cordillera, donde la geografía local, como la altitud y la cercanía al mar, influye en la estabilidad térmica.

La relación entre la geografía y la autocorrelación espacial también se evidencia en la variabilidad observada en ciertos puntos específicos que no presentan una autocorrelación significativa. Estos puntos parecen estar influenciados por condiciones locales como la orientación del terreno, la proximidad a cuerpos de agua o la altitud. Este comportamiento coincide con lo observado en los variogramas, donde la variabilidad en las mediciones puede complicar el ajuste de los modelos. Los valores atípicos detectados también subrayan la importancia de considerar estos factores geográficos y climáticos al analizar la distribución de la temperatura, ya que pueden introducir cambios significativos en los resultados.

En conclusión, el análisis de la autocorrelación espacial y los variogramas ha permitido identificar patrones generales de agrupamiento de temperaturas en la región, al mismo tiempo que ha resaltado la influencia de la variabilidad local y de las características geográficas en estos patrones. Mientras que el Moran’s I global proporciona una visión amplia de la tendencia al agrupamiento, el análisis local y los variogramas ofrecen detalles sobre las diferencias puntuales, como la influencia de la distancia y las condiciones geográficas. En otro contexto, estos valores atípicos podrían haberse descartado si no se comprendiera su relación con factores geográficos y estacionales. Una mayor cantidad de puntos de observación, con una distribución más uniforme, podría mejorar la continuidad de los datos y captar de manera más precisa la variabilidad de la temperatura.

5. Referencias

Bohling, G. (2005). Introduction to geostatistics and variogram analysis. Kansas geological survey, 1(10), 1-20. Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93-115. Recuperado de https://dces.webhosting.cals.wisc.edu/wp-content/uploads/sites/128/2013/08/W4_Anselin1995.pdf ESRI Support. (s.f.). GIS Dictionary: Variogram. Disponible en https://support.esri.com/es-es/gis-dictionary/variogram ArcGIS Pro Help. (s.f.). Modeling a Semivariogram. Disponible en https://pro.arcgis.com/es/pro-app/latest/help/analysis/geostatistical-analyst/modeling-a-semivariogram.htm