Introducción

El dióxido de nitrogeno (NO2) es un gas contaminante que se emite en los procesos de combustión que se llevan a cabo en relación con el tráfico (sobre todo vehículos automóviles, y en especial de motores diésel) y con el transporte en general, así como en instalaciones industriales de alta temperatura y de generación eléctrica. En ambiente urbano, generalmente más del 75% del NO2 en aire ambiente es aportado por el tráfico rodado. (Fuente:https://www.miteco.gob.es/es/calidad-y-evaluacion-ambiental/temas/atmosfera-y-calidad-del-aire/calidad-del-aire/salud/oxidos-nitrogeno.aspx)

Niveles elevados de dióxido de nitrógeno pueden irritar los pulmones y disminuir la función pulmonar, así como disminuir la resistencia a infecciones respiratorias. Y es que la irritación que provoca este contaminante se relaciona con un aumento de la mucosidad de las vías altas respiratorias, lo que puede hacer aumentar las infecciones respiratorias y reagudizar los síntomas de pacientes con enfermedades crónicas respiratorias, asmáticos y alérgicos. De hecho, recientes estudios científicos relacionan la exposición a NO2 con una mayor incidencia de bronquitis, especialmente en mayores e inmunodeprimidos, así como de bronquiolitis en niños. Otros estudios apuntan también a relacionar contaminación atmosférica con un bajo peso al nacer y mayor probabilidad de parto prematuro, por lo que las embarazadas son un colectivo de especial protección ante estos episodios.(Fuente:https://madridsalud.es/dioxido-de-nitrogeno-y-salud/#)

La presente tarea tiene como objetivo determinar en qué zonas de la Península Ibérica 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.

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 actualmente por lo que seleccionamos los valores registrados a día 03/04/2022.(Link para descargar los datos: https://discomap.eea.europa.eu/map/fme/AirQualityUTDExport.htm)

ES_NO2 <- read_csv("C:/Users/sanch/OneDrive/Escritorio/BIA/2021-2022/Datos_Espaciales_y_Espacio_Temporales/practica7/ES_NO2.csv")

Para simplificar la tarea trabajaremos únicamente con las comunidades autónomas que componen la península ibérica. Por lo que, con la ayuda de los nombres que toma la variable “network_name”, eliminamos los datos que se corresponden con Ceuta, Melilla, Canarias y Baleares

unique(ES_NO2$network_name)
##  [1] "CCAA Castilla y Le<f3>n" "CCAA Andaluc<ed>a"      
##  [3] "CCAA Catalu<f1>a"        "EMEP/VAG/CAMP"          
##  [5] "CCAA Murcia"             "CCAA Islas Canarias"    
##  [7] "Ayto Madrid"             "CCAA Com. Valenciana"   
##  [9] "CCAA Asturias"           "CCAA Navarra"           
## [11] "Ayto Zaragoza"           "CCAA Pa<ed>s Vasco"     
## [13] "CCAA Madrid"             "CCAA Islas Baleares"    
## [15] "CCAA Arag<f3>n"          "CCAA Cantabria"         
## [17] "CCAA Cast. la Mancha"    "CCAA Extremadura"       
## [19] "CCAA Galicia"            "CCAA La Rioja"          
## [21] "CA Ceuta"
ES_NO2<- ES_NO2%>%
  filter(!network_name %in% c("CA Ceuta","CCAA Islas Canarias",
                              "CCAA Islas Baleares","EMEP/VAG/CAMP" ))

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$value_numeric) == TRUE )

ES_NO2 <- ES_NO2 %>% distinct(station_code, .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$value_numeric, # 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 23, por lo que, con la aydua de nuevo de la función filter, eliminamos dichos valores.

ES_NO2<- ES_NO2%>%
  filter(!value_numeric > 23 )

boxplot(ES_NO2$value_numeric, # 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"))
network_countrycode network_localid network_name network_namespace network_timezone pollutant samplingpoint_localid samplingpoint_namespace samplingpoint_x samplingpoint_y coordsys station_code station_localid station_name station_namespace value_datetime_begin value_datetime_end value_datetime_inserted value_datetime_updated value_numeric value_validity value_verification station_altitude value_unit
ES NET_ES208A CCAA Castilla y Le<f3>n ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_24089010_8_8 ES.BDCA.AQD -5.5663889 42.57528 EPSG:4979 ES1988A STA_ES1988A LE<d3>N 4 ES.BDCA.AQD 2022-04-03 00:00:00 2022-04-03 01:00:00 2022-04-02 23:49:10 2022-04-03 02:48:44+01:00 8 1 3 814 ug/m3
ES NET_ES209A CCAA Catalu<f1>a ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_17079003_8_8 ES.BDCA.AQD 2.8165473 41.97639 EPSG:4979 ES1999A STA_ES1999A Girona (Escola de M<fa>sica) ES.BDCA.AQD 2022-04-03 14:00:00 2022-04-03 15:00:00 2022-04-03 13:50:47 2022-04-03 16:49:31+01:00 8 1 3 73 ug/m3
ES NET_ES201A CCAA Andaluc<ed>a ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_14021007_8_8 ES.BDCA.AQD -4.7623400 37.89261 EPSG:4979 ES1799A STA_ES1799A LEPANTO ES.BDCA.AQD 2022-04-02 21:00:00 2022-04-02 22:00:00 2022-04-02 20:51:09 2022-04-02 23:51:54+01:00 11 1 3 119 ug/m3
ES NET_ES209A CCAA Catalu<f1>a ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_08019004_8_8 ES.BDCA.AQD 2.2045226 41.40388 EPSG:4979 ES0691A STA_ES0691A Barcelona (el Poblenou) ES.BDCA.AQD 2022-04-02 21:00:00 2022-04-02 22:00:00 2022-04-02 20:49:56 2022-04-02 23:50:08+01:00 20 1 3 3 ug/m3
ES NET_ES214A CCAA Murcia ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_30015001_8_8 ES.BDCA.AQD -1.8686111 38.11472 EPSG:4979 ES1882A STA_ES1882A CARAVACA ES.BDCA.AQD 2022-04-02 21:00:00 2022-04-02 22:00:00 2022-04-02 20:50:16 2022-04-02 23:50:34+01:00 7 1 3 1 ug/m3
ES NET_ES131A Ayto Madrid ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_28079047_8_8 ES.BDCA.AQD -3.6866667 40.39806 EPSG:4979 ES1937A STA_ES1937A MENDEZ ALVARO ES.BDCA.AQD 2022-04-03 13:00:00 2022-04-03 14:00:00 2022-04-03 12:48:57 2022-04-03 15:44:27+01:00 7 1 3 609 ug/m3
ES NET_ES210A CCAA Com. Valenciana ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_46220003_8_8 ES.BDCA.AQD -0.2347222 39.66722 EPSG:4979 ES1185A STA_ES1185A SAGUNT PORT ES.BDCA.AQD 2022-04-03 13:00:00 2022-04-03 14:00:00 2022-04-03 12:48:57 2022-04-03 15:44:27+01:00 1 1 3 10 ug/m3
ES NET_ES214A CCAA Murcia ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_30016018_8_8 ES.BDCA.AQD -1.0647222 37.69361 EPSG:4979 ES1423A STA_ES1423A LA ALJORRA ES.BDCA.AQD 2022-04-03 13:00:00 2022-04-03 14:00:00 2022-04-03 12:48:57 2022-04-03 15:44:28+01:00 6 1 3 80 ug/m3
ES NET_ES131A Ayto Madrid ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_28079036_8_8 ES.BDCA.AQD -3.6452778 40.40806 EPSG:4979 ES1426A STA_ES1426A MORATALAZ ES.BDCA.AQD 2022-04-03 14:00:00 2022-04-03 15:00:00 2022-04-03 13:49:16 2022-04-03 16:48:52+01:00 6 1 3 685 ug/m3
ES NET_ES203A CCAA Asturias ES.BDCA.AQD http://dd.eionet.europa.eu/vocabulary/aq/timezone/UTC NO2 SP_33004020_8_8 ES.BDCA.AQD -5.8989000 43.55030 EPSG:4979 ES0879A STA_ES0879A LLARANES ES.BDCA.AQD 2022-04-03 14:00:00 2022-04-03 15:00:00 2022-04-03 13:49:16 2022-04-03 16:48:52+01:00 3 1 3 8 ug/m3

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$value_numeric,
     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 0 y el 5 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$value_numeric)
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).

knitr::include_graphics('imgcordenadas.png')
Área de cobertura del sistema de coordenadas de referencia EPSG:4258

Área de cobertura del sistema de coordenadas de referencia EPSG: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("samplingpoint_x", "samplingpoint_y"), 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, mapping = aes(col = VALORLOG))+
  theme_void()

También podemos ver los puntos en un mapa más interactivo

mapview(no2.sf$geometry)

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] 83028

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

dim(st_coordinates(no2.sf))
## [1] 408   2
(zi <- st_coordinates(no2.sf)[1,]) #coordenadas del primer punto
##         X         Y 
## -5.566389 42.575278
(zj <- st_coordinates(no2.sf)[2,])  #coordenadas del segundo punto
##         X         Y 
##  2.816547 41.976386
(sep <- dist(st_coordinates(no2.sf)[1:2,]))  #distancia entre ambos puntos
##          1
## 2 8.404302

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

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  2030  14.03169 0.4581156       0       0 var1
## 2  1575  44.43110 0.6298419       0       0 var1
## 3  1963  74.43728 0.6052445       0       0 var1
## 4  1629 102.32738 0.6576204       0       0 var1
## 5  1774 132.50682 0.7416319       0       0 var1
## 6  2069 162.15312 0.7787286       0       0 var1
## 7  2300 191.98751 0.7431689       0       0 var1
## 8  2618 220.97954 0.7386511       0       0 var1
## 9  2891 250.35718 0.7094686       0       0 var1
## 10 3441 280.05191 0.7273462       0       0 var1
## 11 4790 308.88654 0.7613817       0       0 var1
## 12 4234 337.65046 0.7206334       0       0 var1
## 13 3815 367.00214 0.6997176       0       0 var1
## 14 4159 396.83666 0.7781165       0       0 var1
## 15 3912 425.83711 0.7239984       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 25 metros, hasta una distancia máxima (cutoff) de 450 metros.

(ve <- variogram(VALORLOG ~ 1, no2.sf, cutoff = 450, width = 25))
##      np      dist     gamma dir.hor dir.ver   id
## 1  1790  12.29791 0.4502001       0       0 var1
## 2  1349  37.83720 0.5751358       0       0 var1
## 3  1431  62.97588 0.6423770       0       0 var1
## 4  1683  86.50930 0.5957699       0       0 var1
## 5  1374 112.53449 0.7326361       0       0 var1
## 6  1537 137.67641 0.7479800       0       0 var1
## 7  1760 162.75905 0.7569657       0       0 var1
## 8  1857 187.68832 0.7745962       0       0 var1
## 9  2201 212.39862 0.7024023       0       0 var1
## 10 2333 237.32969 0.7336474       0       0 var1
## 11 2631 262.56240 0.7398106       0       0 var1
## 12 3242 288.13588 0.7082900       0       0 var1
## 13 4168 312.45149 0.7749171       0       0 var1
## 14 3572 337.21074 0.7234920       0       0 var1
## 15 3306 361.99831 0.7061122       0       0 var1
## 16 3403 387.81432 0.7426960       0       0 var1
## 17 3571 412.54408 0.7587766       0       0 var1
## 18 3158 437.54490 0.7293715       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.1798

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.7215725
max(vc$dist)
## [1] 440.909
vc <- vc %>% arrange(dist)
head(vc)
##        dist      gamma dir.hor dir.ver   id left right
## 1 0.7215725 0.20202081       0       0 var1  144   106
## 2 0.8278407 0.00000000       0       0 var1  310   118
## 3 0.8656351 0.32880391       0       0 var1  236   183
## 4 0.8789724 0.14525861       0       0 var1   28    14
## 5 0.9307018 0.01662058       0       0 var1  251   109
## 6 0.9366611 0.04138049       0       0 var1  338   133

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.”

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.”

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.44,
          range=100,
          nugget=0.3)

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.4036002   0.0000
## 2   Sph 0.3281108 125.8406
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. A partir de una distancia de 100 metros no existe dependencia.

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, especialmente, en las grandes ciudades como Madrid, Barcelona, Granada, Huesca, Terrasa, LLeida, Ablacete, Cuenca, Murcia, San Sebastian y Bilbao. La Comunidad Valenciana presenta unos valores considerablemente bajos de este gas contaminante, aunque se aprecian valores algo más elevados en Valencia y Elche. En Castilla y León, destaca Valladolid como el área que más NO2 concentra dentro de esta comunidad; en Galicia destaca Lugo; en Asturias destaca Ourense; y en Navarra destaca Pamplona. Además llama la atención que en zonas como Jerez de la frontera y Talavera de la reina se registren valores tan altos de NO2.

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)") +
  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 Barcelona, área en la mayor es la concentración de NO2, además de Madrid, Granada, Oviedo, Valladolid, Bilbao, San Sebastián, Zaragoza y Huesca.

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$value_numeric/max(no2.sf$value_numeric))

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$value_numeric/max(no2.sf$value_numeric))+
  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. Se trata de Barcelona, área en la mayor es la concentración de NO2, además de Madrid, Granada, Oviedo, Valladolid, Bilbao, San Sebastián, Zaragoza y Huesca.

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(value_numeric) ~ 1, locations = no2.sf)
head(idw.cv)
## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -5.566389 ymin: 37.89261 xmax: 2.816547 ymax: 42.57528
## Geodetic CRS:  ETRS89
##   var1.pred var1.var observed    residual zscore fold
## 1 0.6041346       NA 0.903090  0.29895540     NA    1
## 2 0.8742550       NA 0.903090  0.02883503     NA    2
## 3 0.8586460       NA 1.041393  0.18274671     NA    3
## 4 1.1424883       NA 1.301030  0.15854165     NA    4
## 5 0.6696827       NA 0.845098  0.17541532     NA    5
## 6 1.0556766       NA 0.845098 -0.21057860     NA    6
##                     geometry
## 1 POINT (-5.566389 42.57528)
## 2  POINT (2.816547 41.97639)
## 3  POINT (-4.76234 37.89261)
## 4  POINT (2.204523 41.40388)
## 5 POINT (-1.868611 38.11472)
## 6 POINT (-3.686667 40.39806)
thiessen.cv <- krige.cv(formula = log10(value_numeric) ~ 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: -5.566389 ymin: 37.89261 xmax: 2.816547 ymax: 42.57528
## Geodetic CRS:  ETRS89
##   var1.pred var1.var observed    residual zscore fold
## 1 0.6020600       NA 0.903090  0.30103000     NA    1
## 2 1.0791812       NA 0.903090 -0.17609126     NA    2
## 3 1.0413927       NA 1.041393  0.00000000     NA    3
## 4 1.1760913       NA 1.301030  0.12493874     NA    4
## 5 0.7781513       NA 0.845098  0.06694679     NA    5
## 6 1.1461280       NA 0.845098 -0.30103000     NA    6
##                     geometry
## 1 POINT (-5.566389 42.57528)
## 2  POINT (2.816547 41.97639)
## 3  POINT (-4.76234 37.89261)
## 4  POINT (2.204523 41.40388)
## 5 POINT (-1.868611 38.11472)
## 6 POINT (-3.686667 40.39806)
ok.cv <- krige.cv(formula = log10(value_numeric) ~ 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: -5.566389 ymin: 37.89261 xmax: 2.816547 ymax: 42.57528
## Geodetic CRS:  ETRS89
##   var1.pred  var1.var observed   residual     zscore fold
## 1 0.4909367 0.5614679 0.903090  0.4121533  0.5500426    1
## 2 0.6335980 0.6569480 0.903090  0.2694920  0.3324914    2
## 3 0.6337037 0.5152038 1.041393  0.4076889  0.5679883    3
## 4 1.0919282 0.4533970 1.301030  0.2091018  0.3105407    4
## 5 0.7329661 0.6877665 0.845098  0.1121319  0.1352100    5
## 6 0.9676445 0.4349434 0.845098 -0.1225465 -0.1858166    6
##                     geometry
## 1 POINT (-5.566389 42.57528)
## 2  POINT (2.816547 41.97639)
## 3  POINT (-4.76234 37.89261)
## 4  POINT (2.204523 41.40388)
## 5 POINT (-1.868611 38.11472)
## 6 POINT (-3.686667 40.39806)
sk.cv <- krige.cv(formula = log10(value_numeric) ~ 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: -5.566389 ymin: 37.89261 xmax: 2.816547 ymax: 42.57528
## Geodetic CRS:  ETRS89
##   var1.pred  var1.var observed    residual      zscore fold
## 1 0.6028738 0.5606618 0.903090  0.30021622  0.40094397    1
## 2 0.8506969 0.6539086 0.903090  0.05239311  0.06479115    2
## 3 0.7111070 0.5148176 1.041393  0.33028572  0.46032347    3
## 4 1.1122301 0.4533703 1.301030  0.18879992  0.28039823    4
## 5 0.9724773 0.6840403 0.845098 -0.12737927 -0.15401319    5
## 6 0.9669872 0.4349433 0.845098 -0.12188918 -0.18481994    6
##                     geometry
## 1 POINT (-5.566389 42.57528)
## 2  POINT (2.816547 41.97639)
## 3  POINT (-4.76234 37.89261)
## 4  POINT (2.204523 41.40388)
## 5 POINT (-1.868611 38.11472)
## 6 POINT (-3.686667 40.39806)

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.3924133 0.4639246 0.4852647 0.4349050
Media residuos -0.0256372 -0.0528395 -0.0045039 -0.0639362
RMSE 0.4070866 0.3362296 0.3227619 0.3386253

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)

Conclusión

Gracias a los métodos de interpolación hemos obtenido información acerca los niveles de concentración de NO2 en el aire en el territorio penisnular partiendo de los valores recogidos en puntos concretos en los que se encuentran estaciones de medición y obtención de estos datos.

Hemos obtenido que las zonas con mayores concentraciones de NO2 son Barcelona y alrededores, Madrid y alrededores, Granada, Oviedo, Valladolid, Bilbao, San Sebastián, Zaragoza y Huesca.

“Zaragoza y Madrid, las que más empeoran por contaminación en NO2.”

“En lo que respecta a NO2, las grandes ciudades que más mejoraron sus niveles de contaminación fueron Valencia (con una bajada del 15%), Málaga (-6%) y Sevilla (-1%), mientras que en sentido contrario figuran Zaragoza y Madrid, con subidas del 7%.”

“Las urbes grandes intermedias que más mejoraron fueron Murcia (-38%), Alicante (-11%) y Palma de Mallorca (-7%), en tanto que el otro lado de la tabla lo lideran Vigo y Valladolid, con incrementos del 14% y del 15%, respectivamente.”

“En el caso de ciudades de entre 1.000 y 250.000 habitantes, las que más redujeron sus niveles de NO2 fueron Castellón de la Plana (25%), Cartagena (21%) y Logroño (15%). Leganés (9%), Salamanca (14%) y Badajoz (14%) fueron las que más los elevaron.”

“En el grupo de ciudades más pequeñas, las que más redujeron su contaminación por NO2 fueron Gandía (-27%), Palencia (-23%) y Avilés (-21%), y las que más aumentaron fueron Guadalajara (11%), Ciudad Real (31%) y Arrecife (43%).”

“En términos absolutos, las más contaminadas por NO2 fueron Leganés (32,5 μg/m3), Mollet del Vallès (30,6), Coslada (29,7), Madrid (29,2), Tarrasa (28,1), Getafe (28,1), Granollers (28,1), Granada (26,3), Alcalá de Henares (24,5) y Barcelona (24).” “Aquellas con mayor concentración de NO2 son las que conforman áreas urbanas de mayor tamaño, con gran población y en correspondencia con mayores parques de automóviles, como es el caso de Madrid, Barcelona, la costa de Málaga y Zaragoza.”

“En Granada, los altos niveles se deben a un régimen climático muy específico que, de forma parecida a Ourense, favorece la concentración de contaminantes al permanecer estables grandes masas de aire sobre sus áreas urbanas.”

El conocimiento de esta información es relevante a la hora de determinar posibles eventos futuros. De hecho, durante la cuarentena por el COVID-19 fueron varios científicos los que se dedicaron a predecir en que zonas iba a tener mayor impacto la presencia del virus. Esto era posible gracias a la utilización de datos como la edad de la población de la zona, casos de enfermedades similares en años anteriores y, además de otras variables, la calidad del aire. Como bien hemos explicado en la introducción de la tarea, una de las consecuencias para la salud por niveles elevados de dióxido de nitrogeno es la disminución de la resistencia a infecciones respiratorias (como por ejemplo el COVID-19, que se ha demostrado que puede afectar negativamente al sistema respiratorio).

En definitiva, esta información puede ser muy útil y relevante ya no solo en lo que se refiere a sanidad pública sino que, además, para otros aspectos como a la localización de una fábrica, en un estudio de mercado para determinar dónde abrir una tienda, la apertura de parques, etc.