Resumen

El clima de Australia varía ampliamente, aunque la mayor parte de Australia es desértico o semiárido, el 40% del territorio está cubierto por médanos. Solo las esquinas sudestes y sudoestes tienen un clima templado y de suelos moderadamente fértiles. La parte norte del país tiene un clima tropical: parte es bosque lluvioso tropical, parte pastizales, y parte desiertos.

El gran objetivo de este trabajo, es analizar la variable climatológica de temperatura en un enfoque geoestadístico, centrándonos en el perímetro del territorio australiano, el cual es bastante diverso en éste ámbito de estudio.

  1. Introducción

Las temperaturas más calientes australianas pueden alcanzar los 50 °C, mientras las mínimas pueden bajar de cero. El continente no tiene grandes variaciones en altitud y está rodeado mayormente por océanos templados. Así resulta que son extremadamente infrecuentes las bajas temperaturas. Consecuentemente, Australia es un país tan grande que goza de muchas condiciones climáticas: desde el calor del desierto hasta las bajas temperaturas nevadas de las montañas; desde las abundantes lluvias hasta la ausencia total de precipitaciones.

Por estas razones, es de gran interés realizar un análisis geoestadístico de la temperatura, con el fin de pensar que tan correlacionadas se pueden encontrar los valores de temperatura medidos y si es posible ajustar un modelo con el fin de llegar a posibles predicciones para la variable de interés en sitios no muestreados.

  1. Descubrimiento y entendimiento de datos

Este conjunto de datos presenta las locaciones de puntos muestreados en diferentes zonas de Australia, donde se toman valores de temperatura máxima diaria promedio, elevación y coordenadas de georeferenciación en latitud - longitud.

los datos fueron extraídos de la página de Kaggle, la cual es una plataforma web que reune a la comunidad de ciencia de datos más grande del mundo, con más de 536 mil miembros activos en 194 países Temperarura_Australia.

Las variables que se tienen en cuenta para el estudio son las siguientes:

Variable Descripción Tipo Detalle
X Código del registro integer -
lat Coordenadas latitud float -
long Coordenadas longitud float -
elev Elevación de la zona float -
name ciudad-pueblo donde se registró la temperatura string -
Year Año en que se tomaron los datos interger -
Temp Valor de temperatura tomado (grados celcius) float - 
#Lectura librerías
library(leaflet)
library(dplyr)
library(ggplot2)
library(geoR)
library(sf)
library(rgdal) 
library(readxl)
library(raster)
library(readr)
library(kableExtra)
library(ade4)
library(vegan)
library(gstat)
library(sp)

Lectura Base De Datos 

#Lectura Datos
temper <- read.csv("temp.australia.csv")
head(temper)
  • Transformando columnas latitud y longitud a tipo numérico:
temper$lat <- as.numeric(temper$lat)
temper$long <- as.numeric(temper$long)
class(temper$lat)
## [1] "numeric"
class(temper$long)
## [1] "numeric"

Mapa y locaciones muestreadas 

mapa_aus <- leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap)
mapa_aus %>% 
  addMarkers(lat=~lat, lng =~long, data = temper,
                   icon = makeIcon(iconUrl = "temper2.png",
                            iconWidth = 20, iconHeight = 20),
             popup = paste(
               "<h3>Registro climatológico</h3>",
               paste("Temperatura:",
                     round(temper$Temp, 2)),
               paste("Elevación:",
                     round(temper$elev), 2),
               sep = "<br>"
             )
  )

Análisis descriptivo

ggplot(temper,aes(x=long,y=lat,color=Temp))+geom_point()+theme_bw()+
   scale_color_continuous(type = "viridis")

Se puede notar que las temperaturas bajas son muy frecuentadas en la zona sur de Australia y en el norte las temperaturas más cálidas, lo cual puede tener una relación muy significativa con la cercanía a la Antártida.

Datos en formato geoR

# Covertir datos en formato geoR
geodatoslatlong <- as.geodata(temper,coords.col=3:2,data.col=7)
class(geodatoslatlong)
## [1] "geodata"
plot(geodatoslatlong,lowess=TRUE)

gr <- pred_grid(geodatoslatlong$coords, by=20)
points(geodatoslatlong)
points(gr, col=7, pch=19, cex=0.3)

gr0 <- gr[.geoR_inout(gr, geodatoslatlong$coords),]
points(gr0, col=7, pch=19, cex=0.3)

Shapefile Australia

austr <- shapefile("AUS_2021_AUST_GDA2020.shp")
class(austr)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
austr
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 96.81695, 167.998, -43.7405, -9.142163  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## variables   : 6
## names       : AUS_CODE21, AUS_NAME21, CHG_FLAG21, CHG_LBL21,   AREASQKM21,                                        LOCI_URI21 
## value       :        AUS,  Australia,          0, No change, 7688094.9379, http://linked.data.gov.au/dataset/asgsed3/AUS/AUS

Geometría Australia

croq <- geometry(austr)
plot(croq)

  1. Test de autocorrelación espacial de Mantel

Para realizar el test, claramente se deben proyectar los puntos muestreados a UTM, con el fin de llevar las geolocalizaciones al plano, para ésto, se crea una variable llamada zona, la cual especifica la zona en mercator para cada uno de los puntos muestreados. Posteriormente, se genera una función que utiliza a spTransform de la librería rgdal y recorre cada una de las zonas y las coordenadas latitud - longitud para adicionar a los datos las columnas de coordenadas UTM para cada uno de los puntos.

Para la conversión a UTM, Australia cuenta con múltiples zonas de mercator, por ende, se genera una función que recorre las zonas y transforma dichas coordenadas de latitud - longitud a UTM.

Al transformar las coordenadas a UTM, se obtienen valores en metros, pero dado el tamaño de Australia, se decide convertir dichas métricas a kilómetros.

Transformación coordenadas a UTM

get_utm <- function(long, lat, zone, loc){
  points = SpatialPoints(cbind(long, lat), proj4string = CRS("+proj=longlat +datum=WGS84"))
  points_utm = spTransform(points, CRS(paste0("+proj=utm +south=T +zone=",zone[1]," +ellps=WGS84")))
  if (loc == "long") {
    return(coordinates(points_utm)[,1])
  } else if (loc == "lat") {
    return(coordinates(points_utm)[,2])
  }
}

temper %<>% 
  mutate(zone = (floor((long + 180)/6) %% 60) + 1, keep = "all"
  ) %>% 
  group_by(zone) %>% 
  mutate(utm_x = get_utm(long, lat, zone, loc = "long"),
         utm_y = get_utm(long, lat, zone, loc = "lat"))

temper$utm_x <- temper$utm_x/1000
temper$utm_y <- temper$utm_y/1000

Datos finales

head(temper)

Datos representados en el plano 

plot(temper$utm_x,temper$utm_y,ylab="UTM Y",xlab="UTM X")

Realización del test

Se genera Matriz de distancias espaciales entre los lugaes de muestreo

# Matriz de distancias espaciales entre los lugaes de muestreo
geodist <- dist(temper[,11:10], diag=TRUE, upper=TRUE)

Matriz de distancias entre las observaciones registradas en los lugares de muestreo de temperatura

# Matriz de distancias entre las observacioens registradas 
# en los lugares de muestreo de temperatura
dist.temp <- dist(temper[,7], diag=TRUE, upper=TRUE)^2

Se genera un data frame para extraer los valores de las matrices, utilizando la fracción inferior lower

# Creamos un data frame para extraer los valores de las matrices, ya sea la fracción triangular superior o inferior (en este ejemplo, utilizaremos la inferior, es decir: lower)
df2 <- data.frame(Disimilitud=dist.temp[lower.tri(dist.temp)], Distancia=geodist[lower.tri(geodist)])

Gráfico Disimilitud Vs Distancia geográfica (m)

plot <- ggplot(df2, aes(x=Distancia, y=Disimilitud)) + geom_point() +
  theme_test(base_size = 14, base_family = "Times New Roman") +
  labs(x="Distancia geográfica (km)", y="Disimilitud") +
  theme(axis.title = element_text(face="bold"), axis.text = element_text(colour="black"))
plot

Gráfico Disimilitud Vs Distancia geográfica (m) con línea tendencia

plot2 <- ggplot(df2, aes(x=Distancia, y=Disimilitud)) + 
  geom_point() + theme_test(base_size = 14, base_family = "Times New Roman") +
  labs(x="Distancia geográfica (km)", y="Disimilitud") +
  theme(axis.title = element_text(face="bold"), axis.text = element_text(colour="black")) +
  geom_smooth(method=lm)
plot2

Test Mantel

mantel(geodist, dist.temp) # vegan
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = geodist, ydis = dist.temp) 
## 
## Mantel statistic r: 0.7665 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0556 0.0726 0.0842 0.1069 
## Permutation: free
## Number of permutations: 999

Evidentemente se rechaza la hipótesis nula, lo que sugiere que si hay autocorrelación espacial en los datos.

Monte Carlo test

mantel.rtest(geodist, dist.temp, nrepet=999) # ade4
## Monte-Carlo test
## Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
## 
## Observation: 0.7665377 
## 
## Based on 999 replicates
## Simulated p-value: 0.001 
## Alternative hypothesis: greater 
## 
##       Std.Obs   Expectation      Variance 
##  1.857250e+01 -3.274998e-05  1.703583e-03

Con el test de Monte Carlo, se llega a la misma conclusión, que nos entrega evidencia estadística para quedarnos con auocorrelación espacial.

Datos UTM en formato geoR

# Covertir datos en formato geoR
geodatos <- as.geodata(temper,coords.col=10:11,data.col=7)
plot(geodatos,lowess=TRUE)

Estos gráficos, entregan señales muy a favor de la autocorrelación espacial, en el primer gráfico superior izquierdo, se observan clusters y en el gráfico inferior izquierdo una clara tendencia positiva..

Resumen geodatos

summary(geodatos)
## Number of data points: 111 
## 
## Coordinates summary
##        utm_x    utm_y
## min 187.4691 5184.758
## max 785.5417 8830.168
## 
## Distance summary
##         min         max 
##    8.032545 3647.691357 
## 
## Data summary
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  9.815563 20.176630 23.356290 23.928424 27.313460 34.336910
points(geodatos,col = "gray", pt.divide = "equal")

Semivariograma muestral

# Semivariograma muestral
variograma <- variog(geodatos,option="bin",uvec=seq(0,3800,200))
## variog: computing omnidirectional variogram
# Grafico del variograma muestral
plot(variograma)

Test Gráfico (Bandas de confianza)

Destruyendo la correlación espacial a traves de permutaciones aleatorias

# Destruir la correlación espacial a traves de permutaciones aleatorias
geoP2 <- data.frame(geodatos$coords,sample(geodatos$data))
geodatosP2 <- as.geodata(geoP2,coords.col=1:2,data.col=3)
plot(geodatosP2)

Semivariograma muestral

# Semivariograma muestral
variogramaP2 <- variog(geodatosP2,option = "bin",uvec=seq(0,3800,200))
## variog: computing omnidirectional variogram
plot(variograma, col="2")
lines(variogramaP2,type="l")

Bandas de confianza para test gráfico

plot(variograma)
variograma.emv=variog.mc.env(geodatos,obj=variograma)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
lines(variograma.emv,col="blue")

Observaciones

El test gráfico manifiesta un comportamiento esperado, donde la mayoría de puntos caen fuera de las bandas, por lo cual, la variable Temperatura cuenta con autocorrelación espacial.

  1. Ajuste de modelo y predicciones

Considerando que ya se puede ajustar un posible modelo, pues la variable temperatura cuenta con las condiciones de correlación espacial, se procede a ajustar un modelo con el fin de realizar predicciones Kriging en locaciones no muestreadas de Australia.

Determinación de efectos direccionales - anisotropía

plot(variog4(geodatos, trend = "2nd", max.dist = 3647.69), omnidirectional = T)
## 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 variogram

Como se puede observar en el grafico de semivariogramas direccionales, existe anisotropia la cual se da cuando la autocorrelación no es igual en todas las direcciones, en este caso es anisotropia zonal, ya que todos los semivariogramas tienen igual rango pero diferente meseta.

Ajustes

  • Modelo exponencial
exp_fit_fix <- variofit(variograma, cov.model = "exp", fix.nugget = F)
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq  phi    tausq kappa
## initial.value "108.95" "2880" "0"   "0.5"
## status        "est"    "est"  "est" "fix"
## loss value: 1145121.65594699
  • Modelo esférico
sph_fit_fix <- variofit(variograma, cov.model = "sph", fix.nugget = F)
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "72.63" "2880" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 1101980.59861456
  • Modelo gaussiano
gauss_fit_fix <- variofit(variograma, cov.model = "gaussian", fix.nugget = F)
## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq  phi    tausq kappa
## initial.value "145.26" "2304" "0"   "0.5"
## status        "est"    "est"  "est" "fix"
## loss value: 51317.2341546973
  • Comparación de ajuste modelos

Modelo gaussiano:Rojo

plot(variograma,pch=16)
lines(exp_fit_fix,col="green",lwd=4, lty = 1)
lines(sph_fit_fix,col="blue",lwd=4, lty = 3) 
lines(gauss_fit_fix,col="red",lwd=4, lty = 4) 
legend("topleft", c("Gauss","Expo","Spe"),
         col = c("red","blue","green"), lwd = 1, cex = 0.5)

Comparación del SSE de los modelos

Modelo exponencial

(exp_SSQ_fix <- summary(exp_fit_fix)$sum.of.squares)
##    value 
## 520916.5

Modelo esférico

(sph_SSQ_fix <- summary(sph_fit_fix)$sum.of.squares)
##    value 
## 520876.9

Modelo gaussiano

(gauss_SSQ_fix <- summary(gauss_fit_fix)$sum.of.squares)
##    value 
## 11706.03
which.min(list(exp_SSQ_fix,
          sph_SSQ_fix, gauss_SSQ_fix))
## [1] 3
  • Dado que la suma de cuadrados representa la distancia que separa a los valores de los datos de los valores ajustados, podemos inferir que a menor SSE mejor modelo para la prediicción, en este caso el menor SSE lo tiene el modelo gaussiano, éste es el modelo más óptimo para aplicar el krigging, lo que también se evidenció en la gráfica de ajuste con los 3 modelos.

Creación de una grilla de predicción

#xmin:180     xmax: 800    ymin:5000       ymax:9000
prediction_grid <- expand.grid(seq(187.4691, 785.5417, length.out = 100), seq(5184.758, 8830.168, length.out = 100))
plot(prediction_grid)

krig_temp <- krige.conv(geodatos, loc=prediction_grid,
                       krige=krige.control(obj.model=gauss_fit_fix))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(krig_temp,main="Predicción Modelo Gaussiano")
points(geodatos$coords, pch=20) #añadir posiciones datos

contour(krig_temp,add=TRUE,filled = TRUE) #añadir gráfico de contorno

Comparación de modelos

variogramas <- list()
bandas_conf <- list()
ini.vals <- expand.grid(seq(10, 40, l = 10), seq(50000, 150000, l = 10))
model_mco_exps <- list()
model_mco_gaus <- list()
model_mco_spes <- list()
for(i in 1:4){
  variograma <- variog(geodatos, data = geodatos$data, 
                       option = "bin",
                       uvec = seq(0, 3800, length.out = 20),
                       messages = F)
  variogramas[[i]] <- variograma
  
  ## bandas de confianza, si hay o no autocorrelacion 
  variograma.emv <- variog.mc.env(geodatos, data = geodatos$data,
                                  obj = variograma, messages = F)
  bandas_conf[[i]] <- variograma.emv
  
  model_mco_exps[[i]] <- variofit(variogramas[[i]], ini = ini.vals,
                                  cov.model = "exponential",
                                  wei = "npair", min = "optim", messages = F)
  model_mco_gaus[[i]] <- variofit(variogramas[[i]], ini = ini.vals,
                          cov.model = "gaussian", 
                          wei = "npair", min = "optim", nugget = 0, messages = F)
  model_mco_spes[[i]] <- variofit(variogramas[[i]], ini = ini.vals,
                       cov.model = "spheric", fix.nug = TRUE,
                       wei = "npair", min = "optim", messages = F)
}
Titulos <- c("Temperatura")
  plot(variograma, envelope = bandas_conf[[i]], 
       pch = 19, las = 1, main = Titulos,
       xlab = "Distancia (km)", ylab = "Semivarianza")
  lines(model_mco_exps[[i]], col="blue")
  lines(model_mco_gaus[[i]], col="red")
  lines(model_mco_spes[[i]], col="green")
  legend("topleft", c("Gauss","Expo","Spe"),
         col = c("red","blue","green"), lwd = 1, cex = 0.5)

Aqui podemos ratificar como el modelo gaussiano es el más optimo para la predicción, además de ser el que mejor se ajusta a los puntos del variograma, tambien verificamos que existe autocorrelación espacial a ver como los puntos del semivariograma se encuentran en su mayoria fuera de las bandas de confianza.

raster.ko <- list(map_image_gaus.ko, map_image_exps.ko, map_image_spes.ko)
vars.ko <- list(map_image_gausvar.ko, map_image_expsvar.ko, map_image_spesvar.ko)
nombres_modelos <- c("Modelo gaussiano", "Modelo exponencial",
                     "Modelo esférico")
modelos <- list(model_mco_gaus, model_mco_exps, model_mco_spes)
xvs <- c()
labels <- c()
maps.plot <- function(i, j, plot.var = F) {
  #par(mfrow = c(1, 2))
    label <- paste(Titulos[i], "\n",nombres_modelos[j])
    if (!plot.var) {
      plot(raster.ko[[j]][[i]],
         main = paste(label, "(predicción)"), cex.main = 0.8,
         col = viridis::inferno(20))
    } else {
      plot(vars.ko[[j]][[i]], main = paste(label, "(desviación estándar)"),
         cex.main = 0.8, col = viridis::inferno(20))
    }
    xv <- xvalid(geodatos, data = geodatos$data,
               model = modelos[[j]][[i]])$error^2
    xvs <- append(xvs, mean(xv))
    #labels <- append(labels, str_replace(label, " \\n ", " "))
}

Modelo Gaussiano

maps.plot(1, 1)

Desviación modelo Gaussiano

maps.plot(1, 1, T)

Modelo exponencial

maps.plot(1, 2)

Desviación modelo exponencial

maps.plot(1, 2, T)

Modelo esférico

maps.plot(1, 3)

Desviación modelo esférico

maps.plot(1, 3, T)

*Una vez más, se encuentra que el mejor modelo es el gaussiano, pues los modelos esférico y exponencial, cuentan con menos variación únicamente en los puntos muestreados.

Conclusiones

Predicción

summary(geodatoslatlong)
## Number of data points: 111 
## 
## Coordinates summary
##       long    lat
## min 113.67 -43.49
## max 153.47 -10.58
## 
## Distance summary
##       min       max 
##  0.222036 39.948813 
## 
## Data summary
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  9.815563 20.176630 23.356290 23.928424 27.313460 34.336910
variogramalatlong=variog(geodatoslatlong,option = "bin")
## variog: computing omnidirectional variogram
gauss_fit_fix_latlong <- variofit(variogramalatlong, cov.model = "gaussian", fix.nugget = F)
## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi     tausq  kappa
## initial.value "36.29" "12.29" "4.84" "0.5"
## status        "est"   "est"   "est"  "fix"
## loss value: 130720.717504019
min(geodatoslatlong$coords[,1])  #x
## [1] 113.67
max(geodatoslatlong$coords[,1])
## [1] 153.47
min(geodatoslatlong$coords[,2])
## [1] -43.49
max(geodatoslatlong$coords[,2])
## [1] -10.58
geodatos_grid=expand.grid(lat=seq(113.67,153.47,l=100),long=seq(-43.49,-10.58,l=100))

plot(geodatos_grid)
points(temper[,c(3,2)],col="red",pch=16)

krig_temp_latlong <- krige.conv(geodatoslatlong, loc=geodatos_grid,
                       krige=krige.control(obj.model=gauss_fit_fix_latlong))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,1))
image(krig_temp_latlong, main="kriging Predict", xlab="Lat", ylab="Long")

contour(krig_temp_latlong,main="kriging Predict", add=TRUE, drawlabels=TRUE,filled=TRUE)

Desviación

image(krig_temp_latlong, main="kriging StDv Predicted",val=sqrt(krig_temp_latlong$krige.var), xlab="East", ylab="North")

contour(krig_temp_latlong,main="kriging StDv Predict",val=sqrt(krig_temp_latlong$krige.var), add=TRUE, drawlabels=TRUE,filled=TRUE)

Mapa final Kriging temperatura de Australia

mat_imagen=data.frame(geodatos_grid,krig_temp_latlong$predict)
temp_m1_raster=rasterFromXYZ(mat_imagen)

temp_m1_raster_austr=mask(temp_m1_raster,austr)
plot(temp_m1_raster_austr)
plot(austr,add=TRUE)

crs(temp_m1_raster_austr)=crs(austr)
crs(temp_m1_raster_austr)=CRS("+init=epsg:4326")

pal <- colorNumeric(c("blue", "green", "yellow", "orange", "red"), values(temp_m1_raster_austr),
  na.color = "transparent")


leaflet() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addRasterImage(temp_m1_raster_austr,colors = pal)
  • Como pudimos observar entre más al sur sea la ubicación en Australia, es más baja la temperatura, dado que la antártida está realmente cerca, y entre más al norte sea la ubicación, las temperaturas tienden a ser más altas, puesto que en dichas zonas se encuentran terrenos más áridos y zonas desérticas.