Descripción de los datos
Cokriging para las variables NO2, O3, y NOX. La variable de principal riesgo es ozono (O3), asà que se usan las otras dos como covariables espaciales. DÃa 2020/01/16 A las 17 horas.
datos <- read.csv("Air_polution_cdmx_2020_01_16_17h.csv")
datos <- datos[c("Estacion",
"X",
"Y",
"NO2",
"O3",
"NOX")]
pander::pander((datos))
| AJU |
482901 |
2117907 |
NA |
50 |
NA |
| AJM |
478188 |
2130946 |
5 |
51 |
8 |
| ATI |
473346 |
2164689 |
20 |
70 |
20 |
| CAM |
482180 |
2152665 |
23 |
83 |
24 |
| CCA |
481502 |
2136931 |
6 |
46 |
8 |
| CUA |
469366 |
2141275 |
9 |
56 |
10 |
| CUT |
479189 |
2180751 |
13 |
75 |
14 |
| FAC |
474444 |
2154232 |
31 |
71 |
35 |
| HGM |
484020 |
2146380 |
25 |
81 |
31 |
| IZT |
487647 |
2143367 |
21 |
61 |
23 |
| LLA |
495842 |
2164872 |
14 |
81 |
15 |
| LPR |
487650 |
2160000 |
NA |
64 |
NA |
| MER |
487445 |
2147815 |
23 |
79 |
26 |
| MGH |
478716 |
2145543 |
NA |
66 |
NA |
| MON |
510196 |
2151776 |
9 |
80 |
10 |
| MPA |
498809 |
2123036 |
1 |
50 |
NA |
| NEZ |
497038 |
2144394 |
8 |
62 |
9 |
| PED |
478557 |
2136817 |
6 |
56 |
6 |
| SAG |
496819 |
2159801 |
16 |
80 |
16 |
| SFE |
472393 |
2140390 |
9 |
66 |
10 |
| TAH |
498890 |
2128098 |
3 |
47 |
3 |
| TLA |
478535 |
2159383 |
26 |
64 |
30 |
| TLI |
481421 |
2167509 |
18 |
86 |
20 |
| UAX |
489113 |
2134517 |
NA |
54 |
NA |
| UIZ |
492241 |
2140751 |
7 |
45 |
7 |
| VIF |
489875 |
2173664 |
19 |
80 |
22 |
| XAL |
491355 |
2159031 |
27 |
80 |
27 |
Matrices de coregionalización.
Matriz definida positiva para el modelo Esférico.
mat1 <- cbind(c(30, 30, 30),
c(30, 50, 30),
c(30, 30, 35))
#matriz definida positiva "cercana"
mat1 <- data.frame(as.matrix(nearPD(mat1)$mat))
names(mat1) <- c("NO2", "O3", "NOX")
row.names(mat1) <- c("NO2", "O3", "NOX")
pander::pander(mat1)
| NO2 |
30 |
30 |
30 |
| O3 |
30 |
50 |
30 |
| NOX |
30 |
30 |
35 |
Matriz definida positiva para el modelo efecto Hueco.
mat2 <- cbind(c(13.02, 24.5, 18.739),
c(24.58, 46.4, 35.36),
c(18.73, 35.36, 26.95))
mat2 <- data.frame(as.matrix(nearPD(mat2)$mat))
names(mat2) <- c("NO2", "O3", "NOX")
row.names(mat2) <- c("NO2", "O3", "NOX")
pander::pander(mat2)
| NO2 |
13.02 |
24.54 |
18.73 |
| O3 |
24.54 |
46.4 |
35.36 |
| NOX |
18.73 |
35.36 |
26.96 |
Definición de objeto en gstat
Semivariogramas univariados
vgmno2 <- vgm(psill = mat1[1, 1],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 1],
model = "Hol",
range = 2294))
vgmo3 <- vgm(psill = mat1[2, 2],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 2],
model = "Hol",
range = 2294))
vgmnox <- vgm(psill = mat1[3, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[3, 3],
model = "Hol",
range = 2294))
Semivarogramas cruzados (Bivariados)
vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 2],
model = "Hol",
range = 2294))
vgmno2_nox <- vgm(psill = mat1[1, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 3],
model = "Hol",
range = 2294))
vgmno3_nox <- vgm(psill = mat1[2, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 3],
model = "Hol",
range = 2294))
gstat
remove_na <- function(frame, vari_) {
# Remove na from sp object
datos1 <- frame
bool <- !is.na(datos1@data[vari_])
datos1@data <- datos1@data[bool, ]
datos1@coords <- datos1@coords[bool, ]
return(datos1)
}
coordinates(datos) <- ~ X + Y
g_st <- gstat(NULL,
id = "NO2",
formula = NO2 ~ X + Y,
model = vgmno2,
data = remove_na(datos, "NO2"))
g_st <- gstat(g_st,
id = "O3",
formula = O3 ~ Y,
model = vgmo3,
data = remove_na(datos, "O3"))
g_st <- gstat(g_st,
id = "NOX",
formula = NOX ~ Y,
model = vgmnox,
data = remove_na(datos, "NOX"))
#Cruzados
g_st <- gstat(g_st,
id = c("NO2", "O3"),
model = vgmno2_o3)
g_st <- gstat(g_st,
id = c("NO2", "NOX"),
model = vgmno2_nox)
g_st <- gstat(g_st,
id = c("O3", "NOX"),
model = vgmno3_nox)
pander::pander(do.call(rbind, g_st$model)[, 1:3])
| NO2.1 |
Hol |
13.02 |
2294 |
| NO2.2 |
Sph |
30 |
6096 |
| O3.1 |
Hol |
46.4 |
2294 |
| O3.2 |
Sph |
50 |
6096 |
| NOX.1 |
Hol |
26.96 |
2294 |
| NOX.2 |
Sph |
35 |
6096 |
| NO2.O3.1 |
Hol |
24.54 |
2294 |
| NO2.O3.2 |
Sph |
30 |
6096 |
| NO2.NOX.1 |
Hol |
18.73 |
2294 |
| NO2.NOX.2 |
Sph |
30 |
6096 |
| O3.NOX.1 |
Hol |
35.36 |
2294 |
| O3.NOX.2 |
Sph |
30 |
6096 |
Estimación del semivariograma
plot(variogram(g_st),
model = g_st$model,
pl = T,
xlab = "Distancias",
ylab = "Semivarianza")
##Mapas de predicción de O3 con las covariables espaciales NO2 y NOX
prediction_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 100000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
geom_tile() +
scale_fill_viridis_c() +
theme_void()
return(plot)
}
variance_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
var <- paste(variable, ".var", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
geom_tile() +
scale_fill_viridis_c(option = "inferno",
direction = -1) +
theme_void()
return(plot)
}
cv_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
var <- paste(variable, ".var", sep = "")
aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
geom_tile() +
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
return(plot)
}
pl1 <- prediction_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\JUAN FRANCO\Documents\septiembre25\CursoEspacial2020_II\Co_kri_gstat\SP\mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl2 <- variance_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\JUAN FRANCO\Documents\septiembre25\CursoEspacial2020_II\Co_kri_gstat\SP\mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl3 <- cv_plot(g_st, "O3",
"SP/mpiosutm.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\JUAN FRANCO\Documents\septiembre25\CursoEspacial2020_II\Co_kri_gstat\SP\mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
ggplotly(pl1)