Interpolación espacial y temporal

Maestría en Estadística Aplicada

Roberto Trespalacios

• Estadística espacio temporal
  • Modelo general
  • Función de covarianza espacio temporal
  • Covarianza basada en predicción espatio Temporal
  • Ejemplos de aplicación (ozono en vias trenes y velocidad del viento)

• Un aspecto importante del análisis estadístico es descomponer variación en componentes específicos.
• Esto generalmente puede ser realizado con procesos espacio-temporales también, obteniendo un componente general de variación del modelo para \( \{Z(s; t): s\in D, t \in T\} \) donde \[ Z(s; t) = \mu(s; \theta_t)+ \gamma(t; \theta_s)+ \kappa(s;t;\theta_{s;t}) +\delta(s; t) \] donde los índices de espacio y tiempo pueden representar continuo o procesos discretos y están en \( D_s \times T_t \subset \mathbb{R}^2 \times \mathbb{R} \).
• Por lo tanto, este modelo separa la variabilidad espacial, temporal y espacio-temporal en componentes aditivos, aquí, la generalidad sustancial surge de permitir que los parámetros varíen de espacio y/o tiempo también.

• La tendencia espacial en nuestro modelo está representado por \( \mu(s; \theta_t) \), que puede tener parámetros que varían con eltiempo y el espacio.
• Del mismo modo, \( \gamma(t; \theta_s) \) representa una tendencia temporal que puede tener parámetros espacialmente variables.
• Más complicada es la dependencia espacio-temporal (por ejemplo, como sugerido por un proceso dinámico) típicamente se representa por \( \kappa(s;t;\theta_{s;t}) \).
• El último término, \( \delta(s; t) \), representa variación(error) espaciotemporal a microescala que es a menudo modelado como un proceso sin correlación media cero en \( D_s \times T_t \subset \mathbb{R}^2 \times \mathbb{R} \)

• La esencia del modelado espacio-temporal descriptivo concierne a la función de covarianza asociada.
• Decir estamos interesados en la covarianza espacio-temporal función para el proceso \( Z(·) \), \( Cov(Z(s_1; t_1), Z(s_2; t_2)) \).
• No cualquier función puede ser una función de covarianza; un la función de covarianza es una función no negativa-definida (y viceversa), y típicamente es suficiente para el función para ser positivo.
• Por lo general, se consideran tales funciones que también son estacionarias.

• En particular, la estacionariedad de segundo orden a menudo se considera, en la cual la expectativa de \( Z(s; t) \) es constante en todo el dominio de interés y la covarianza es solo una función de la Lag espacial y temporal, que, en el espacio continuo (\( \mathbb{R}^d \)) y tiempo continuo (\( \mathbb{R} \)) viene dado por \[ Cov(Z(s_1; t_1), Z(s_2; t_2))=C(s_1-s_2;t_1-t_2) \] para \( s_1,t_1 \in \mathbb{R}^d \) y \( s_2,t_2 \in \mathbb{R} \)
• La función de correlación espacio-temporal estacionaria asociado con \( C (·; ·) \) es entonces \[ \rho(h;\tau) =C(h;\tau)/C(0;0); \qquad h\in \mathbb{R}^d, \tau \in \mathbb{R} \] donde \( h = s_1-s_2 \) y \( \tau = t_1 -t_2 \), representan espacial y retrasos temporales, respectivamente.
• Tenga en cuenta que uno puede tener espacialidad espacial pero no estacionariedad temporal, y viceversa.
• También es útil considerar la noción de isotropía espacial, que corresponde a \[ Cov(Z(s_1; t_1), Z(s_2; t_2))=C(||s_1-s_2||;t_1,t_2) \]

• Se dice que un proceso aleatorio \( Z(·; ·) \) tiene función de covarianza espacio-temporal separable si, para todo \( s_1,t_1 \in \mathbb{R}^d \) y \( s_2,t_2 \in \mathbb{R} \) \[ Cov(Z(s_1; t_1), Z(s_2; t_2)) = C^{(s)}(s_1, s_2) \cdot C^{(t)}(t_1, t_1) \]

• donde \( C^{(s)} \) y \( C^{(t)} \) son covarianzas espaciales y temporales funciones, respectivamente.

• Una consecuencia obvia de la condición de separabilidad (aunque poco realista) de funciones de covarianza espacio-temporales es dada por producto de la covarianza espacial y temporal individual funciones.

• El conocimiento de la función de covarianza puede permitir para realizar una predicción espacio-temporal análoga a 'Kriging' para procesos espaciales.
• Estamos interesados en predecir \( Z(s_0; t_0) \) dado un conjunto de datos \[ Z_Y =(Z(s_1,t_1), Z(s_2,t_2),\dots, Z(s_n,t_n))' \] que son observaciones tomadas en el espacio-tiempo conocido “Ubicaciones”.
• Supongamos que hay un proceso oculto (“verdadero”) de interés, \( \{Z(s; t): s \in D_s \subset \mathbb{R}^d, t \leqslant 0\} \), que no es observable debido a un error de medición y, posiblemente, “Faltante”. Escribamos \[ Z=Z_Y+\delta \] donde \( E(\delta) = 0 \), \( Cov(\delta) = \sigma^2_{\varepsilon}I \), y \( Z_Y = (Z(s_1; t_1),\dots,Z(s_n; t_n))' \).

• Buscamos predecir \( Z(s_0; t_0) \) con el predictor lineal, \( \lambda_0'Z + k \), donde \( \lambda_0 \), es un n-vector de coeficientes que deben ser estimados (por simplicidad, asumimos que \( E[Z_Y(s; t)] =0 \).
• Entonces, \( k = 0 \) y minimizamos con respecto a \( \lambda_0 \), el error de predicción cuadrático medio, \[ E[(Z_Y(s_0;t_0)-\lambda_0Z)^2] \]

• Podemos probar que el estimador de \( \lambda_0 \) es \( \hat{\lambda}_0C(s_0;t_0)C_Z^{-1} \) y por lo tanto, tenemos el simple espacio-temporal predictor de kriging: \[ \hat{Z}_Y =C(s_0;t_0)'C^{-1}_zZ \] donde \( C_Z =Cov(Z) \) y \( C(s_0; t_0)'= Cov(Z_Y(s_0; t_0), Z) \). El espacio-temporal asociado El error estándar de kriging simple (s.e.) es: \[ \sigma^2_{k}(s_0;t_0)=[Var(Z_Y(s_0,;t_0))-C(s_0;t_0)'C^{-1}_ZC(s_0,t_0)] \]

• Definimos el semivariograma espacio temporal empírico \( \hat{\gamma}(h_k,\tau_p) \), como promedio de incrementos al cuadrado para una clase \( (h_k,\tau_p) \), \[ \hat{\gamma}(h_k,\tau_p)=\frac{1}{2N(h_k,\tau_p)}\sum_{(s,t)\in N(h_K,\tau_p)}[Z(s_i,t_i)-Z(s_j,t_j)]^2 \]

Paquetes de para datos espacio-temporales

• La libreríagstat puede realizar kriging espacio-temporal explotando las funcionalidades del librería spacetime, que fue desarrollado por el mismo equipo que gstat.
• En spacetime tenemos dos formas de representar datos espacio-temporales: formatos STFDF y STIDF.
• El primero (STFDF) representa objetos con una cuadrícula completa de espacio tiempo. En otras palabras, en esta categoría se incluyen objetos como la grilla de estaciones meteorológicas.
• El objeto espacio-temporal se crea utilizando las \( n \) ubicaciones de las estaciones meteorológicas y los \( m \) intervalos de tiempo de sus observaciones. La cuadrícula espacio-temporal es de tamaño \( n \times m \).
• El segundo (STIDF) son los que vamos a usar para nuestro ejemplo de trabajo.
• Estos son objetos espacio-temporales no estructurados, donde tanto el espacio como el tiempo cambian dinámicamente.

• En este caso tenemos datos recopilados sobre tranvías que se mueven por la ciudad de Zurich. Esto significa que la ubicación de los sensores no es uniforme en toda la ventana de muestreo.

• La creación de objetos STIDF es bastante simple, solo tenemos que desensamblar el marco de datos que tenemos en un componente espacial, temporal y variables del problema, y luego fusionarlos para crear el objeto STIDF.

#librerias
library(gstat)
library(sp)
library(spacetime)
library(raster)
library(rgdal)
library(rgeos) 
library(readr)

• Hay un par de problemas por resolver antes de poder sumergirnos en la interpolación kriging.
• Lo primero que tenemos que hacer es, por supuesto, configurar el directorio de trabajo y cargar el archivo ozono.csv

ozono  <- read_csv("~/Desktop/curso_glm_utb/geo-estadistica/pres.geo/sec_6_geo/ozono.csv")

• Necesitamos transformar el objeto secundario en un objeto espacial, y luego cambié su proyección en UTM para que el variograma se calcule en metros y no en grados.
• Estos son los pasos necesarios para lograr todo esto:

#SpatialPointsDataFrame
coordinates(ozono)=~LON+LAT
projection(ozono)=CRS("+init=epsg:4326")

#Transformacion en proyeccion Mercator
ozone.UTM <- spTransform(ozono,CRS("+init=epsg:3395")) 

Ahora tenemos el objeto ozone.UTM, que es un SpatialPointsDataFrame con coordenadas en metros.

• Ahora, debemos crear el objeto SpatialPoints, con las ubicaciones de los sensores en un momento dado:

ozoneSP <- SpatialPoints(ozone.UTM@coords,CRS("+init=epsg:3395")) 

• Esto es simple de hacer con la función SpatialPoints en la librería sp.
• Esta función toma dos argumentos, el primero es una matriz o un data.frame con las coordenadas de cada punto.
• En este caso, utilicé las coordenadas de SpatialPointsDataFrame que creamos antes, que se proporcionan en un formato de matriz.

• Resulta que tenemos un par de duplicados, que debemos eliminar.
• Podemos hacer eso directamente en las dos líneas de código que necesitaríamos para crear los datos y el componente temporal para el objeto STIDF.

ozoneDF <- data.frame(PPB=ozone.UTM$ozone_ppb[-dupl[,2]]) 

• En esta línea final, se crea un data.frame con solo una columna, llamada PPB, con observaciones de ozono en parte por billón.

• Como puede ver, se eliminan los puntos duplicados al excluir las filas del objeto ozone.UTM con los índices incluidos en una de las columnas del objeto dupl.

• Podemos usar el mismo procedimiento para la parte temporal

ozoneTM <- as.POSIXct(ozone.UTM$TIME[-dupl[,2]],tz="CET") 

Combinación de los objetos: ozoneSP, ozoneDF y ozoneTM de clase STIDF:

timeDF <- STIDF(ozoneSP,ozoneTM,data=ozoneDF) 

• Este es el archivo que vamos a usar para calcular el variograma y realizar la interpolación espacio-temporal.
• Podemos verificar los datos brutos contenidos en el objeto STIDF utilizando la versión espacio-temporal de la función spplot, que es stplot:

stplot(timeDF) 

• Básicamente, la versión espacio-temporal del variograma incluye diferentes retardos temporales.
• Por lo tanto, con lo que terminamos no es un solo variograma sino una serie, que podemos trazar usando la siguiente línea

plot(var,map=TRUE) 

O también, podemos visualizarlo en tres dimensiones con el parámetro wireframe

plot(var,wireframe=TRUE) 

• Como se hizo en el el método kriging para dos dimensiones, en este punto debemos ajustar un modelo a nuestro variograma.
• Para hacerlo, utilizaremos la función vgmST y fit.StVariogram, que son las coincidencias espaciotemporales para vgm y fit.variogram.

• A continuación, vemos el código para algunos posibles modelos.

• Los modelos de variograma, en gstat tienen 5 opciones: separable(separable), suma de producto(productSum), métrica(metric), suma métrica(sumMetric) y suma simple(simpleSumMetric).

• Puede encontrar más información para adaptarse a este modelo, incluidas todas las ecuaciones que se presentan a continuación, en (Gräler et al., 2015).

Existen diferentes tipos de modelos separables. Presentamos el siguiente ejemplo:

\[ C(h,u) = \frac{\sigma^2}{(|h|^{2\lambda}+1)^\tau} \exp\Big\{ - \frac{c\|u\|^{2\lambda}}{(|h|^{2\lambda}+1)^{\beta\lambda}} \Big\} \]

donde \( \tau \) determina la suavidad de la correlación temporal, mientras que \( \lambda \) determina la suavidad de la correlación espacial, \( c \) determina la fuerza de la correlación espacial, y \( \beta \in (0,1) \) determina la fuerza de la interacción espacio/tiempo.

• Podemos comparar visualmente todos los modelos que instalamos utilizando wireframes de la siguiente manera:

plot(var,list(separable_Vgm, prodSumModel_Vgm, metric_Vgm, sumMetric_Vgm, SimplesumMetric_Vgm),all=T,wireframe=T) 

• El parámetro más importante a tener en cuenta para seleccionar el mejor modelo es sin duda el MSE.
• Al observar estos, está claro que el mejor modelo es la suma de la métrica, ya que este tiene el error más pequeño; este modelo lo usaremos para la interpolación kriging.

• Dado que estamos realizando una interpolación espacio-temporal, está claro quequeremos estimar nuevos valores tanto en el espacio como en el tiempo.
• Por esta razón, necesitamos crear una grilla de predicción espacio-temporal.
• En este caso, e debe descargar el archivo vectorial carreteras_zu.shp de la red de carreteras para el área alrededor de Zurich, luego la recorté para que coincidiera con la extensión de mi área de estudio, y luego creé la grilla espacial.

carreteras <- shapefile("~/Desktop/curso_glm_utb/geo-estadistica/pres.geo/sec_6_geo/datos_shp_vias_zu/vias_zu.shp") 

Luego cambiemos la proyección de este objeto de a UTM, de modo que sea comparable con lo que usé hasta ahora.

carr.UTM <- spTransform(carreteras,CRS("+init=epsg:3395")) 

Esto nos muestra la red de carreteras alrededor de las ubicaciones donde se recopilaron los datos.

Veamos el resultado con las siguientes lineas de código

plot(carr.UTM, cex=.2, pch=1)
plot(ozone.UTM,col="blue",add=T,cex=.5, pch=1) 

donde las carreteras carr.UTM están en negro y los puntos de datos están representados en rojo. Con esta selección ahora puedo usar la función spsample para crear una cuadrícula aleatoria de puntos a lo largo de las líneas de la carretera. Esto genera la siguiente grilla.

sp.grid.UTM <- spsample(carr.UTM,n=1500,type="random") 

• Como mencioné, ahora necesitamos agregar un componente temporal a esta grilla. Podemos hacer eso de nuevo usando el paquete espacio-tiempo. Primero tenemos que crear un vector(secuencia) de fecha/hora usando la función seq

tm.grid <- seq(as.POSIXct('2011-12-12 06:00 CET'),as.POSIXct('2011-12-14 09:00 CET'),length.out=5) 

• Esto crea un vector con 5 elementos (length.out = 5), con valores POSIXct entre las dos fechas/tiempos proporcionados.
• En este caso, estamos interesados en crear un dataframe espacio-temporales, ya que aún no tenemos datos para él.
• Por lo tanto, podemos usar la función STF para fusionar datos espaciales y temporales en una cuadrícula espacio-temporal, así

grid.ST <- STF(sp.grid.UTM,tm.grid) 

Esto se puede usar como datos nuevos en la función kriging.

Kriging

• Este es el paso más fácil en todo el proceso.
• Ahora que ya hemos creado el marco de datos espaciotemporales, calculamos el mejor modelo de variograma y creamos la grilla de predicción espacio-temporal.
• Todo lo que tenemos que hacer ahora es una simple llamada a la función krigeST y el variograma que mejor nos resultó(sumMetric_Vgm) para realizar la interpolación

pred_krig <- krigeST(PPB~1, data=timeDF, modelList=sumMetric_Vgm, newdata=grid.ST) 

Podemos graficar los resultados nuevamente usando la función stplot

stplot(pred_krig)