Análisis espacio tiempo del NOx por medio de datos funcionales

library(fda)
library(fda.usc)
library(geoR)
library(splines)
library(colorspace)
library(GenSA)
library(gstat)
library(sp)
library(SpatFD)
library(data.table)
library(rgdal)
library(dplyr)
library(sf)
library(raster)
library(ggplot2)
library(stringr)
library(scatterplot3d)
library(rgeos)
library(spacetime)
library(lattice)
library(lubridate)
library(rstudioapi)
library(leaflet)
library(RColorBrewer)
library(leaflet.esri)

Preparación de los datos

La idea de usar datos funcionales en este caso indexados por el tiempo no se da como una necesidad por términos de capacidad de computación si no como un ejercicio para evaluar la capacidad de predicción de esta metodología. Para poder llevar a cabo la metodología, aunque no es necesario por la misma dinamica de los datos funcionales, se eliminan las estaciones que no cuentan con los 5 años de observaciones:

m <- getData(name = "GADM", country = "Spain", level = 0)

Warning in getData(name = "GADM", country = "Spain", level = 0): getData will be removed in a future version of raster
. Please use the geodata package instead

m <- m %>%
  st_as_sf() %>%
  st_cast("POLYGON") %>%
  mutate(area = st_area(.)) %>%
  arrange(desc(area)) %>%
  slice(1)

Warning in st_cast.sf(., "POLYGON"): repeating attributes for all
sub-geometries for which they may not be constant

m <- m %>% st_transform(25830)
ggplot(m) + geom_sf() + theme_bw() + coord_sf(datum = st_crs(m))

m <- as_Spatial(m[[4]])

airePM10 <- fread("./España/DataExtract.csv", header=T, dec=".", sep=",")
names(airePM10) <- gsub("\\s", "_", names(airePM10))
airePM10 <- airePM10[Year %in% c(2012,2013,2014,2015,2016)]
variables <- names(airePM10)[c(1,4,6,7,8,11,12,18,19)]
airePM10 <- airePM10[,variables, with=FALSE]
airePM10 <- airePM10[!is.na(Air_Pollution_Level),]
eliminarPM10 <- which(duplicated(airePM10,by=c("Year","Air_Quality_Station_EoI_Code")))
airePM10 <- airePM10[-eliminarPM10,]
airePM10Sp <- airePM10
coordinates(airePM10Sp) <-  ~Longitude + Latitude
proj4string(airePM10Sp) <- CRS("+proj=longlat")
airePM10Sp <- spTransform(airePM10Sp,CRS("+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
filtro <- gIntersects(m,airePM10Sp,byid = T)

Warning in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func): spgeom1 and spgeom2
have different proj4 strings

airePM10Sp <- airePM10Sp[which(filtro),]
airePM10 <- data.frame(airePM10Sp) %>% data.table()

airePM10 <- airePM10[order(Air_Quality_Station_EoI_Code,Year),]
filtro <- airePM10[,(A=length(Country)),by=c("Air_Quality_Station_EoI_Code")]
filtro <- filtro[V1==5,]$Air_Quality_Station_EoI_Code 
airePM10 <- airePM10[Air_Quality_Station_EoI_Code %in% filtro,]

Como ya se mencionó antes tenemos sólo 5 años observados para todas y cada una de las estaciones en nuestros datos. Por tanto, se pueden construir los datos funcionales indexados por el tiempo. Para esto y por como lo solicitan las funciones de los paquetes fda.usc y fda es necesario guardar los datos dentro de una matriz. En nuestro caso cada columna representa una estación de medición y cada fila representan las observaciones en los años 2012 al 2016

Construcción de los datos funcionales

mairePM10<- matrix(airePM10$Air_Pollution_Level,nrow=5,ncol=length(filtro),byrow = T, dimnames = list(rownames=c(2012,2013,2014,2015,2016),colnames=filtro))
datosf=fdata(mairePM10,argvals=1:nrow(mairePM10))

Warning in fdata(mairePM10, argvals = 1:nrow(mairePM10)): The transposed data
is returned

plot(datosf)

Una vez construidos los datos datos funcionales pasamos a hacer el ajuste de los mismos por medio de B-SPLINES. Se usaron 7 bases como polinomios de grado 3

nbasis <-10
yearange <- c(1,nrow(mairePM10))
lambda=0.000001
harmaccelLfd <- vec2Lfd(c(1,length(filtro)), yearange)
yearbasis_Bsplines <- create.bspline.basis(yearange,nbasis,norder=4)
plot(yearbasis_Bsplines)

PM10_fdPar_Bspline <- fdPar(fdobj=yearbasis_Bsplines,Lfdobj=harmaccelLfd,lambda)
PM10_fd_Bspline  <- smooth.basis(argvals=1:nrow(mairePM10),mairePM10,PM10_fdPar_Bspline)
par(mfrow=c(1,2))
plot(PM10_fd_Bspline)

[1] "done"

plot(datosf)

¿Cómo quedaron ajustadas las curvas?

PM10_fd_Bspl=PM10_fd_Bspline$fd
# plot(PM10_fd_Bspl)
# lines(PM10_fd_Bspl,col=rainbow(10),lwd=2,lty=1)
# filtro[c(50,180,225,290)]
par(mfrow=c(2,2))
plotfit.fd(mairePM10[,c(50,180,225,290)], argvals=1:nrow(mairePM10),PM10_fd_Bspl[c("ES1117A", "ES1672A", "ES1852A", "ES2014A")],lwd=1,ylab=" ")

Multiple plots:  Click in the plot to advance to the next plot

ACP para datos funcionales

Con el fin de poder llevar los datos a un escalar y poder desde ahí aplicar la metodología Kriging es necesario llevar a cabo el ACP para datos funcionales y así reducir las funciones a las mínimas necesarias para poder proyectar las curvas sobre estas y así calcular los scores:

PM10_FPC_bspl=pca.fd(PM10_fd_Bspl,centerfns=T,nharm=3)
par(mfrow=c(1,3))
plot.pca.fd(PM10_FPC_bspl)

puntaje=PM10_FPC_bspl$scores
colnames(puntaje)=c("f1","f2","f3")
rownames(puntaje)=filtro
puntajes=as.data.frame(puntaje)
coordinates(puntajes) <- airePM10[which(!duplicated(airePM10$Air_Quality_Station_EoI_Code)),.(Longitude,Latitude)]

Ajuste de los semivariogramas para cada uno de los scores del ACP funcional. De acá en adelante ya se pueden manejar como si fueran escalares pues las curvas están reducidas en valores unicos escalares:

f1.vgm = variogram(f1~1,puntajes,cutoff=200000)
# plot(f1.vgm)
f1.fit <- fit.variogram(f1.vgm,model=vgm(60,"Gau",25000,15))
plot(f1.vgm,f1.fit)

f2.vgm = variogram(f2~1,puntajes,cutoff=200000)
# plot(f2.vgm)
f2.fit <- fit.variogram(f2.vgm,model=vgm(30,"Wav",50000,15))
plot(f2.vgm,f2.fit)

f3.vgm = variogram(f3~1,puntajes,cutoff=200000)
# plot(f3.vgm)
f3.fit <- fit.variogram(f3.vgm,model=vgm(20,"Wav",50000,10))
plot(f3.vgm,f3.fit)

Paquete SpatFD

De acá en adelante y para poder realizar el ejercicio se usará el paquete SpatFD que ya tiene todas las funciones programadas para poder llevar a cabo el ejercicio completo. La siguiente función realiza el proceso de ajuste de las curvas para los datos funcionales y ademas lleva a cabo el ACP funcional y extrae cada uno de los scores solicitados, en este caso solicitamos 3.

SFD_PM10 <- SpatFD(mairePM10, coords =airePM10[which(!duplicated(airePM10$Air_Quality_Station_EoI_Code)),.(Longitude,Latitude)], basis = "Bsplines", nbasis = 9,norder=4, lambda = 0.00002, nharm=3)

Para los modelos del semivariograma se usarán los ajustados anteriormente:

modelos <- list(vgm(60,"Gau",25000,15),
                vgm(30,"Wav",50000,15),
                vgm(20,"Wav",50000,10))

Creamos de nuevo un objeto con los datos que se van a predecir:

set.seed(41684)
puntos <- sp::spsample(m, n = 5000, type = "regular")
coordinates(puntos) ~ x1 + x2

coordinates(puntos) ~ x1 + x2

colnames(puntos@coords) <- c("Longitude", "Latitude")
puntos <- spTransform(puntos,CRS("+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))

Warning: PROJ support is provided by the sf and terra packages among others

plot(m)
plot(puntos, main="Malla de puntos",add=T)

puntos <- data.frame(puntos)

Se lleva a cabo la predicción Kriging de cada uno de los puntos generados anteriormente. Esta predicción se hace bajo los scores de las funciones:

KS_SFD_PM10_sc <- KS_scores_lambdas(SFD_PM10, puntos ,method = "scores", model = modelos)

Using first variable by default

Using fill.all = TRUE by default

[using simple kriging]
[using simple kriging]
[using simple kriging]

plot(recons_fd(KS_SFD_PM10_sc), main="Reconstrucción de las curvas")

[1] "done"

Predicción para los tiempos 1.5, 2.5 y 3.5

ggmap_KS(KS_SFD_PM10_sc,
         map_path = map,
         window_time = c(1.5,2.5,3.5),
         method = "scores",
         zmin = 6,
         zmax = 120)

Using fill.all = TRUE by default

[using simple kriging]
[using simple kriging]
[using simple kriging]

[[1]]

[[2]]

[[3]]