Análisis espacio temporal de la temperatura y precipitación en la península de Yucatán
Parcial 1 Teoría de Estadística Espacial
Camila Alejandra Pinzón Caro
2025-09-12
Planteamiento del problema
La península de Yucatán se extiende entre el Golfo de México al oeste y norte, y el Mar Caribe al este. Los sitios arqueológicos que en esta se encuentran, como Chichén Itza y Uxmal (designados como Patrimonio de la Humanidad por la UNESCO) sumado a las playas de Cancún hacen de esta península uno de los puntos más turísticos del país. Esto resalta la importancia de analizar sus condiciones climáticas, pues es una variable a tener en cuenta al momento de planear una visita al territorio.
El clima de Yucatán es uno de los más calurosos de México; esto se debe a la ubicación geográfica. Las altas temperaturas azotan el estado durante todo el año, dando a la zona un clima semihúmedo y caluroso. Aproximadamente el 85% del territorio recibe entre 24° C y 30° C durante todo el año. El estado de Yucatán puede llegar a los 42° C en verano.
Descripción de los datos
Los datos usados provienen de la información histórica de las estaciones climatológicas convencionales que conforman la Red Nacional de la Comisión Nacional de Agua (CONAGUA). Se tomó una muestra de 33 estaciones ubicadas al sureste de México en la península de Yucatán, concretamente en los estados de Yucatán, Quintana Roo y Campeche. Esta zona es una de las más famosas del país por sus playas en el golfo de México y las ruinas mayas, resaltando la importancia de analizar las condiciones climáticas de esta.
Se seleccionaron las variables Temperatura máxima diaria y Precipitación acumulada diaria desde 2002-01-01 hasta 2024-12-31, resultando en 8401 datos. Aunque se cuenta con datos históricos anteriores, se toma este periodo de tiempo para no tener tantos datos faltantes en la base de datos.
Los datos se pueden consultar en: https://smn.conagua.gob.mx/es/climatologia/informacion-climatologica/informacion-estadistica-climatologica
Las unidades de las variables son:
- Temperatura máxima diaria: ºC
- Precipitación acumulada diaria: mm
Para este primer análisis solamente espacial, se toman los datos recolectados del día 2018-06-16.
Carga de librerías
library(readxl)
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-6 (built on 2025-08-29) is now loaded
## --------------------------------------------------------------library(sp)
library(sf)## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(tmap)
library(rnaturalearth)
library(ggplot2)
setwd("/Users/macbook/Desktop/Universidad/Maestria Estadistica/Estadistica Espacial")Lecura de datos
base_PRECIP <- read_excel("base_PRECIP_2.xlsx", sheet = "Sheet1")
base_TMAX <- read_excel("base_TMAX_2.xlsx", sheet = "Sheet1")
ubicaciones <- read_excel("base_TMAX_2.xlsx", sheet = "Ubicaciones")## New names:
## • `` -> `...1`fecha <- as.POSIXct("2018-06-16", tz="UTC")
dato_TMAX <- base_TMAX[base_TMAX$FECHA==fecha,]
dato_PRECIP <- base_PRECIP[base_PRECIP$FECHA==fecha,]
df_clima <- cbind(Temp = as.numeric(dato_TMAX[-1]), Precip = as.numeric(dato_PRECIP[-1]))
df_clima <- as.data.frame(df_clima)
head(df_clima)## Temp Precip
## 1 30 9.0
## 2 32 28.0
## 3 34 29.0
## 4 32 36.2
## 5 30 40.0
## 6 30 40.0par(mfrow=c(1, 2))
boxplot(df_clima$Temp, main = "Boxplot Temperatura máxima")
boxplot(df_clima$Precip, main = "Boxplot Precipitación acumulada")
Convertir a coordenadas planas
# Crear un objeto SpatialPoints
ori_coords <- SpatialPoints(cbind(ubicaciones$LONGITUD, ubicaciones$LATITUD),
proj4string = CRS("+proj=longlat"))
# Convertirlo a UTM
utm_coords <- spTransform(ori_coords, CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
plot(utm_coords, axes = TRUE, main = "Coordenadas UTM")
utm_coords_df <- as.data.frame(utm_coords)
colnames(utm_coords_df) <- c("Easting", "Northing")Mapa de la zona de interés
# Convertir tus puntos de sp → sf
utm_coords_sf <- st_as_sf(utm_coords)
# 2. Descargar el polígono de México
mexico <- ne_countries(scale = "medium", country = "Mexico", returnclass = "sf")
# 3. Transformar México a UTM zona 16N
mexico_utm <- st_transform(mexico, 32616)
# 4. Graficar
ggplot() +
geom_sf(data = mexico_utm, fill = "gray90", color = "black") +
geom_sf(data = utm_coords_sf, color = "red", size = 2) +
theme_minimal()
ggplot() +
geom_sf(data = mexico_utm, fill = "gray90", color = "black") +
geom_sf(data = utm_coords_sf, color = "red", size = 2) +
coord_sf(xlim = range(st_coordinates(utm_coords_sf)[,1]) + c(-60000, 15000),
ylim = range(st_coordinates(utm_coords_sf)[,2]) + c(-30000, 30000)) +
theme_minimal()
Análisis de tendencia
df_clima <- cbind(utm_coords_df, df_clima)
head(df_clima)## Easting Northing Temp Precip
## 1 212518.7 2214240 30 9.0
## 2 173715.3 2213196 32 28.0
## 3 162137.3 2210460 34 29.0
## 4 168065.1 2229804 32 36.2
## 5 172492.9 2234861 30 40.0
## 6 173684.1 2236378 30 40.0pairs(df_clima[,-4])
pairs(df_clima[,-3])
Los graficos de dispersión dejan ver que no se cuenta con ningún tipo de tendencia en los datos ni de Temperatura ni de Precipitación. Es por esta razón que se puede asumir que las medias de los procesos no dependen de la ubicación espacial, luego son de media constante.
# Convertimos los datos a un objeto de la clase geodata
temp.g <- as.geodata(df_clima, data.col = 3)
precip.g <- as.geodata(df_clima, data.col = 4)
# Resumen para la variable Temperatura
summary(temp.g)## Number of data points: 33
##
## Coordinates summary
## Easting Northing
## min 162137.3 2210460
## max 296321.7 2285992
##
## Distance summary
## min max
## 1067.368 145794.046
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 26.00000 28.50000 30.00000 30.93636 33.00000 38.00000plot(temp.g)
# Resumen para la variable Precipitación
summary(precip.g)## Number of data points: 33
##
## Coordinates summary
## Easting Northing
## min 162137.3 2210460
## max 296321.7 2285992
##
## Distance summary
## min max
## 1067.368 145794.046
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 10.00000 25.00000 24.74545 36.20000 55.00000plot(precip.g)
Una vez finalizado este análisis de tendencia, concluimos que dado que la media del proceso espacial es constante y que podemos asumir que la función de covarianza \(C(\mathbf{h})\) existe y es función única del vector de separación \(\mathbf{h}\) tenemos un proceso estacionario débil.
Explorar anisotropía
Un supuesto importante a revisar es el de isotropía. Este es, corroborar si la función de covarianza o semivarianza dependen únicamete de la magnitud de la separación \(\mathbf{h}\). Para esto se llevan a cabo estimaciones empíricas de cuatro semivariogramas direccionales, como se observa en las siguientes gráficas.
var_direc_temp <- variog4(as.geodata(df_clima, data.col = 3))## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramplot(var_direc_temp, main = "Semivariogramas direccionales para la Temperatura máxima diaria")
var_direc_precip <- variog4(as.geodata(df_clima, data.col = 4))## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramplot(var_direc_precip, main = "Semivariogramas direccionales para la Precipitación acumulada diaria")
Note que los semivariogramas se ven diferentes depeniendo de la dirección en la que se consideren, en este sentido los semivariogramas serían anisotrópicos.
Estimación del semivariograma
El semivariograma se define como la varianza de la variable incrementos. Este se estima empíricamente como \[\hat{\gamma}(\mathbf{h})=\frac{1}{2\mid N(\mathbf{h})\mid}\sum_{N(\mathbf{h})}[Z(\mathbf{s+h})-Z(\mathbf{s})]^2\]
vari_temp <- variog(as.geodata(df_clima, data.col = 3), max.dist = 100000)## variog: computing omnidirectional variogramvari_precip <- variog(as.geodata(df_clima, data.col = 4), max.dist = 100000)## variog: computing omnidirectional variogram#plot(vari_temp)
plot(vari_precip, main = "Semivariograma empírico para la Precipitación acumulada diaria")
rob_vari_temp <- variog(as.geodata(df_clima, data.col = 3),
estimator.type = "modulus", max.dist = 100000)## variog: computing omnidirectional variogramrob_vari_precip <- variog(as.geodata(df_clima, data.col = 4),
estimator.type = "modulus", max.dist = 100000)## variog: computing omnidirectional variogram#plot(rob_vari_temp)
plot(rob_vari_precip, main = "Semivariograma empírico robusto para la Precipitación acumulada diaria")
La función de semivarianza empírica robusta muestra un comportamiento más sinusoidal que la estimada por el método clásico, además de un menor efecto pepita. En ambos casos no se observa que la función de semivarianza tienda concretamente hacia un punto de silla \(C(\mathbf{0})=\sigma^2\), aunque si se evidencia su comportamiento creciente.