df1_sf <- st_as_sf(df1, coords = c("LONGITUD", "LATITUD"), crs = 4326) # WGS84
df1_utm <- st_transform(df1_sf, crs = 32630) # UTM Zona 30N para MadridProyecto Geoestadística. 2025-II
Estadística Espacial
Objetivo General
Estimar y analizar la distribución espacial de la temperatura y la humedad en la ciudad de Madrid mediante técnicas de geoestadística, con el fin de modelar su comportamiento espacial y generar predicciones en ubicaciones no muestreadas.
Objetivos Especificos
Analizar la relación entre la distribución espacial de la temperatura y la humedad con la presencia de ríos o canales en la ciudad de Madrid, para identificar posibles patrones de influencia hidrológica sobre estas variables climáticas.
Evaluar la influencia de la cobertura vegetal, en particular árboles y zonas verdes, sobre la distribución espacial de la temperatura y la humedad en la ciudad de Madrid.
Modelar la estructura espacial de la temperatura y la humedad mediante semivariogramas y ajustar modelos teóricos que permitan comprender la dependencia espacial de estas variables en la ciudad de Madrid.
Descripción de los datos
Los Datos utilizados en este proyecto fueron obtenidos del Portal de datos abiertos del Ayuntamiento de Madrid. Se emplearon registros correspondientes al año 2023, recopilados por un total de 26 estaciones meteorológicas distribuidas en la ciudad. Los datos utilizados en este proyecto fueron obtenidos del Portal de datos abiertos del Ayuntamiento de Madrid Se emplearon registros correspondientes al año 2023, recopilados por un total de 26 estaciones meteorológicas distribuidas en la ciudad. Las variables disponibles en el conjunto de datos son las siguientes:
Consideraciones de las mediciones
No todas las 26 estaciones meteorológicas registran las mismas variables, y en algunos casos, incluso aquellas que sí lo hacen pueden no contar con mediciones disponibles en determinados días. Junto a cada registro diario de una variable, se incluye una columna de validación que indica si la medición fue efectivamente tomada por la estación correspondiente: “V” si fue registrada ese día y “F” si no lo fue.
Existen estaciones meteorológicas que, aunque presentan valores diarios distintos validados con ‘V’, comparten exactamente la misma ubicación geográfica. Para este trabajo, se optó por realizar una ligera modificación en la ubicación de una de las estaciones que se encontraban demasiado próximas (J.D.M Vallecas 2), con el fin de evitar conflictos en el análisis espacial.
Instrumentos de medición
Higrómetro: Instrumento de medición para la Humedad relativa:
Sonda PT-100: Instrumento de medición de la Temperatura.
Temperatura
La primera variable analizada en este estudio es la temperatura, correspondiente al día 1 de enero de 2023. Cabe destacar que no todas las estaciones meteorológicas de Madrid registran información sobre esta variable; en total, 24 estaciones reportaron datos de temperatura para esa fecha.
- La variable temperatura en la ciudad de Madrid es medida en grados centígrados
1. Transformación de coordenadas UTM (Coordenadas planas)
CRS (Coordinate Reference System) significa el sistema de coordenadas inicial que en este caso es WGS84 el más común.
st_transform() : Es la función de transformación de coordenadas geograficas a coordenadas planas.
CRS = 32630 es la zona UTM para adecuada para la ciudad de Madrid.
2. Gráficos de interacción
Gráficos de Interacción: Coordenadas Geográficas
Gráficos de Interacción: Coordenadas planas
Los patrones espaciales no cambian al transformar las coordenadas, pero se recomienda usar sistemas de coordenadas planas (como UTM) para mayor precisión en análisis locales.
3. Matriz de Correlación
| Este | Norte | Temperatura | |
|---|---|---|---|
| Este | 1.00 | -0.26 | 0.19 |
| Norte | -0.26 | 1.00 | -0.09 |
| Temperatura | 0.19 | -0.09 | 1.00 |
Parece no haber un patron claro de ajuste a los datos de temperatura.
4. Gráfico Interactivo de Tendencia con Ploty
5. Modelo de Regresión
El mejor modelo encontrado fue el siguiente: \[ \hat{T} = \beta_0 + \beta_1 \cdot \text{Este} + \beta_2 \cdot \text{Norte} + \beta_3 \cdot (\text{Este} \cdot \text{Norte}) + \beta_4 \cdot \text{Este}^2 + \beta_5 \cdot \text{Norte}^2 + \varepsilon \]
Call:
lm(formula = D01 ~ Este + Norte + Este:Norte + I(Este^2) + I(Norte^2),
data = df1)
Residuals:
Min 1Q Median 3Q Max
-3.9779 -0.3249 0.4127 0.6572 1.6037
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.600e+05 3.534e+05 -1.302 0.2095
Este 1.951e-01 1.101e-01 1.772 0.0933 .
Norte 1.863e-01 1.540e-01 1.210 0.2420
I(Este^2) -5.718e-08 2.847e-08 -2.009 0.0598 .
I(Norte^2) -1.922e-08 1.674e-08 -1.148 0.2661
Este:Norte -3.228e-08 2.209e-08 -1.461 0.1612
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.45 on 18 degrees of freedom
Multiple R-squared: 0.3223, Adjusted R-squared: 0.1341
F-statistic: 1.712 on 5 and 18 DF, p-value: 0.1828| Métrica | Valor |
|---|---|
| R² | 0.322 |
| R² Ajustado | 0.134 |
Un resultado interesante es que el término \(Este^2\) resultó significativo incluso con un nivel de significancia relativamente alto (α>0,5) Sin embargo, tanto el \(R^2\) como el \(R^2\) ajustado fueron bajos, lo que sugiere que no se logró detectar una tendencia clara, ni lineal ni polinómica.
| Grados de Libertad | Suma de Cuadrados | Media Cuadrática | Estadístico F | Valor P | |
|---|---|---|---|---|---|
| Este | 1 | 1.939 | 1.939 | 0.922 | 0.350 |
| Norte | 1 | 0.080 | 0.080 | 0.038 | 0.847 |
| I(Este^2) | 1 | 10.683 | 10.683 | 5.080 | 0.037 |
| I(Norte^2) | 1 | 0.814 | 0.814 | 0.387 | 0.542 |
| Este:Norte | 1 | 4.490 | 4.490 | 2.135 | 0.161 |
| Residuals | 18 | 37.854 | 2.103 | NA | NA |
Justificación de la exclusión de la tendencia espacial
El modelo anterior se utilizó únicamente con fines exploratorios para evidenciar la baja presencia de tendencia espacial en los datos. Sin embargo, no se considera un modelo parsimonioso, y varios de sus coeficientes no resultaron estadísticamente significativos. Por lo tanto, se optó por no incorporar una estructura de tendencia en el análisis final para Temperatura .
6. Mapa interactivo con mapview
Para esta visualización, se identificó que dos estaciones meteorológicas compartían exactamente las mismas coordenadas espaciales con otras dos estaciones, a pesar de registrar valores distintos de temperatura. Para evitar problemas en el análisis geoestadístico, se optó por modificar ligeramente sus ubicaciones, desplazándolas a puntos cercanos a su posición original.
Se observa que las zonas periféricas de Madrid, especialmente aquellas cercanas a áreas rurales o espacios verdes, presentan los valores más bajos de temperatura. En contraste, las zonas urbanas muestran temperaturas superiores en al menos 4 grados Celsius con respecto a estas áreas, lo que evidencia un marcado efecto de isla de calor urbana.
7. Variograma
En el semivariograma de temperatura no se usará residuales, los gráficos que se muestran acontinuación no evidencian un cambio entre sin agregar residuos o con agregar residuos.
variog: computing omnidirectional variogram
Gráfico del variograma empírico
Para este gráfico sólo se usó los puntos del semivariograma que tenian 8 parejas o más.
8. ¿Anisotropía?
A continuación se muestra graficos del semivariograma a diferentes angulos para observar patrones en diferentes direcciones
En el los gráficos de 0° y de 195° no podemos identificar patrones, es similar a un ruido blanco.
En los gráficos de 45° y 90° se puede identificar patrones crecientes.
Conclusión hay anisotropía
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)
9. Estimación de la semivarianza
Como podermos ver en el siguiente gráfico los mejores candidatos respecto a nuestro semivariograma son el modelo wave y circular
Estimación usando Wave
Dado el gráfico del semivariograma, y observando que la semivarianza parece decaer hasta cero en el origen, no incluiremos un efecto pepita en nuestro modelo.
| Modelo | Sigma_Sq | Phi | Tausq | PracticalRange |
|---|---|---|---|---|
| Wave | 2.7 | 1520 | 0 | 4481.232 |
Estimación usando modelo Circular
Dado el gráfico del semivariograma, y observando que la semivarianza parece decaer hasta cero en el origen, no incluiremos un efecto pepita en nuestro modelo Circular
| Modelo | Sigma_Sq | Phi | Tausq | PracticalRange |
|---|---|---|---|---|
| Circular | 3.6 | 8613.57 | 0 | 8613.57 |
variofit: covariance model used is wave
variofit: weights used: cressie
variofit: minimisation function used: optim variofit: covariance model used is circular
variofit: weights used: cressie
variofit: minimisation function used: optim 
Humedad Relativa
1. Descripción:
Se seleccionan únicamente los registros que corresponden a:
La magnitud 86 (humedad relativa)
Año 2023
Mes de Octubre
Día 1 con datos válidos (
V01 == "V")La Humedad relativa es medida de 0 a 100
2. Gráficos de interacción
Gráficos de Interacción: Coordenadas Geográficas
Gráficos de Interacción: Coordenadas planas
Los patrones espaciales no cambian al transformar las coordenadas, pero se recomienda usar sistemas de coordenadas planas (como UTM) para mayor precisión en análisis locales.
3. Matriz de correlación
| Este | Norte | Humedad | |
|---|---|---|---|
| Este | 1.00 | -0.31 | -0.51 |
| Norte | -0.31 | 1.00 | 0.34 |
| Humedad | -0.51 | 0.34 | 1.00 |
4. Gráfico 3D usando scatterplot
Descripción:
Visualiza la humedad relativa en un espacio tridimensional para identificar tendencias espaciales complejas.
5. Modelo de regresión
Descripción:
El modelo de regresión para la captura de tendencia en los datos para la humedad tiene un gran problema, es altamente influenciado por el rio Manzanares que está en la zona Oeste en dirección de norte a sur de Madrid.
Para captar esta información creó una variable dummy para para las estaciones que estuvieran a 200 metros o menos. estas estaciones fueron las siguientes:
\[ \hat{D01} = \beta_0 + \beta_1 \cdot \text{Este} + \beta_2 \cdot \text{Norte} + \beta_3 \cdot (\text{Este} \cdot \text{Norte}) + \beta_4 \cdot \text{Este}^2 + \beta_5 (zona\ ribiera) +\varepsilon \]
Call:
lm(formula = D01 ~ Este + Norte + Este:Norte + I(Este^2) + zona_ribera,
data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.3800 -2.6208 0.1328 1.8879 7.5739
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.113e+05 1.384e+05 3.694 0.00197 **
Este -1.334e+00 3.383e-01 -3.942 0.00117 **
Norte -9.692e-02 2.876e-02 -3.370 0.00390 **
I(Este^2) 3.947e-07 8.562e-08 4.610 0.00029 ***
zona_riberaurbana -6.946e+00 2.869e+00 -2.421 0.02772 *
Este:Norte 2.198e-07 6.520e-08 3.371 0.00389 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.403 on 16 degrees of freedom
Multiple R-squared: 0.7586, Adjusted R-squared: 0.6831
F-statistic: 10.05 on 5 and 16 DF, p-value: 0.0001699| Métrica | Valor |
|---|---|
| R² | 0.7586 |
| R² Ajustado | 0.6831 |
El modelo presentó un \(R^2\) de 0.7586 y un \(R^2\) ajustado de 0.6831, lo que indica que aproximadamente el 68% de la variabilidad en la variable Humedad puede explicarse por las variables incluidas en el modelo. Aunque el ajuste no es perfecto, estos valores sugieren la presencia de una tendencia espacial significativa, principalmente explicada por el término cuadrático de \(Este^2\)
Se destaca que el término cuadrático \(Este^2\) fue altamente significativas \(0.001\)), lo que evidencia una curvatura en la estructura espacial a lo largo del eje Este. Además, tanto el efecto lineal de \(Este\) como la interacción \(Este:Norte\) resultaron significativos, lo que indica una posible dependencia espacial combinada entre ambas coordenadas.
Finalmente, la variable dummy también fue significativa, lo cual refuerza la existencia de un cambio sistemático asociado a las estaciones cercanas al río.
| Suma de Cuadrados | Grados de Libertad | Media Cuadrática | Estadístico F | Valor P | |
|---|---|---|---|---|---|
| Este | 1 | 331.072 | 331.072 | 17.078 | 0.001 |
| Norte | 1 | 44.654 | 44.654 | 2.303 | 0.149 |
| I(Este^2) | 1 | 307.400 | 307.400 | 15.857 | 0.001 |
| zona_ribera | 1 | 71.150 | 71.150 | 3.670 | 0.073 |
| Este:Norte | 1 | 220.323 | 220.323 | 11.365 | 0.004 |
| Residuals | 16 | 310.175 | 19.386 | NA | NA |
Resumen tabla ANOVA
Todos los terminos fueron significativos a un nivel de significancia del 5%
Gracias al valor relativamente alto del \(R^2\) ajustado y a la significancia estadística de todos los coeficientes del modelo, se puede concluir que este logró capturar adecuadamente la estructura espacial presente en los datos.
Mapa interactivo con mapview
Descripción:
Visualiza espacialmente los valores de humedad directamente sobre un mapa con colores.
Podemos ver que al igual con la temperatura nuestros valores mas atípicos son los mas cercanos al rio adyacente a la ciudad (Rio Manzanares).
6. Análisis geoestadístico de los residuos usando con geoR
Descripción:
Se crea un objeto geoestadístico para:
- Ver la distribución espacial de los datos
- Evaluar los residuos del modelo
- Estimar el semivariograma, que mide la dependencia espacial de los datos
Observaciones
- La transformación aplicada ha removido efectivamente las tendencias sistemáticas presentes en los datos, particularmente en la coordenada, Este (X), donde se observaba la mayor influencia.
Este Norte D01 V01 ESTACIÓN LONGITUD LATITUD
1 446059.8 4472326 40 V J.M.D. Moratalaz -3.635637 40.39979
2 439745.8 4466916 51 V J.M.D. Villaverde -3.709525 40.35063
3 437242.9 4477104 67 V Centro Mpal. De Acústica -3.740000 40.44222
4 444327.1 4479330 46 V J.M.D. Hortaleza -3.656667 40.46278
5 439149.8 4480909 48 V Peñagrande -3.717881 40.47663
6 440884.8 4475930 37 V J.M.D.Chamberí -3.696950 40.43191
MUNICIPIO ESTACION MAGNITUD ANO MES zona_ribera residuales
1 79 102 86 2023 10 urbana -2.7279204
2 79 103 86 2023 10 urbana 1.9694771
3 79 106 86 2023 10 ribera 5.6491770
4 79 107 86 2023 10 urbana -0.1714482
5 79 108 86 2023 10 urbana 1.6430503
6 79 109 86 2023 10 urbana -7.3800141
7. Variograma
variog: computing omnidirectional variogram
| Lag | Distancia Media | Semivarianza | Frecuencia |
|---|---|---|---|
| 1 | 815.094 | 15.876 | 5 |
| 2 | 2445.283 | 19.035 | 27 |
| 3 | 4075.472 | 24.823 | 44 |
| 4 | 5705.661 | 29.409 | 33 |
| 5 | 7335.850 | 17.842 | 44 |
| 6 | 8966.039 | 13.235 | 31 |
| 7 | 10596.227 | 20.299 | 20 |
| 8 | 12226.416 | 8.609 | 7 |
| 9 | 13856.605 | 13.395 | 10 |
| 10 | 15486.794 | 9.722 | 4 |
| 11 | 17116.983 | 7.572 | 4 |
¿Anisotropía?
A continuación se muestra graficos del semivariograma a diferentes angulos para observar patrones en diferentes direcciones
- En los todos los gráficos se puede observar una ligera tendencia decendente de manera muy parecida en todos los angulos.
Por ende no parece haber Anisotropía sin embargo aun hay una tendencia en los datos no capturada por la regresión
par(mfrow=c(2,2))
vari_0 <- variog(df1G,
trend = "2nd",
max.dist = 15000,
dir = 0)variog: computing variogram for direction = 0 degrees (0 radians)
tolerance angle = 22.5 degrees (0.393 radians)vari_45 <- variog(df1G,
trend = "2nd",
max.dist = 15000,
dir = pi / 4)variog: computing variogram for direction = 45 degrees (0.785 radians)
tolerance angle = 22.5 degrees (0.393 radians)vari_90 <- variog(df1G,
trend = "2nd",
max.dist = 15000,
dir = pi / 2)variog: computing variogram for direction = 90 degrees (1.571 radians)
tolerance angle = 22.5 degrees (0.393 radians)vari_135 <- variog(df1G,
trend = "2nd",
max.dist = 15000,
dir = 3 * pi / 4)variog: computing variogram for direction = 135 degrees (2.356 radians)
tolerance angle = 22.5 degrees (0.393 radians)par(mfrow = c(2, 2),
mar = c(3, 3, 1, 1),
mgp = c(2, 1, 0))
plot(vari_0, main = "vari 0")
plot(vari_45, main = "vari 45")
plot(vari_90, main = "vari 90")
plot(vari_135, main = "vari 135")
8. Estimación de la semivarianza
En esta sección comparamos distintos modelos teóricos para ajustar la semivarianza espacial observada. Esto nos permitirá seleccionar un modelo adecuado para realizar kriging o interpolación geoestadística.
Como podemos ver en el siguiente gráfico, los modelos wave, circular, y gaussian parecen ajustarse mejor a la forma del semivariograma empírico.
variofit: covariance model used is gaussian
variofit: weights used: cressie
variofit: minimisation function used: optim variofit: covariance model used is wave
variofit: weights used: cressie
variofit: minimisation function used: optim variofit: covariance model used is circular
variofit: weights used: cressie
variofit: minimisation function used: optim 
En el caso de la humedad, resulta más difícil determinar con claridad la presencia de un efecto pepita. No obstante, al ajustar los modelos teóricos sin considerar dicho efecto, se obtiene un buen ajuste a la semivarianza empírica. Por lo tanto, se optará por no incluir el efecto pepita en el modelo final.
Discusiones
Para este análisis, fue relevante que la función
variogde geoR no permita especificar explícitamente una tendencia, ya que persiste la duda sobre la existencia o no de isotropía en la humedad relativa. Dada esta limitación, podría ser conveniente considerar el uso del paquete gstat en futuros análisis, ya que ofrece mayor flexibilidad para modelar tendencias y explorar estructuras espaciales más complejas.Resulta llamativo que, si bien existen estaciones meteorológicas ubicadas muy próximas entre sí, también se identifican amplias zonas urbanas de Madrid que carecen de estaciones en varios kilómetros a la redonda. Esta distribución desigual puede afectar la representatividad espacial de los datos y limita la precisión de los análisis geoestadísticos en ciertas áreas de la ciudad especialmente cerca a el municipio de Coslada.
Los valores más bajos de temperatura registrados el 1 de enero de 2023, así como los valores más altos de humedad relativa observados el 1 de octubre del mismo año, se localizaron en las cercanías de la estación situada junto al río Manzanares (Estacion Meteorológica-Centro Mpal. De Acústica). Este patrón sugiere una posible influencia del río sobre estas variables climáticas.
Conclusiones
Para la temperatura y para la humedad, la semivarianza muestra un comportamiento que sugiere convergencia de la varianza. Esto se debe a que los modelos teóricos que mejor se ajustan a la semivarianza empírica son los modelos wave, los cuales se caracterizan por oscilar alrededor de un valor constante.
Para la temperatura en la ciudad de Madrid no hubo captura de tendencia apartir de métodos de regresión, además resultó anisotropica, lo que significa que no solo basta con conocer cuanto están separados un par de puntos, sino tambien se necesita conocer la orientación de dicha distancia, esto se conoce como un campo aleatorio anisotrópico
La humedad relativa parece comportarse de manera isotrópica, ya que los patrones observados es debido a la ausencia de tendencia capturada en la regresión, por lo tanto que no es necesario considerar la orientación de la distancia. Esto indica que estamos frente a un campo aleatorio aparentemente isotrópico
Se observa una relación significativa entre la humedad y la coordenada X, lo que sugiere la presencia de una tendencia espacial en dirección Este-Oeste. Esta variación podría atribuirse a factores geográficos como la proximidad del río Manzanares y la presencia de la zona montañosa, ambos ubicados en el oeste de la ciudad de Madrid.
No hay efecto pepita en ninguna de nuestras estimaciones del semivariograma.
Kriging Temperatura.
Estimación del Semivariograma
Modelo Teórico propuesto:
Modelo = “Wave”
Sill = 2.5
range = 4300
Aunque este modelo fue planteado como ajuste inicial, no logró converger durante el proceso iterativo de ajuste con gstat. Por esta razón, se utilizó como modelo de partida para la estimación mediante mínimos cuadrados generalizados (GLS), permitiendo así obtener los parámetros ajustados de manera más estable.
y el ajuste por mínimos cuadrados ponderados quedo de la siguiente manera:
Modelo = “Wave”
Sill = 0.6711084
Range = 4300
Kriging temperatura 1 de enero de 2023
Se observan varios aspectos importantes: los valores más altos de temperatura representados por los colores amarillo y verde se concentran en la zona más céntrica de Madrid, que corresponde a las áreas más urbanizadas de la ciudad. En contraste, los colores azul y morado indican temperaturas más bajas, predominando en las zonas periféricas y verdes del municipio, lo que resalta la diferencia térmica entre el Madrid urbano y el Madrid rural. Además, se aprecia que las áreas cercanas al río Manzanares también presentan temperaturas relativamente más bajas seguramente por su correlación con la humedad presente cerca al rio.
[using ordinary kriging]Como se mencionó anteriormente, las estaciones de monitoreo están mayoritariamente concentradas en las zonas urbanas de Madrid. En contraste, las áreas periféricas del municipio presentan una menor densidad de estaciones, lo que se traduce en una mayor varianza y error estándar en las predicciones para esas zonas. Esta falta de cobertura se evidencia claramente en el siguiente mapa.
Validación del modelo.
Las métricas de validación cruzada para el modelo de kriging univariado muestran un comportamiento aceptable, teniendo en cuenta que el modelo solo utiliza la variable temperatura Temperatura como insumo, sin incorporar variables auxiliares. En particular, la correlación de 0.41 entre los valores observados y los predichos representa una correlación mediana, lo que indica que el modelo logra capturar parcialmente la estructura espacial de la variable.
Este valor es coherente con las limitaciones propias de un modelo univariado, especialmente en zonas donde la densidad de estaciones es baja o hay presencia de factores ambientales no modelados directamente, como la influencia del río Manzanares, cuyas condiciones de humedad local pueden estar afectando los patrones térmicos y generando mayor variabilidad local que el modelo no logra explicar completamente.
Aun así, el error medio absoluto (MAE = 1.03) y el error cuadrático medio (RMSE = 1.39) se encuentran dentro de un rango razonable para datos ambientales de este tipo, lo que refuerza la utilidad del modelo como una primera aproximación para entender el comportamiento espacial de la temperatura en el municipio.
| ME | MAE | RMSE | MSPE | COR |
|---|---|---|---|---|
| -0.036 | 1.029 | 1.395 | 1.947 | 0.411 |