Niveles de concentración de NO2 en el aire en la España peninsular en 2018

Introducción

El presente documento tiene como objetivo mostar el proceso que se ha seguido para determinar en qué zonas de la España peninsular se concentra un nivel mayor de NO2 y en que partes la presencia de este gas contaminante es reducido. Para ello, utilizaremos una base de datos que contiene los valores de esta sustancia recogidos en estaciones de medición situadas en diferentes puntos del país. Nos ayudaremos del mapa de España a partir del cuál crearemos una maya que nos permitirá interpolar los valores para obtener los niveles de concentración correspondientes a todas las zonas de la Península y dar respuesta a la cuestión que queremos resolver.

El proceso se ha hecho para el año 2018. Para la obtención de los mapas que se corresponden con los años 2017, 2016, 2015 y 2014 se ha seguido el mismo proceso.

Cargamos las librerias

librerias <- c("tidyverse",
               "dplyr",
               "readxl",
               "data.table",
               "knitr",
               "tables",
               "kableExtra",
               "sp",
               "sf",
               "gstat",
               "osmdata",
               "ggmap",
               "stars",
               "gstat",
               "mapview",
               "gridExtra"
               )

for (pkg in librerias)
  {
  if (!(pkg %in% installed.packages()))
      install.packages(pkg, repos = "http://cran.r-project.org")
 
  library(pkg, character.only = TRUE)
}

Cargamos y limpiamos los datos

Los datos con los que vamos a trabajar han sido obtenidos de la European Environmental Agency (EEA). Los estados miembros europeos proporcionan medidas de la calidad del aire a esta agencia. En este caso, nos interesa conocer el nivel de concentración de NO2 en el aire en la Península Ibérica en el año 2018 .(Link para descargar los datos:https://www.eea.europa.eu/data-and-maps/data/air-pollutant-concentrations-at-station/air-pollutant-concentrations-2018/air-pollutant-concentrations-2018-compared )

Estos datos han sido tratados antes de ser cargados en R. Se han seleccionado las observaciones que se corresponden con España y se han eliminado los datos atribuibles a las provincias de Ceuta, Melilla, Islas Baleares e Islas Canarias.

setwd("C:/Users/sanch/OneDrive/Escritorio/BIA/2021-2022/Datos_Espaciales_y_Espacio_Temporales/trabajo_final")

ES_NO2 <- read_excel("C:/Users/sanch/OneDrive/Escritorio/BIA/2021-2022/Datos_Espaciales_y_Espacio_Temporales/trabajo_final/NO2_por_anyos.xlsx", 
                     sheet = "2018")

NA y datos duplicados

También eliminamos los NA y los datos que puedan estar duplicados. Este último paso es importante ya que al no hacerlo, los datos no eran del todo correctos y, además de tardar una importante cantidad de tiempo en ejecutar los diferentes trozos de código, no era posible llevar a cabo la interpolación con el Kriging Ordinario. Esto se debe a que los datos están medidos en diferentes momentos del tiempo por lo que hay valores duplicados

ES_NO2<- ES_NO2%>%
  filter(!is.na(ES_NO2$AirPollutionLevel) == TRUE )

ES_NO2 <- ES_NO2 %>% distinct(Longitude, .keep_all = TRUE)

Detección de valores atípicos (outliers)

Para identidicar valores atípicos visualmente, representamos un diagrama de caja y bigotes (boxplot). Para ello hacemos uso de la función boxplot:

boxplot(ES_NO2$AirPollutionLevel, # Datos
        horizontal = FALSE, # Horizontal o vertical
        lwd = 2, # Lines width
        col = "cadetblue3", # Color
        xlab = "",  # Etiqueta eje X
        ylab = "",  # Etiqueta eje Y
        main = "Boxplot valores NO2", # Título
        notch = F, # Añade intervalos de confianza para la mediana
        border = "black",  # Color del borde del boxplot
        outpch = 25,       # Símbolo para los outliers
        outbg = "red",   # Color de los datos atípicos
        whiskcol = "black", # Color de los bigotes
        whisklty = 2,      # Tipo de línea para los bigotes
        lty = 1) # Tipo de línea (caja y mediana)

Podemos apreciar la existencia de valores atipicos representados en color rojo. A continuación procedemos a eliminarlos y volver a hacer la representación para comprobar que ya no existen.

Vemos que se consideran valores outliers aquellos que son mayores que 44, por lo que, con la aydua de nuevo de la función filter, eliminamos dichos valores.

ES_NO2<- ES_NO2%>%
  filter(!AirPollutionLevel > 44 )

boxplot(ES_NO2$AirPollutionLevel, # Datos
        horizontal = FALSE, # Horizontal o vertical
        lwd = 2, # Lines width
        col = "cadetblue3", # Color
        xlab = "",  # Etiqueta eje X
        ylab = "",  # Etiqueta eje Y
        main = "Boxplot valores NO2", # Título
        notch = F, # Añade intervalos de confianza para la mediana
        border = "black",  # Color del borde del boxplot
        outpch = 25,       # Símbolo para los outliers
        outbg = "red",   # Color de los datos atípicos
        whiskcol = "black", # Color de los bigotes
        whisklty = 2,      # Tipo de línea para los bigotes
        lty = 1) # Tipo de línea (caja y mediana)

Ya no existen valores atípicos en nuestra base de datos. Veamos como queda el dataset de trabajo mostrando las 10 primeras líneas:

kbl(head(ES_NO2,10)) %>%
  kable_styling(c("striped","bordered"))

…1	Country	AirQualityStationEoICode	AQStationName	AirPollutant	AirPollutionLevel	UnitOfAirpollutionLevel	AirQualityStationType	AirQualityStationArea	Longitude	Latitude	Altitude
Coruña	Spain	ES1960A	A CABANA	NO2	7.301683	ug/m3	Industrial	suburban	-8.252200	43.4916	33
Zaragoza	Spain	ES0588A	ABANTO	NO2	17.935069	ug/m3	Industrial	suburban	-3.074140	43.3205	136
Alava	Spain	ES1544A	AGURAIN	NO2	12.072474	ug/m3	Background	suburban	-2.393700	42.8490	594
Alicante	Spain	ES1968A	ALACANT - RABASSA	NO2	14.089008	ug/m3	Industrial	suburban	-0.513889	38.3511	20
Alicante/Alacant	Spain	ES1635A	ALACANT-EL PLÁ	NO2	21.510279	ug/m3	Traffic	urban	-0.471944	38.3594	45
Alicante/Alacant	Spain	ES1915A	ALACANT-FLORIDA-BABEL	NO2	23.783469	ug/m3	Background	urban	-0.506667	38.3403	55
Zaragoza	Spain	ES1418A	ALAGÓN	NO2	16.934297	ug/m3	Traffic	suburban	-1.143330	41.7628	235
Albacete	Spain	ES1535A	ALBACETE	NO2	13.538145	ug/m3	Background	suburban	-1.852100	38.9793	686
Valencia	Spain	ES1911A	ALBALAT DELS TARONGERS	NO2	7.267711	ug/m3	Industrial	suburban	-0.336667	39.7053	320
Sevilla	Spain	ES1640A	ALCALÁ DE GUADAIRA	NO2	17.138363	ug/m3	Background	urban	-5.833720	37.3425	68

Representación de histogramas para ver la distribución de la variable de interés

Vemos la distribución de los valores que toma la variable de interés:

hist(ES_NO2$AirPollutionLevel,
     col= "cadetblue3",
     main="Histograma de NO2",
     xlab = "NO2 (ug/m3)",
     ylab = "Frecuencia")

La gran mayoría de los valores de NO2 se concentran entre el 5 y el 15. Creamos una nueva variable en la base de datos que contendrá el logaritmo neperiano de esta variable para que su distribución se asemeje a una normal y de esta manera evitar futuros problemas a la hora de llevar a cabo la interpolación. Asimismo, eliminamos aquellos valores que puedan ser negativos, ya no solo para no incurrur en dificultades sino porque no tiene mucho sentido que el nivel de concenctración de NO2 en el aire sea negativo, como poco, debería ser 0.

ES_NO2$VALORLOG = log(ES_NO2$AirPollutionLevel)
hist(ES_NO2$VALORLOG,
     col= "cadetblue3",
     main="Histograma de NO2 transformada con log Ln",
     xlab = "NO2 (ug/m3)",
     ylab = "Frecuencia")

ES_NO2<- ES_NO2%>%
  filter(!VALORLOG < 0 )

hist(ES_NO2$VALORLOG,
     col= "cadetblue3",
     main="Histograma de NO2 transformada con log Ln",
     xlab = "NO2 (ug/m3)",
     ylab = "Frecuencia")

Creación objeto espacial

Convertimos nuestra base de datos en un objeto espacial. Utilizamos para ello el sistema de coordenadas de referencia “EPSG:4258” (https://epsg.io/4258).

Nuestra base de datos utilza el sistema de coordenadas de referencia “EPSG:4979” que es el que tiene un área de uso mundial, por eso le indicamos que tome ese crs a la hora de convertir la base de datos en un objeto espacial y posteriormente lo transformamos.

library(sf)
crs <- st_crs("EPSG:4258")
no2.sf <- st_as_sf(ES_NO2, coords = c("Longitude", "Latitude"), crs = "EPSG:4979") %>%
  st_transform(crs)

Si representamos los puntos de datos podemos apreciar el mapa de España:

plot(no2.sf$geometry,asp=1,pch=1, col="cadetblue4")

Cargamos y trabajamos con el mapa de España

Una vez tenemos los datos, vamos a utilizar el mapa de España agrupando por las diferentes Comunidades Autónomas. Dicho mapa nos resultará útil para la creación del grid a partir del cuál interpolaremos más adelante. Asimismo, utilizaremos el contorno del mismo para dar contexto a los valores obtenidos en la interpolación y así poder identificar con mayor facilidad a que zona se corresponde cada dato obtenido.

Para la obtención del plano hacemos uso de la libreria mapSpain y, como bien hemos explicado anteriormente, utilizaremos el de la Península Ibérica por lo que eliminamos las observaciones cuyo código de comundiad autónoma se corresponde con alguno de los siguientes: 04, 05, 18, 19. En otras palabras, eliminamos aquellos datos que se atribuyen a Ceuta, Melilla, Canarias y Baleares.

library(mapSpain)
mapa_ca<-mapSpain::esp_get_prov()
mapa_ca<- mapa_ca%>%
  filter(!codauto %in% c("04","05","18","19")) %>%
  group_by(codauto) %>%
  summarise()

plot(mapa_ca$geometry)

Transformamos la variable del mapa en un objeto espacial y le asignamos el sistema de coordenads de referencia con el que vamos a trabajar que es el mismo que hemos asignado al objeto que contiene los datos, el “EPSG:4258”.

mapa.sf <- st_as_sf(mapa_ca, coords = geometry) 
st_crs(mapa.sf) <- crs



ggplot() + 
  geom_sf(data = mapa.sf) +  
  geom_sf(data = no2.sf, color="cadetblue4")+#, mapping = aes(col = VALORLOG))+
  labs(title = "Localización de las estaciones de medición de NO2",
       subtitle="España 2018",
       caption = "Fuente: European Environment Agency")+
  theme_void()

Cálculo de la semivarianza

La dependencia espacial local significa que cuanto más cerca estén los puntos en el espacio geográfico, más cerca también lo están en el espacio de atributos. Esto se llama autocorrelación. Una variable puede correlacionarse consigo misma, dependiendo la fuerza de la correlación de la distancia de separación entre puntos. También es importante concocer el concepto de semivarianza, que es una medida de la similitud que existe entre observaciones situadas a una determinada distancia. Cada par de puntos, estará separado por una distancia en el espacio geográfico y una semivarianza en el espacio de atributos. Fuentes: · https://www.redalyc.org/pdf/360/36014577008.pdf
· https://www.uv.es/pegivir/BIA/sesion_practica_7.html

Podemos calcular el número de pares de puntos que hay en el conjunto de datos:

n <- length(no2.sf$VALORLOG)
n*(n-1)/2

## [1] 84666

Calculamos la distancia y la semivarianza para los primeros dos puntos del dataset:

dim(st_coordinates(no2.sf))

## [1] 412   2

(zi <- st_coordinates(no2.sf)[1,]) #coordenadas del primer punto

##       X       Y 
## -8.2522 43.4916

(zj <- st_coordinates(no2.sf)[2,])  #coordenadas del segundo punto

##        X        Y 
## -3.07414 43.32050

(sep <- dist(st_coordinates(no2.sf)[1:2,]))  #distancia entre ambos puntos

##          1
## 2 5.180886

La distancia entre ambos puntos es de 8.4 metros.

(gamma <- 0.5 * (no2.sf$VALORLOG[1] - no2.sf$VALORLOG[2])^2) # semivarianza, unidades log(mg kg-1)^2

## [1] 0.4037887

Variograma experimental

Ahora graficaremos el variograma experimental, teniendo en cuenta las concentraciones de NO2 (en log). Recordemos que el variograma es una representación gráfica que agrupa los pares de puntos por intervalos de distancia, calculando la media de todas las semivarianzas para cada intervalo. En otras palabras, el variograma mide la variación existente entre la distancia y la variabilidad. Por defecto, la función variogram() usa un valor de cutoff igual a la tercera parte de la diagonal del boundg box del conjunto de puntos, y lo divide en 15 partes para obtener el valor de width. La fórmula que usaremos es VALORLOG~1. El 1 quiere decir que VALORLOG (la variable dependiente) solo depende de sí misma (“autocorrelación”). Es decir, la presencia de NO2 en el aire no depende de otra variable.

(ve <- variogram(VALORLOG~1, no2.sf)) 

##      np      dist     gamma dir.hor dir.ver   id
## 1  1924  14.63110 0.1422311       0       0 var1
## 2  1536  43.95134 0.4773076       0       0 var1
## 3  1876  74.30104 0.5319810       0       0 var1
## 4  1738 102.85488 0.7385145       0       0 var1
## 5  1964 132.45715 0.7033644       0       0 var1
## 6  2211 162.03678 0.6236147       0       0 var1
## 7  2413 192.13888 0.6564020       0       0 var1
## 8  2705 220.66761 0.6169243       0       0 var1
## 9  2967 250.38872 0.5641465       0       0 var1
## 10 3506 279.95074 0.6463895       0       0 var1
## 11 4470 308.81773 0.5519221       0       0 var1
## 12 4188 338.07885 0.6182939       0       0 var1
## 13 3913 367.02457 0.6308508       0       0 var1
## 14 4237 396.59295 0.6309867       0       0 var1
## 15 3709 425.71211 0.5901305       0       0 var1

plot(ve, plot.numbers = TRUE, pch = 20,
     xlab = "Distancia",
     ylab = "Semivarianza",
     col= "cadetblue3")

Ahora repetimos el paso anterior pero esta vez en lugar de usar los valores por defecto, fijamos los intervalos (width) en una distancia de 30 metros, hasta una distancia máxima (cutoff) de 500 metros.

(ve <- variogram(VALORLOG ~ 1, no2.sf, cutoff = 500, width = 30))

##      np      dist     gamma dir.hor dir.ver   id
## 1  1951  14.83899 0.1435275       0       0 var1
## 2  1568  44.77696 0.4899381       0       0 var1
## 3  1930  75.61863 0.5205324       0       0 var1
## 4  1778 105.10569 0.7632675       0       0 var1
## 5  2013 135.10713 0.6763529       0       0 var1
## 6  2285 165.19115 0.6349109       0       0 var1
## 7  2486 195.87667 0.6670568       0       0 var1
## 8  2801 224.95162 0.5908225       0       0 var1
## 9  3062 255.25967 0.5982555       0       0 var1
## 10 3818 285.55265 0.6004536       0       0 var1
## 11 4558 314.66515 0.5729544       0       0 var1
## 12 4312 345.09318 0.6442019       0       0 var1
## 13 4030 375.41595 0.6095890       0       0 var1
## 14 4218 404.85771 0.6367131       0       0 var1
## 15 3589 434.62169 0.5660850       0       0 var1
## 16 3641 465.80155 0.6092482       0       0 var1
## 17 3146 490.33361 0.5514807       0       0 var1

plot(ve, plot.numbers = TRUE, pch = 20,
     xlab = "Distancia",
     ylab = "Semivarianza",
     col= "cadetblue3")

Calculamos la diagonal del bounding box

min <- st_point(c(st_bbox(no2.sf)[1],st_bbox(no2.sf)[2]))
max <- st_point(c(st_bbox(no2.sf)[3],st_bbox(no2.sf)[4]))
st_distance(min,max)

##          [,1]
## [1,] 14.17983

Pueden representarse los valores de la semivarianza frente a las distancias, obteniendo una nube de puntos (nube del variograma) como esta:

plot(variogram(VALORLOG ~ 1, no2.sf, cloud=TRUE), pch = 20,
     xlab = "Distancia",
     ylab = "Semivarianza",
     col= "cadetblue3")

Se puede ver que a medida que aumenta la distancia, los valores que toma la semivarianza se van “despegando” del fondo (dejan de concentrarse en los valores más bajos) y van tomando valores más altos. No obstante, esta nube aporta, en principio, poca información, pero puede resumirse agrupando los pares de puntos por intervalos de distancia, y calculando la media de todas las semivarianzas en cada intervalo, como veremos más adelante.

También podemos mostrar con la función variogram() los pares de puntos que están más próximos/lejanos entre sí:

vc <- variogram(VALORLOG ~ 1, no2.sf, cloud=TRUE)
min(vc$dist)

## [1] 0.7222379

max(vc$dist)

## [1] 440.9152

vc <- vc %>% arrange(dist)
head(vc)

##        dist       gamma dir.hor dir.ver   id left right
## 1 0.7222379 0.051326991       0       0 var1  185   146
## 2 0.8298581 0.030524949       0       0 var1  137    83
## 3 0.8782850 0.059199794       0       0 var1  181    93
## 4 0.9318389 0.002764745       0       0 var1  156    52
## 5 0.9947059 0.048959306       0       0 var1  183    82
## 6 1.0165444 0.032729850       0       0 var1  401   278

Y verificar que distancia y semivarianza tienen una relación positiva:

ggplot(data = vc, aes(x=dist, y=gamma)) +
  geom_point(color="cadetblue3") +
  geom_smooth(method='lm', formula=y~x,
   col="cadetblue4")+
  xlab("Distancia") + ylab("Gamma") + 
  theme_classic()

vemos que hay una relacion positiva entre la distancia y gamma. Una pendiente positiva nos dice que a mayor distancia mayor es la variabilidad mientras que a menos distacnia, más parecidos son los puntos.

Determinación de efectos direccionales: Anisotropía-Isotropía

Un aspecto importante a tener en cuenta a la hora de hacer una interpolación espacial es la determinación de efectos direccionales. En caso de existir anistropía espacial, influirá bastante la dirección de la interpolación, pues, la estructura espacial será diferente en distintas direcciones. Por el contrario, en caso de exisitir isotpía espacial, la dirección en la que se haga la interpolación no influirá en los resultados obtenidos ya que la estructura espacial será la misma en diferentes direcciones. No se procederá, por tanto, de la misma manera en un caso que en otro.

“En caso de detectar anisotropía espacial se definiría una elipse cuyos ejes mayor y menor reflejasen adecuadamente la misma. En caso de no detectar anisotropía se optaría por el modelo isotrópico,en el que sólo se tendría en cuenta la variación de los valores en función de la distancia.”

• https://dialnet.unirioja.es/descarga/articulo/4757146.pdf

Los gases son considerados materiales isótropos.

“Un material se considera isotrópico si sus propiedades no varían con la dirección. Los materiales isotrópicos, por lo tanto, tienen un módulo elástico, un coeficiente de Poisson, una conductividad térmica, un coeficiente de expansión térmica, etc. idénticos en todas direcciones. El término isotérmico algunas veces se utiliza para indicar materiales sin direcciones preferidas para coeficientes de expansión térmica.”

• https://help.solidworks.com/2013/spanish/solidworks/cworks/c_isotropic_orthotropic_materials.htm#.

No obstante, vamos a ajustar un variograma en diferentes grados de dirección, concretamente en 0º,45º,90º y 135º para comprobar si nuestra variable de estudio es isótropa.

Para ello, usamos la función variogram de nuevo pero con la peculariedad de que esta vez añadimos el argumento alpha= en el que le damos un vector con las direcciones en plano en grados positivos. Posteriormente, dibujamos los variogramas obtenidos para cada una de las direcciones. Esto lo haremos con la función ggplot y, dentro de esta, usaremos facet_wrap para que muestre cada uno de los variogramas.

vgm.aniso <- variogram(VALORLOG ~ 1, no2.sf, alpha = c(0, 45, 90, 135))


(no2.aniso.plot<- ggplot(aes(x = dist, y = gamma), data = vgm.aniso)+ # dibujamos los variogramas
                 geom_point(color="cadetblue4")+ # queremos que muestre los puntos y en ese color
                 facet_wrap(~dir.hor)+ # queremos que divida en función de la dirección
                 scale_y_continuous(limits = c(0, 1.3))+ # ajustamos el eje de las y
                 xlab("Distancia") + ylab("Gamma") + # ponemos nombre a los ejes
                 labs(title = "Determinación de efectos direccionales") + # ponemos titulo
                 theme_bw()+
                 theme(strip.background = element_rect(fill="cadetblue3"))) # queremos pintar los títulos de ese color

La correlación espacial se mantiene igual de fuerte en las diferentes direcciones, por lo que confirmamos que el gas, en este caso el NO2, es un material isótropo.

Variograma empírico

Ahora queremos conocer cómo se relacionan las semivarianzas respecto de las distancias, en promedio. Para saberlo, debemos comparar nuestro viariograma experimiental con los distintos variogramas empíricos existentes. Primero, veamos los diferentes tipos:

print(show.vgms())

Una vez tenemos claro qué modelo se asemeja más a nuestro variograma experimental, lo generamos. En nuestro caso, usaremos el modelo esférico: “Sph”. Para ello, debemos introducir los parámetros del que conforman el variograma empírico, es decir, sill, range y nugget:

• sill: representa la máxima variabilidad en ausencia de dependencia espacial.
• range: representa la máxima distancia a partir de la cual existe dependencia espacial. Es el valor en el que se alcanza la máxima varianza, o a partir del cual ya presenta una tendencia asintótica.
• nugget: conforme la distancia tiende a cero, el valor de la semivarianza tiende a este valor. Representa aquella variabilidad que no puede explicarse mediante la estructura espacial.

vm <- vgm(model="Sph",
          psill=0.51,
          range=95,
          nugget=0.1)

plot(ve, pl=T, model=vm, pch=20, xlab = "distance", ylab = "semivariance",
     col= "cadetblue3")

Vemos que a partir de una determinada distancia, la variabilidad no cambia, la distancia ya no influye en la variabilidad de los valores.

Pero puede que no veamos de forma clara los valores, o no hayamos detectado bien el modelo del variograma empírico. Esto mismo lo podemos hacer de forma automática:

(vmf <- fit.variogram(ve, vm))

##   model     psill    range
## 1   Nug 0.0000000  0.00000
## 2   Sph 0.6418862 93.74466

plot(ve, pl=T, model=vmf, pch=20, xlab = "distance", ylab = "semivariance")

Podemos afirmar, por tanto, que nuestro conjunto de datos sigue un modelo esférico.

Creamos el GRID

Para interpolar primero debemos decidir dónde lo queremos hacer. Normalmente esto se suele realizar en un grid (una red o cuadrícula) que cubre el área de interés. En nuestro caso, queremos interpolar en el área de España, por lo que, partiendo del contorno del país, creamos un grid uniforme con 10km de celdas de grid sobre la región:

st_bbox(mapa.sf) %>% # seleccionamos el contorno del mapa de España
  st_as_stars(mapa.sf = 10000) %>% # creamos un grid uniforme con 10km de celdas de grid sobre España
  st_crop(mapa.sf) -> grd # creamos el objeto grd que contendrá el grid
plot(grd, main="GRID", col="cadetblue3")

Interpolación

Ahora vamos a practicar la interpolación espacial, es decir, calcular valores desconocidos a partir de una muestra. El resultado que vamos a obtener se conoce habitualmente por superficie estadística, una superficie “continua” (capa raster) con valores interpolados a partir de otros conocidos.

Interpolación Polígonos de Thiessen

Este método de interpolación es considerado como método no-geoestadístico. Al igual que el IDW no requiere ningún modelado especial de relaciones espaciales. En este se producen cambios abruptos de bordes y solo utiliza un punto para cada predicción.

thiessen <- krige(VALORLOG ~ 1, no2.sf, grd, nmax = 1)

## [inverse distance weighted interpolation]

plot(thiessen)

Esta representación gráfica permite hacernos una idea de las zonas en las que el nivel de NO2 es mayor. No obstante, no resulta muy visual y hay zonas en las que puede ser complicado determinar a que comunidad autónoma pertenece. Por tanto, representamos gráficamente con ggplot la predicción obtenida y añadiremos una capa con el contorno del mapa de España y sus comunidades autónomas así como una paleta de colores:

ggplot() + geom_stars(data = thiessen, aes(fill = var1.pred, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción Poligonos de Thiessen log NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

Resulta difícil dar una buena interpretación de la predicción obtenida ya que la presencia de tantas figuras geométricas (polígonos) dentro del mapa hace que sea un poco caótico y algo impreciso.

Interpolación IDW

La interpolación mediante distancia inversa ponderada (IDW) determina los valores desconocidos a partir de una combinación ponderada de un conjunto de puntos de muestra. Este método presupone que el valor de la variable observada disminuye conforme aumenta la distancia respecto del punto muestreado.

Quizá este sea el método de interpolación más simple ya que no requiere ningún modelado especial de relaciones espaciales.

Utilizamos la función idw:

i = idw(VALORLOG~1, no2.sf, grd)

## [inverse distance weighted interpolation]

ggplot() + geom_stars(data = i, aes(fill = var1.pred, x = x, y = y)) + # predicción IDW
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) + # contorno del mapa de España
  xlab("Longitud") + ylab("Latitud") + # damos nombre a los ejes aunque no se verán 
  labs(title = "Predicción IDW log NO2 (ug/m3)") + # ponemos un título al mapa
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+ # añadimos color y ponemos el fondo blanco (transparente)
  theme_void() # utilizamos este tema por lo que no aparecerán los nombres de los ejes

Observamos que los niveles de concentración de NO2 en el aire más elevados se situan en Madrid y Barcelona, los centros urbanos más grandes de España. Estas provincias destacan además por concentrar valores elevados de este gas contaminante en todas las áreas. También se aprecian concentraciones elevadas en ciudades como Valencia, Bilbao, Girona y Ourense, así como algunas áreas puntuales en la Comunidad Autónoma de Sevilla.

Interpolación Kriging Ordinario

Otro método de interpolación muy utilizado es el kriging ordinario (ordinary kriging). La principal ventaja que tiene este método respecto al anterior es que los pesos asignados a cada valor de la muestra son los óptimos, teniendo en cuenta el variograma (IDW solo tiene en cuenta la distancia).

ok <- krige(formula = VALORLOG ~ 1,
            locations = no2.sf,
            newdata = grd,
            model = vmf)

## [using ordinary kriging]

ggplot() + geom_stars(data = ok, aes(fill = var1.pred, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción Kriging Ordinario log NO2 (ug/m3) 03/04/2022") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

Observando el gráfico podemos identificar que las zonas con mayores concentraciones de NO2 están definidas de una manera más homogénea. Se trata de Madrid, área en la mayor es la concentración de NO2, además de Barcelona, Valencia, Bilbao, Girona y Ourense, así como algunas áreas puntuales en la Comunidad Autónoma de Sevilla.

En este caso, también podemos representar las predicciones de las varianzas:

ggplot() + geom_stars(data = ok, aes(fill = var1.var, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción varianzas Kriging Ordinario log NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

Se puede observar que, en general, las zonas con predicciones de las varianzas más bajas están en las proximidades al mar, es decir, en las zonas costeras de las comunidades autónomas que son colindantes con el mar mediterráneo, el mar cantábrico y el océano atlántico.
También cabe mencionar Madrid y algunas zonas del interior de Andalucía, Extremadura, Aragón, Castilla-La Mancha y Castilla y León.

Puede ser útil visualizar las predicciones y las observaciones en un mismo gráfico. Para ello añadimos una capa con las observaciones las cuales estarán representadas con un círculo cuya grandaria representa el nivel de concentración. A mayor nivel de concentración, mayor es el tamaño del círculo.

La línea con la capa que añadimos es: geom_sf(data = no2.sf, pch=1,col="white", cex=12*no2.sf$AirPollutionLevel/max(no2.sf$AirPollutionLevel))

ggplot() + geom_stars(data = ok, aes(fill = var1.pred, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  geom_sf(data = no2.sf, pch=1,col="white",                                                   cex=5*no2.sf$AirPollutionLevel/max(no2.sf$AirPollutionLevel))+
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción Kriging Ordinario log NO2 (ug/m3) + NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

Interpolación Kriging Simple

Una vez hemos hecho la interpolación con el método Kriging Ordinario, llevar a cabo un Kriging Simple es bastante fácil. Este utiliza los mismos argumentos utilizados en el OK con la diferencia de que tiene un argumento adicional, beta.

Si usamos la ayuda de lafunción krige nos dice que si no hay variables independientes (como es el caso) “el modelo solo contiene un intercepto y beta debería ser la media del Kriging simple”. En nuestro caso, como estamos trabajando con el logaritmo usaremos una media igual a 1. Obtenemos lo siguiente:

sk<-krige(VALORLOG~ 1, 
      loc=no2.sf,        
      newdata=grd,     
      model = vmf,     
      beta = 1)  

## [using simple kriging]

ggplot() + geom_stars(data = sk, aes(fill = var1.pred, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción Kriging Simple log NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

Observando el gráfico podemos identificar que las zonas con mayores concentraciones de NO2 están definidas de una manera más homogénea.

Validación cruzada

Otra de las características del método kriging es que permite usar el mismo conjunto de datos tanto para modelar y predecir como para evaluar dichas predicciones. Esto se llama validación cruzada. Consiste en predecir los valores en las localizaciones observadas a partir de los datos generados anteriormente. De este modo podemos comparar las predicciones con la realidad. Esta técnica la podemos aplicar también a otro tipo de interpolaciones, como la IDW, mediante la función krige.cv() del paquete gstat:

idw.cv <- krige.cv(formula = log10(AirPollutionLevel) ~ 1, locations = no2.sf)
head(idw.cv)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -8.2522 ymin: 38.3403 xmax: -0.471944 ymax: 43.4916
## Geodetic CRS:  ETRS89
##   var1.pred var1.var observed    residual zscore fold                  geometry
## 1  1.113660       NA 0.863423 -0.25023693     NA    1   POINT (-8.2522 43.4916)
## 2  1.104284       NA 1.253703  0.14941934     NA    2  POINT (-3.07414 43.3205)
## 3  1.167544       NA 1.081796 -0.08574748     NA    3    POINT (-2.3937 42.849)
## 4  1.366619       NA 1.148880 -0.21773881     NA    4 POINT (-0.513889 38.3511)
## 5  1.255236       NA 1.332646  0.07741023     NA    5 POINT (-0.471944 38.3594)
## 6  1.169909       NA 1.376275  0.20636578     NA    6 POINT (-0.506667 38.3403)

thiessen.cv <- krige.cv(formula = log10(AirPollutionLevel) ~ 1, locations = no2.sf, nmax = 1)
head(thiessen.cv)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -8.2522 ymin: 38.3403 xmax: -0.471944 ymax: 43.4916
## Geodetic CRS:  ETRS89
##   var1.pred var1.var observed    residual zscore fold                  geometry
## 1 1.1189932       NA 0.863423 -0.25557025     NA    1   POINT (-8.2522 43.4916)
## 2 0.9336456       NA 1.253703  0.32005746     NA    2  POINT (-3.07414 43.3205)
## 3 1.1026303       NA 1.081796 -0.02083402     NA    3    POINT (-2.3937 42.849)
## 4 1.3762752       NA 1.148880 -0.22739478     NA    4 POINT (-0.513889 38.3511)
## 5 1.3762752       NA 1.332646 -0.04362916     NA    5 POINT (-0.471944 38.3594)
## 6 1.1488804       NA 1.376275  0.22739478     NA    6 POINT (-0.506667 38.3403)

ok.cv <- krige.cv(formula = log10(AirPollutionLevel) ~ 1, locations = no2.sf, model = vmf)
head(ok.cv)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -8.2522 ymin: 38.3403 xmax: -0.471944 ymax: 43.4916
## Geodetic CRS:  ETRS89
##   var1.pred   var1.var observed    residual     zscore fold
## 1  1.102084 0.02678944 0.863423 -0.23866066 -1.4581393    1
## 2  1.048951 0.04555412 1.253703  0.20475240  0.9593237    2
## 3  1.194732 0.19248550 1.081796 -0.11293549 -0.2574136    3
## 4  1.357330 0.02464658 1.148880 -0.20844919 -1.3277671    4
## 5  1.241641 0.06419993 1.332646  0.09100518  0.3591690    5
## 6  1.183676 0.02427743 1.376275  0.19259935  1.2360997    6
##                    geometry
## 1   POINT (-8.2522 43.4916)
## 2  POINT (-3.07414 43.3205)
## 3    POINT (-2.3937 42.849)
## 4 POINT (-0.513889 38.3511)
## 5 POINT (-0.471944 38.3594)
## 6 POINT (-0.506667 38.3403)

sk.cv <- krige.cv(formula = log10(AirPollutionLevel) ~ 1, locations = no2.sf, model = vmf, beta=1)
head(sk.cv)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -8.2522 ymin: 38.3403 xmax: -0.471944 ymax: 43.4916
## Geodetic CRS:  ETRS89
##   var1.pred   var1.var observed    residual     zscore fold
## 1  1.103176 0.02678602 0.863423 -0.23975263 -1.4649044    1
## 2  1.049311 0.04555377 1.253703  0.20439242  0.9576408    2
## 3  1.196345 0.19247839 1.081796 -0.11454886 -0.2610957    3
## 4  1.356889 0.02464606 1.148880 -0.20800820 -1.3249720    4
## 5  1.242696 0.06419693 1.332646  0.08995017  0.3550135    5
## 6  1.184243 0.02427657 1.376275  0.19203266  1.2324843    6
##                    geometry
## 1   POINT (-8.2522 43.4916)
## 2  POINT (-3.07414 43.3205)
## 3    POINT (-2.3937 42.849)
## 4 POINT (-0.513889 38.3511)
## 5 POINT (-0.471944 38.3594)
## 6 POINT (-0.506667 38.3403)

Comparamos las predicciones con los datos reales, para los dos tipos de interpolación:

par(mfrow=c(2,2))
plot(var1.pred~observed, thiessen.cv, pch = 20, main="Thiessen", col="cadetblue3", reset = FALSE)
plot(var1.pred~observed, idw.cv, pch = 20, main="IDW", col="cadetblue3", reset = FALSE)
plot(var1.pred~observed, ok.cv, pch = 20, main="Kriging ordinario", col="cadetblue3",reset = FALSE)
plot(var1.pred~observed, sk.cv, pch = 20, main="Kriging simple", col="cadetblue3")

Da la sensación de que en el gráfico del Kriging Ordinario los puntos están ligeramente más agrupados que en el de IDW.

Podemos calcular el coeficiente de correlación entre valores reales y predicciones (1 equivale a relación perfecta):

cor.idw.cv <- cor(idw.cv$var1.pred, idw.cv$observed)

cor.thiessen.cv <- cor(thiessen.cv$var1.pred, thiessen.cv$observed)

cor.ok.cv <- cor(ok.cv$var1.pred, ok.cv$observed)

cor.sk.cv <- cor(sk.cv$var1.pred, sk.cv$observed)

Y también podemos calcular el error de la interpolación (la media de los residuos). Cuanto más pequeño, mejor:

mean.idw.resid <- mean(idw.cv$residual)

mean.thiessen.resid <- mean(thiessen.cv$residual)

mean.ok.resid <- mean(ok.cv$residual)

mean.sk.resid <- mean(sk.cv$residual)

Y también podemos calcular el error cuadrático medio (RMSE):

rmse.idw <- sqrt(sum(idw.cv$residual^2)/length(idw.cv$residual))

rmse.thiessen <- sqrt(sum(thiessen.cv$residual^2)/length(thiessen.cv$residual))

rmse.ok <- sqrt(sum(ok.cv$residual^2)/length(ok.cv$residual))

rmse.sk <- sqrt(sum(sk.cv$residual^2)/length(sk.cv$residual))

Volcamos toda esta iformación en una tabla y comparamos los dos métodos para ver cuál de los dos seleccionamos:

col1<-cbind(cor.thiessen.cv,cor.idw.cv,cor.ok.cv,cor.sk.cv)
col2<-cbind(mean.thiessen.resid,mean.idw.resid,mean.ok.resid,mean.sk.resid)
col3<-cbind(rmse.thiessen,rmse.idw,rmse.ok,rmse.sk)

tabla<-rbind(col1,col2,col3)
rownames(tabla) <- c("Correlación", "Media residuos", "RMSE")
colnames(tabla) <- c("Polígonos de Thiessen", "IDW", "Kriging Ordinario","Kriging Simple")

kbl(tabla) %>%
  kable_styling(c("striped","bordered"))%>%
  column_spec(1:5, bold = T)%>%
  add_header_above(c("DATOS VALIDACIÓN CRUZADA" = 5),background = "blue",color = "white")

DATOS VALIDACIÓN CRUZADA
	Polígonos de Thiessen	IDW	Kriging Ordinario	Kriging Simple
Correlación	0.6263509	0.6787736	0.6948987	0.6910236
Media residuos	-0.0289108	-0.0856011	-0.0252730	-0.0302377
RMSE	0.2850853	0.2562169	0.2401582	0.2415531

Los datos obtenidos muestran que el método Kriging Ordinario ofrece mejores predicciones que el resto de métodos utilizados. No sólo la correlación es mayor sino que, además, tanto la media de los resiudos como el error cuadrático medio son más bajos. Por lo que nos quedamos con la predicción que proporciona el método de Kriging Ordinario

ok.pred<-ggplot() + geom_stars(data = ok, aes(fill = var1.pred, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción OK log NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()

ok.var<-ggplot() + geom_stars(data = ok, aes(fill = var1.var, x = x, y = y)) + 
  geom_sf(data = st_cast(mapa.sf, "MULTILINESTRING")) +
  xlab("Longitud") + ylab("Latitud") +
  labs(title = "Predicción varianzas OK log NO2 (ug/m3)") +
  scale_fill_gradientn(colours = sf.colors(20), na.value=NA)+
  theme_void()


grid.arrange(ok.pred, ok.var, nrow = 1)

---