Discusión del Problema

En la actualidad es de gran importancia e interés mantener un registro y un control de los datos que corresponden a la calidad del aire en las grandes ciudades.
En este informe se pretende mostrar el análisis de datos georreferenciados tomados en la ciudad de Seúl, Corea del Sur. Nos enfocaremos en las mediciones realizadas y registradas en 25 estaciones el día 23 de marzo del 2017 a las 12:00 p.m. Las variables que tomaremos para este estudio corresponden al NO2 y a PM2.5.

Dióxido de Nitrogeno

En el caso del Dióxido de Nitrógeno (NO2) este se genera principalmente por el tráfico de vehículos y por las emisiones de diversas industrias. Además, una alta concentración de dicho compuesto se relaciona con un aumento en el material particulado. Es un contaminante que afecta la calidad de vida de las personas, en especial población vulnerable, al incidir sobre el sistema respiratorio. Ocasionando así, irritación en los pulmones, disminuir la resistencia a infecciones pulmonares, mayor riesgo de padecer bronquitis, etc. También se ha relacionado con un bajo peso en los recien nacidos y en aumentar el riesgo de un parto prematuro.

Partículas en suspensión (PM2.5)

Son partículas en el aire que son menores a una medida de 2.5 micras. En gran medida se usan ya que están relacionadas con el aumento de emisiones de la quema de combustibles tipo Diesel. Además, al ser considerablemente pequeñas (100 veces menor en tamaño a un cabello humano) tienen una gran capacidad de alojarse e impactar el sistema respiratorio. Pueden estar conformadas por polvo, cenizas, hollín, partículas metálicas, cemento, polen, etc. Tienen efectos en la función pulmonar y recientemente se relacionan con problemas de tipo cardiovascular. Producen latidos irregulares, asma grave, infartos de miocardio no mortales. Estas partículas también contribuyen de manera negativa al medio ambiente, pues pueden hacer que los lagos se vuelvan ácidos, una disminución de nutrientes en los suelos perjudicando así la agricultura en las regiones.

El aumento de este material particulado es provocado también por fenómenos meteorológicos. En Corea del Sur, se pueden ver aumentos de dicho contaminante por el denominado “polvo amarillo” proveniente de los desiertos de China y Mongolia, lo cual ocasiona que las mediciones de estas variables se vean considerablemente más grandes.

Con esto en mente, es evidente la necesidad de tener información precisa que ayude a observar la problemática en materia de calidad de aire y permita que los organismos correspondientes tomen decisiones acertadas para disminuir el impacto que tienen estas variables en la calidad del aire y por ende, en la salud de las personas en su ciudad.

Análisis Geoestadístico

setwd('C:/Users/AxRen/Documents/UNAL/Espacial') ##Carpeta de trabajo

Haremos uso de las siguientes librerías:

library(readxl)
library(readr)
library(rgeos)
library(geoR)
library(sf)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(maps)
library(tools)
library(dplyr)
library(fields)
library(akima)
library(sp)
library(rgdal)
library(viridis)
library(scatterplot3d)
library(plotly)
library(cowplot)
library(pander)
library(gstat)
library(tidyverse)

Leyendo la base de datos

ContaminacionSeul<-read_excel("Aire_seul.xlsx")
head(ContaminacionSeul)
## # A tibble: 6 x 11
##   `Measurement date`  `Station code` Address      Latitude Longitude   SO2   NO2
##   <dttm>                       <dbl> <chr>           <dbl>     <dbl> <dbl> <dbl>
## 1 2017-03-23 12:00:00            101 19, Jong-ro~     37.6      127. 0.006 0.038
## 2 2017-03-23 12:00:00            102 15, Deoksug~     37.6      127. 0.006 0.031
## 3 2017-03-23 12:00:00            103 136, Hannam~     37.5      127. 0.006 0.037
## 4 2017-03-23 12:00:00            104 215, Jinheu~     37.6      127. 0.005 0.025
## 5 2017-03-23 12:00:00            105 32, Segeomj~     37.6      127. 0.008 0.023
## 6 2017-03-23 12:00:00            106 10, Poeun-r~     37.6      127. 0.006 0.03 
## # ... with 4 more variables: O3 <dbl>, CO <dbl>, PM10 <dbl>, PM2.5 <dbl>
str(ContaminacionSeul)
## tibble [25 x 11] (S3: tbl_df/tbl/data.frame)
##  $ Measurement date: POSIXct[1:25], format: "2017-03-23 12:00:00" "2017-03-23 12:00:00" ...
##  $ Station code    : num [1:25] 101 102 103 104 105 106 107 108 109 110 ...
##  $ Address         : chr [1:25] "19, Jong-ro 35ga-gil, Jongno-gu, Seoul, Republic of Korea" "15, Deoksugung-gil, Jung-gu, Seoul, Republic of Korea" "136, Hannam-daero, Yongsan-gu, Seoul, Republic of Korea" "215, Jinheung-ro, Eunpyeong-gu, Seoul, Republic of Korea" ...
##  $ Latitude        : num [1:25] 37.6 37.6 37.5 37.6 37.6 ...
##  $ Longitude       : num [1:25] 127 127 127 127 127 ...
##  $ SO2             : num [1:25] 0.006 0.006 0.006 0.005 0.008 ...
##  $ NO2             : num [1:25] 0.038 0.031 0.037 0.025 0.023 0.03 0.032 0.036 0.031 0.033 ...
##  $ O3              : num [1:25] 0.048 0.054 0.036 0.049 0.046 0.052 0.043 0.054 0.044 0.044 ...
##  $ CO              : num [1:25] 0.7 0.6 0.5 0.6 0.5 0.6 0.6 0.7 0.5 0.5 ...
##  $ PM10            : num [1:25] 87 90 84 84 97 89 94 101 94 96 ...
##  $ PM2.5           : num [1:25] 56 51 51 42 41 69 62 64 43 36 ...

Nos quedamos con las columnas correspondientes a Latitud, Longitud, NO2 y PM2.5

ContaminacionSeul<-data.frame(ContaminacionSeul[,c(4,5,7,11)])
head(ContaminacionSeul)
##   Latitude Longitude   NO2 PM2.5
## 1 37.57202  127.0050 0.038    56
## 2 37.56426  126.9747 0.031    51
## 3 37.54003  127.0049 0.037    51
## 4 37.60982  126.9348 0.025    42
## 5 37.59374  126.9497 0.023    41
## 6 37.55558  126.9056 0.030    69

Observamos los puntos correspondientes a las coordenadas

deci_coord = SpatialPoints(cbind(ContaminacionSeul$Longitude,
                                 ContaminacionSeul$Latitude),
                           proj4string = CRS("+proj=longlat"))
class(deci_coord)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
plot(deci_coord, axes = TRUE, main = "Lat-Long Coordinates")

Transformamos las coordenadas a coordenadas en el plano XY

CRS_UTM_S = (c("+proj=utm +zone=52 +datum=WGS84 
                  +units=m +no_defs"))

utm_coord = spTransform(deci_coord, CRS(CRS_UTM_S))
utm_coord_df = as.data.frame(utm_coord)
plot(utm_coord, axes = TRUE, main = "UTM Coordinates", col = "red")

Con lo cual, nos quedamos con los datos

ContaminacionSeul1<-data.frame(ContaminacionSeul,utm_coord_df)
head(ContaminacionSeul1)
##   Latitude Longitude   NO2 PM2.5 coords.x1 coords.x2
## 1 37.57202  127.0050 0.038    56  323822.7   4160202
## 2 37.56426  126.9747 0.031    51  321125.3   4159399
## 3 37.54003  127.0049 0.037    51  323733.5   4156653
## 4 37.60982  126.9348 0.025    42  317718.5   4164531
## 5 37.59374  126.9497 0.023    41  318988.8   4162718
## 6 37.55558  126.9056 0.030    69  315002.2   4158569

Mapa para la variable NO2

Importamos el archivo shp que contiene la información del mapa de Seúl

shp_map <- read_sf("seoul_municipalities.shp")
shp_map1<-shp_map[5]

Creamos el mapa

map_variable_seul <- function(datos, variable, map) {
  
  plot1 <- ggplot() +
    geom_sf(data = map, size = 0.3) +
    geom_point(data = datos, aes_string(x = "Longitude", y = "Latitude",
                                        color = variable)) +
    scale_color_viridis_c(direction = -1) +
    labs(
      x = "",
      y = "",
      title = paste(variable, "map Seoul")
    ) +
    theme_void() +
    theme(
      plot.title = element_text(size = 14,
                                face = "bold",
                                colour = "black"),
      plot.margin = unit(c(1, 1, 1.5, 1.2), "cm"))
  
  return(plot1)
  
}

pl1 <- map_variable_seul(ContaminacionSeul,
                           "NO2",
                           shp_map1)

ggplotly(pl1)

Mapa para la variable PM2.5

pl2 <- map_variable_seul(ContaminacionSeul,
                         "PM2.5",
                         shp_map1)

ggplotly(pl2)

Gráficos para observar la tendencia

simple_scatter_plot <- function(datos, variable1, variable2) {
  
  plot1 <- ggplot(as.data.frame(datos),
                  aes_string(x = variable1, y = variable2,
                             color = as.factor(1))) +
    geom_point() +
    scale_colour_viridis_d() +
    labs(
      x = variable1,
      y = variable2
    ) +
    theme_light() +
    theme(legend.position = "none")
  
  return(plot1)
  
}

Para la variable NO2

pl1 <- simple_scatter_plot(ContaminacionSeul1,
                           "coords.x1",
                           "NO2")

pl2 <- simple_scatter_plot(ContaminacionSeul1,
                           "coords.x2",
                           "NO2")

cowplot::plot_grid(pl1, pl2)

Podemos observar claramente una tendencia que posiblemente puede ser lineal o cuadrática en el gráfico de NO2 vs. las coordenadas en Y. Exploraremos esto más adelante.

Para la variable PM2.5

pl3 <- simple_scatter_plot(ContaminacionSeul1,
                           "coords.x1",
                           "PM2.5")

pl4 <- simple_scatter_plot(ContaminacionSeul1,
                           "coords.x2",
                           "PM2.5")

cowplot::plot_grid(pl3, pl4)

Aquí observamos cierta tendencia de la variable PM2.5 respecto a ambas coordenadas. Exploraremos esto con modelos de regresión que nos permitan quitar dicha tendencia.

Ajustando Modelos Variable NO2

Modelo 1

model1NO2 <- lm(NO2 ~ coords.x1 + coords.x2, data =ContaminacionSeul1 )
pander::pander(summary(model1NO2))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.203 0.7723 8.032 5.521e-08
coords.x1 1.285e-07 1.584e-07 0.8112 0.4259
coords.x2 -1.493e-06 1.891e-07 -7.899 7.278e-08
Fitting linear model: NO2 ~ coords.x1 + coords.x2
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 0.005289 0.7474 0.7244

Modelo 2

model2NO2 <- lm(NO2 ~ coords.x2, data =ContaminacionSeul1)
pander::pander(summary(model2NO2))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.053 0.7441 8.134 3.219e-08
coords.x2 -1.447e-06 1.789e-07 -8.088 3.557e-08
Fitting linear model: NO2 ~ coords.x2
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 0.005249 0.7398 0.7285

Modelo 3

model3NO2 <- lm(NO2 ~ coords.x2+I(coords.x2^2), data =ContaminacionSeul1)
pander::pander(summary(model3NO2))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 750.4 432.2 1.736 0.09654
coords.x2 -0.0003594 0.0002078 -1.729 0.09781
I(coords.x2^2) 4.303e-11 2.499e-11 1.722 0.09909
Fitting linear model: NO2 ~ coords.x2 + I(coords.x2^2)
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 0.005038 0.7707 0.7499

Ajustando Modelos Variable PM2.5

Modelo 1

model1 <- lm(PM2.5 ~ coords.x1 + coords.x2, data =ContaminacionSeul1 )
pander::pander(summary(model1))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 6412 1577 4.066 0.0005131
coords.x1 -0.000494 0.0003235 -1.527 0.1409
coords.x2 -0.001491 0.000386 -3.862 0.0008445
Fitting linear model: PM2.5 ~ coords.x1 + coords.x2 Vemos que la coordenada en X no está siendo significativa en el modelo. Trabajaremos sin tenerla presente y más adelante se intentará ajustar de otra forma.
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 10.8 0.5096 0.465

Modelo 2

model2 <- lm(PM2.5 ~ coords.x2, data =ContaminacionSeul1) 
pander::pander(summary(model2))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 6990 1574 4.44 0.0001881
coords.x2 -0.001668 0.0003786 -4.405 0.000205
Fitting linear model: PM2.5 ~ coords.x2
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 11.11 0.4576 0.434

Modelo 3

model3 <- lm(PM2.5 ~ coords.x2+I(coords.x2^2), data =ContaminacionSeul1)
pander::pander(summary(model3))
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 714413 962493 0.7423 0.4658
coords.x2 -0.3418 0.4628 -0.7386 0.468
I(coords.x2^2) 4.089e-08 5.564e-08 0.735 0.4701
Fitting linear model: PM2.5 ~ coords.x2 + I(coords.x2^2)
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
25 11.22 0.4706 0.4225

En el caso de la variable NO2 nos quedamos con el modelo 3, en donde sus términos son significativos y tenemos \(R^2\) altos.
En el caso de la variable PM2.5 nos quedamos con el modelo 2, pues de igual forma tenemos que sus términos son significativos para el modelo y además tenemos \(R^2\) altos.

Además, observando los residuales podemos ver estacionariedad en estos, que era lo esperado al remover la tendencia.

geo_data <- data.frame(ContaminacionSeul1, 
                       residuales_NO2 = residuals(model3NO2),residuales_PM25=residuals(model2))

pl1 <- simple_scatter_plot(geo_data,
                           "coords.x1",
                           "residuales_NO2")

pl2 <- simple_scatter_plot(geo_data,
                           "coords.x2",
                           "residuales_NO2")

cowplot::plot_grid(pl1, pl2)

pl3 <- simple_scatter_plot(geo_data,
                           "coords.x1",
                           "residuales_PM25")

pl4 <- simple_scatter_plot(geo_data,
                           "coords.x2",
                           "residuales_PM25")

cowplot::plot_grid(pl3, pl4)

Aquí podemos ver que para la variable NO2 se ha eliminado la tendencia. Sin embargo, los residuales de la variable PM2.5 vs. la coordenada en X presentan un comportamiento en el que se puede apreciar la tendencia para dos grupos de datos. Se intentará modelar de esta forma más adelante y observaremos si proporciona una mejor estimación.

Detección de Anisotropía NO2

.

c<-as.geodata(geo_data)
str(c)
## List of 2
##  $ Latitude                  , Longitude                 : num [1:25, 1:2] 37.6 37.6 37.5 37.6 37.6 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:25] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:2] "Latitude" "Longitude"
##  $ data                      : num [1:25] 0.038 0.031 0.037 0.025 0.023 0.03 0.032 0.036 0.031 0.033 ...
##  - attr(*, "class")= chr "geodata"
df_reg=c
df_reg$data=model3NO2$residuals


vari_0=variog(df_reg, estimator.type = "modulus", dir=0)
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_45=variog(df_reg, estimator.type = "modulus", dir=pi/4)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_90=variog(df_reg, estimator.type = "modulus", dir=pi/2)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_135=variog(df_reg, estimator.type = "modulus", dir=3*pi/4)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
par(mfrow=c(1,2))
plot(vari_0,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*0)))
plot(vari_90,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*90)))

par(mfrow=c(1,2))
plot(vari_45,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*45)))

plot(vari_135,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*135)),xlab="h",ylab="Semivarianza direccional")

Detección de Anisotropía PM2.5

.

c2<-as.geodata(geo_data)
df_reg2=c2
df_reg2$data=model2$residuals


vari_0_2=variog(df_reg2, estimator.type = "modulus", dir=0)
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_45_2=variog(df_reg2, estimator.type = "modulus", dir=pi/4)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_90_2=variog(df_reg2, estimator.type = "modulus", dir=pi/2)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_135_2=variog(df_reg2, estimator.type = "modulus", dir=3*pi/4)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
par(mfrow=c(1,2))
plot(vari_0_2,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*0)))
plot(vari_90_2,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*90)))

par(mfrow=c(1,2))
plot(vari_45_2,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*45)))

plot(vari_135_2,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*135)),xlab="h",ylab="Semivarianza direccional")

Modelamiento del Variograma variable NO2

Creación de objeto gstat para visualizar el variograma muestral

coordinates(geo_data) <-  ~ coords.x1 + coords.x2

g_obj <- gstat(id = "NO2",
               formula = NO2 ~ coords.x2+I(coords.x2^2),
               data = geo_data)
var_s <- variogram(g_obj,cutoff=15000)
var_s <- var_s[-1, ]  
plot(var_s, pl = T)

variograma muestral valores iniciales modelo 1

vgm_model1 <- vgm("Gau",
                  psill = 4.0e-05,
                  range = 10000)
plot(var_s, vgm_model1, pl = T)

Variograma ajustado modelo 1

fitted_vgm1 <- fit.variogram(object = var_s,
                             model = vgm_model1)
plot(var_s, fitted_vgm1, pl = T)

variograma muestral valores iniciales modelo 2

vgm_model2 <- vgm("Wav",
                  psill = 2.3e-05,
                  range = 2500)
plot(var_s, vgm_model2, pl = T)

Variograma ajustado modelo 2

fitted_vgm2 <- fit.variogram(object = var_s,
                             model = vgm_model2)
plot(var_s, fitted_vgm2, pl = T)

variograma muestral valores iniciales modelo 3

vgm_model3 <- vgm("Exp",
                  psill = 4e-05,
                  range = 9000)
plot(var_s, vgm_model3, pl = T)

Variograma ajustado modelo 3

fitted_vgm3 <- fit.variogram(object = var_s,
                             model = vgm_model3)
plot(var_s, fitted_vgm3, pl = T)

Luego de ajustar diferentes modelos podemos observar todas las propuestas juntas

VGauLine<-variogramLine(fitted_vgm1,15000)
VWavLine<-variogramLine(fitted_vgm2,15000)
VExpLine<-variogramLine(fitted_vgm3,15000)

ggplot(mapping=aes(dist,gamma))+
  geom_point(data=var_s)+
  geom_line(data = VGauLine,aes(color="Gaussiano"))+
  geom_line(data = VWavLine,aes(color="Wave"))+
  geom_line(data = VExpLine,aes(color="Exponencial"))

Realizamos la respectiva comparación de modelos

var_s$fitted1 <- variogramLine(fitted_vgm1, dist_vector = var_s$dist)$gamma
var_s$fitted2 <- variogramLine(fitted_vgm2, dist_vector = var_s$dist)$gamma
var_s$fitted3 <- variogramLine(fitted_vgm3, dist_vector = var_s$dist)$gamma

cmr <- function(obs, fit) {
  return(mean((obs - fit)**2))
}

tabla <- data.frame(rbind(fitted_vgm1,
                          fitted_vgm2,
                          fitted_vgm3),
                    CMR =
                      c(cmr(var_s$gamma, var_s$fitted1),
                        cmr(var_s$gamma, var_s$fitted2),
                        cmr(var_s$gamma, var_s$fitted3)))


pander::pander(tabla[c(1:3,10)])
model psill range CMR
Gau 2.681e-05 2081 7.685e-11
Wav 2.441e-05 2500 6.752e-11
Exp 2.622e-05 1161 7.779e-11

Se concluye que el modelo 1 (Modelo Wave) es el que presenta menor CMR por lo tanto será una primera opción para realizar el mapa de predicciones.

Predicciones NO2 modelo Wave

g_obj <- gstat(id = "NO2",
               formula = NO2 ~ coords.x2+I(coords.x2^2),
               model = fitted_vgm2, # Modelo ajustado 2
               data = geo_data)


map <- readOGR("seoul_municipalities.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\AxRen\Documents\UNAL\Espacial\seoul_municipalities.shp", layer: "seoul_municipalities"
## with 25 features
## It has 6 fields
## Integer64 fields read as strings:  ESRI_PK
class(map)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
new <- sp::spsample(map, n = 50000, type = "regular")
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
utm_coord2 = spTransform(new, CRS(CRS_UTM_S))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))

predic <- predict(g_obj, newdata = utm_coord2)
## [using universal kriging]
plot1<-ggplot()+
  geom_sf(aes(col=NO2.pred),data=predic)+
  scale_color_carto_c(palette='Prism', na.value = '#111111')


plot1

Podemos observar que hacia el sur de la capital del país Sur Coreano es en donde se concentran mayores cantidades de NO2.

Varianza

plot12<-ggplot()+
  geom_sf(aes(col=NO2.var),data=predic)+
  scale_color_carto_c(palette = 6, direction = -1)


plot12

De igual forma, tenemos una varianza muy pequeña cerca a los puntos en donde se encuentran las estaciones de medición y una variabilidad mayor a medida que nos alejamos de ellas.

Coeficiente de Variación

prediction<-as.data.frame(predic)
aux <- abs(sqrt(prediction["NO2.var"]) / abs(prediction["NO2.pred"]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
predic<-st_as_sf(prediction)

plot13<-ggplot()+
  geom_sf(aes(col=cv),data=predic)+
  scale_color_carto_c(palette = 2, direction = -1)


plot13

Ahora, dado que en los puntos en el variograma se observan datos que pueden ser atípicos, realizamos una propuesta robusta usando el método de Cressie.

Semivariograma con Cressie para el modelo lineal

g_objm2 <- gstat(id = "NO2",
                 formula = NO2 ~ coords.x2,
                 data = geo_data)

var_sm2c0 <- variogram(g_objm2,cutoff=15000,cressie=TRUE)
var_sm2c0
##    np      dist        gamma dir.hor dir.ver  id
## 1   1  1664.240 8.652615e-10       0       0 NO2
## 2   5  2539.361 1.776900e-05       0       0 NO2
## 3  14  3635.185 4.247265e-05       0       0 NO2
## 4  14  4500.973 1.472806e-05       0       0 NO2
## 5  12  5508.812 1.841832e-05       0       0 NO2
## 6  16  6430.014 7.788009e-06       0       0 NO2
## 7  18  7486.294 4.815550e-05       0       0 NO2
## 8  22  8473.032 1.149517e-05       0       0 NO2
## 9  22  9556.700 2.245542e-05       0       0 NO2
## 10 14 10501.684 2.041345e-05       0       0 NO2
## 11 20 11485.576 2.056315e-05       0       0 NO2
## 12 25 12605.682 3.175828e-05       0       0 NO2
## 13 17 13454.064 4.081152e-05       0       0 NO2
## 14 15 14548.032 2.166546e-05       0       0 NO2
var_sm2c <- var_sm2c0[-1, ]  
plot(var_sm2c, pl = T)

Aquí podemos ver que la diferencia no es mucha salvo para una corrección que hace para el dato 14 y 18 y también en algunos puntos inferiores. Realizaremos diferentes ajustes para observar cuál puede ser el más acertado

VGauLine<-variogramLine(fitted_vgm1m2c,15000)
VWavLine<-variogramLine(fitted_vgm2m2c,15000)
VExpLine<-variogramLine(fitted_vgm3m2c,15000)
VWav2Line<-variogramLine(fitted_vgm4m2c,15000)
VPenLine<-variogramLine(fitted_vgm5m2c,15000)
VPerLine<-variogramLine(fitted_vgm6m2c,15000)
VHolLine<-variogramLine(fitted_vgm7m2c,15000)
VSteLine<-variogramLine(fitted_vgm8m2c,15000)


ggplot(mapping=aes(dist,gamma))+
  geom_point(data=var_sm2c)+
  geom_line(data = VGauLine,aes(color="Gaussiano"))+
  geom_line(data = VWavLine,aes(color="Wave"))+
  geom_line(data = VExpLine,aes(color="Exponencial"))+
  geom_line(data = VWav2Line,aes(color="Wave2"))+
  geom_line(data = VPenLine,aes(color="Pentaspherical"))+
  geom_line(data = VPerLine,aes(color="Periodic"))+
  geom_line(data = VHolLine,aes(color="Hole"))+
  geom_line(data = VSteLine,aes(color="Ste Mat"))

Mostramos todas las propuestas juntas y miremos el CMR para los distintos modelos

var_sm2c$fitted1 <- variogramLine(fitted_vgm1m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted2 <- variogramLine(fitted_vgm2m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted3 <- variogramLine(fitted_vgm3m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted4 <- variogramLine(fitted_vgm4m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted5 <- variogramLine(fitted_vgm5m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted6 <- variogramLine(fitted_vgm6m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted7 <- variogramLine(fitted_vgm7m2c, dist_vector = var_sm2c$dist)$gamma
var_sm2c$fitted8 <- variogramLine(fitted_vgm8m2c, dist_vector = var_sm2c$dist)$gamma

cmr <- function(obs, fit) {
  return(mean((obs - fit)**2))
}

tablam2 <- data.frame(rbind(fitted_vgm1m2c,
                            fitted_vgm2m2c,
                            fitted_vgm3m2c,
                            fitted_vgm4m2c,
                            fitted_vgm5m2c,
                            fitted_vgm6m2c[2,],
                            fitted_vgm7m2c,
                            fitted_vgm8m2c[2,]),
                      CMR =
                        c(cmr(var_sm2c$gamma, var_sm2c$fitted1),
                          cmr(var_sm2c$gamma, var_sm2c$fitted2),
                          cmr(var_sm2c$gamma, var_sm2c$fitted3),
                          cmr(var_sm2c$gamma, var_sm2c$fitted4),
                          cmr(var_sm2c$gamma, var_sm2c$fitted5) ,
                          cmr(var_sm2c$gamma, var_sm2c$fitted6),
                          cmr(var_sm2c$gamma, var_sm2c$fitted7),
                          cmr(var_sm2c$gamma, var_sm2c$fitted8)))

ordm2c<-tablam2 %>%
  select(model,psill,range,kappa,CMR) %>%
  arrange(desc(CMR))

pander::pander(ordm2c)
  model psill range kappa CMR
22 Ste 0 5000 0.5 1.721e-10
21 Per 3.542e-06 7146 0.5 1.585e-10
11 Hol 2.489e-05 1221 0.5 1.568e-10
3 Exp 2.496e-05 510.6 0.5 1.441e-10
1 Gau 2.586e-05 1900 0.5 1.435e-10
5 Pen 2.561e-05 4161 0.5 1.428e-10
2 Wav 2.425e-05 2000 0.5 1.409e-10
4 Wav 2.36e-05 2500 0.5 1.289e-10

Observando valores más bajos para los modelos de onda y procedemos a realizar el mapa de predicciones con dicho modelo.

Mapa de Predicciones variable NO2 con Cressie

g_obj_m2 <- gstat(id = "NO2",
               formula = NO2 ~ coords.x2,
               model = fitted_vgm4m2c, # Modelo ajustado 1
               data = geo_data)
predicm2c <- predict(g_obj_m2, newdata = utm_coord2)
## [using universal kriging]
plot1m2c<-ggplot()+
  geom_sf(aes(col=NO2.pred),data=predicm2c)+
  scale_color_carto_c(palette='Prism', na.value = '#111111')

plot1m2c

Varianza

plot12m2c<-ggplot()+
  geom_sf(aes(col=NO2.var),data=predicm2c)+
  scale_color_carto_c(palette = 6, direction = -1)


plot12m2c

De igual forma que en el caso anterior se observa claramente que obtenemos varianzas pequeñas alrededor de las estaciones de medición.

Variograma para PM2.5

.

geo_dataPM <- data.frame(ContaminacionSeul1, residuales_PM = residuals(model2))

coordinates(geo_dataPM) <-  ~ coords.x1 + coords.x2

g_objPM <- gstat(id = "PM2.5",
               formula = PM2.5 ~ coords.x2,
               data = geo_dataPM)

var_sPM <- variogram(g_objPM,cutoff=15000)
var_sPM <- var_sPM[-1, ]  
plot(var_sPM, pl = T)

variograma muestral valores iniciales modelo 1 Matern

vgm_model1PM <- vgm("Mat",
                  psill = 110,
                  range = 3000)
plot(var_sPM, vgm_model1PM, pl = T)

Variograma ajustado modelo 1 Matern

fitted_vgm1PM <- fit.variogram(object = var_sPM,
                             model = vgm_model1PM)

plot(var_sPM, fitted_vgm1PM, pl = T)

variograma muestral valores iniciales modelo 2 Onda

vgm_model2PM <- vgm("Wav",
                  psill = 130,
                  range = 5000)
plot(var_sPM, vgm_model2PM, pl = T)

Variograma ajustado modelo 2 Onda

fitted_vgm2PM <- fit.variogram(object = var_sPM,
                             model = vgm_model2PM)
plot(var_sPM, fitted_vgm2PM, pl = T)

variograma muestral valores iniciales modelo 3 Exponencial

vgm_model3PM <- vgm("Exp",
                  psill = 150,
                  range = 9000)
plot(var_sPM, vgm_model3PM, pl = T)

Variograma ajustado modelo 3 Exponencial

fitted_vgm3PM <- fit.variogram(object = var_sPM,
                             model = vgm_model3PM)
plot(var_sPM, fitted_vgm3PM, pl = T)

variograma muestral valores iniciales modelo 4 Onda (diferentes parametros iniciales)

vgm_model4PM <- vgm("Wav",
                    psill = 110,
                    range = 3500)
plot(var_sPM, vgm_model4PM, pl = T)

Variograma ajustado modelo 4 Onda

fitted_vgm4PM <- fit.variogram(object = var_sPM,
                               model = vgm_model4PM)
plot(var_sPM, fitted_vgm4PM, pl = T)

Podemos verlos todos juntos

VMatLinePM<-variogramLine(fitted_vgm1PM,15000)
VWavLinePM<-variogramLine(fitted_vgm2PM,15000)
VExpLinePM<-variogramLine(fitted_vgm3PM,15000)
VWav2LinePM<-variogramLine(fitted_vgm4PM,15000)


ggplot(mapping=aes(dist,gamma))+
  geom_point(data=var_sPM)+
  geom_line(data = VMatLinePM,aes(color="Matern"))+
  geom_line(data = VWavLinePM,aes(color="Wave"))+
  geom_line(data = VExpLinePM,aes(color="Exponencial"))+
  geom_line(data = VWav2LinePM,aes(color="Wave2"))

var_sPM$fitted1 <- variogramLine(fitted_vgm1PM, dist_vector = var_sPM$dist)$gamma
var_sPM$fitted2 <- variogramLine(fitted_vgm2PM, dist_vector = var_sPM$dist)$gamma
var_sPM$fitted3 <- variogramLine(fitted_vgm3PM, dist_vector = var_sPM$dist)$gamma
var_sPM$fitted4 <- variogramLine(fitted_vgm4PM, dist_vector = var_sPM$dist)$gamma


cmr <- function(obs, fit) {
  return(mean((obs - fit)**2))
}

tabla <- data.frame(rbind(fitted_vgm1PM,
                          fitted_vgm2PM,
                          fitted_vgm3PM,
                          fitted_vgm4PM),
                    CMR =
                      c(cmr(var_sPM$gamma, var_sPM$fitted1),
                        cmr(var_sPM$gamma, var_sPM$fitted2),
                        cmr(var_sPM$gamma, var_sPM$fitted3),
                        cmr(var_sPM$gamma, var_sPM$fitted4)))


pander::pander(tabla[c(1:3,10)])
model psill range CMR
Mat 110.9 962.5 1094
Wav 118.1 5000 1645
Exp 110.9 966 1094
Wav 110 3500 1075

Caso en que se elige el modelo 4 de Onda. Vemos que el CMR más bajo es es que usa el semivariograma teórico de Onda

Predicciones de PM2.5 con el modelo wave2

g_objPM <- gstat(id = "PM2.5",
                 formula = PM2.5 ~ coords.x2,
                 model = fitted_vgm4PM, 
                 data = geo_dataPM)

map <- readOGR("seoul_municipalities.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\AxRen\Documents\UNAL\Espacial\seoul_municipalities.shp", layer: "seoul_municipalities"
## with 25 features
## It has 6 fields
## Integer64 fields read as strings:  ESRI_PK
new <- sp::spsample(map, n = 50000, type = "regular")
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
utm_coord2 = spTransform(new, CRS(CRS_UTM_S))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))
predic <- predict(g_objPM, newdata = utm_coord2)
## [using universal kriging]
plot1w<-ggplot()+
  geom_sf(aes(col=PM2.5.pred),data=predic)+
  scale_color_carto_c(palette='Prism', na.value = '#111111')
plot1w

De igual forma que en el caso anterior podemos ver como hay concentraciones más altas del material particulado PM2.5 hacia el sur de la ciudad. Además, se aprecia una especie de circunferencia con valores altos (círculo amarillo) hacia la zona de la frontera que comparten Yeong-deungpo-gu y Yang-cheong-gu.

Varianza

.

plot12w<-ggplot()+
  geom_sf(aes(col=PM2.5.var),data=predic)+
  scale_color_carto_c(palette = 3, direction = -1)
plot12w

Coeficiente de variación

prediction<-as.data.frame(predic)
aux <- abs(sqrt(prediction["PM2.5.var"]) / abs(prediction["PM2.5.pred"]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
predic<-st_as_sf(prediction)

plot13w<-ggplot()+
  geom_sf(aes(col=cv),data=predic)+
  scale_color_carto_c(palette = 3, direction = -1)


plot13w

De igual manera vemos una variabilidad baja en torno a las estaciones de medición y una variabilidad no muy alta en zonas que se alejan(menos del 20%). Sin embargo, hacia los extremos podemos ver una variabilidad cercana al 40%.

Referencias