Paquetes Utilizados
Para trabajar con datos espaciales (y espacio-temporales), usamos el paquete gstat, que incluye funcionalidad para kriging, entre otras muchas cosas.
install.packages("dplyr")
install.packages("sp")
install.packages("gstat")
install.packages("ggplot2")
install.packages("scales")
install.packages("magrittr")
install.packages("gridExtra")
library(sp)
library(dplyr)
library(gstat)
library(ggplot2)
library(scales)
library(magrittr)
library(gridExtra)
Conjunto de Datos
Los datos que estamos utilizamos son el conjunto de datos del proyecto meuse, que viene con el paquete sp.
# packages for manipulation & visualization
suppressPackageStartupMessages({
library(dplyr) # for "glimpse"
library(ggplot2)
library(scales) # for "comma"
library(magrittr)
})
data(meuse)
glimpse(meuse) # This is like a transposed version of print(): columns run down the page, and data runs across.
Por ejemplo, podemos inspeccionar visualmente como varia el zinc en el dominio de interes asignando el tamaño del punto en forma proporcional a la concentracion:
meuse %>% as.data.frame %>%
ggplot(aes(x, y)) + geom_point(aes(size=zinc), color="blue", alpha=3/4) +
ggtitle("Zinc Concentration (ppm)") + coord_equal() + theme_bw()
Implementacion
En terminos generales, hay algunos pasos a seguir en kriging de un conjunto de datos. Suponiendo que los datos se almacenan en un data frame, entonces uno debe:
- Convertir el data frame en spatial data frame (SPDF).
- Ajustar un modelo variograma a los datos.
- Realizar Kriging los datos segun el variograma hallado.
1) Convertir a un SPDF
En este momento, los datos de meuse estan un data frame:
class(meuse)
str(meuse)
Para convertirlo en un SPDF, primero debemos especificar cuales son las columnas contiene las coordenadas de los datos. Esto se realiza utilizando la notacion de formula de R de la siguiente manera:
coordinates(meuse) <- ~ x + y
class(meuse)
str(meuse)
meuse@coords
Comentarios Iniciales
Aqui vemos que suceden un par de cosas. Primero, el data frame se convierte en un SPDF. La clase SPDF esta estructurada de manera que los datos ahora se distinguen claramente de las coordenadas; se representa como un objeto S4, donde sus datos / atributos se almacenan en diferentes “ranuras”. En particular, para los objetos SPDF, hay cinco ranuras diferentes:
- data
- coords.nrs
- Coords
- bbox
- y proj4string.
Se puede acceder a cada una de estas ranuras directamente mediante el operador @ o mediante las funciones de ayuda que simplifican la sintaxis para acceder a ellas, por ejemplo,
# access various slots of the SPDF
bbox(meuse)
meuse@bbox
coordinates(meuse) %>% glimpse
Y, por supuesto, ambas formas de acceder a las ranuras dan los mismos resultados:
identical( bbox(meuse), meuse@bbox )
identical( coordinates(meuse), meuse@coords )
Algunas veces ciertas funciones (como con ggplot2) requieren el uso de un data frame y no un SPDF. En tales casos, se pueden forzar manualmente los datos a un data frame para retener la informacion de coordenadas, en lugar de simplemente acceder a la ranura de datos:
meuse@data %>% glimpse
meuse %>% as.data.frame %>% glimpse
Una vez forzados, los datos se pueden canalizar (pipe) a la siguiente funcion que desee utilizar.
2) Elaborando un variograma
Para realizar kriging, primero se debe tener un modelo variograma, desde el cual se pueden interpolar los datos. Hay un par de pasos involucrados:
-
Calcular el variograma muestra. Esto se hace con la funcion variograma.
-
Ajustar un modelo al variograma muestra.
Por ejemplo, un variograma podria ajustarse con el siguiente codigo:
lzn.vgm <- variogram(log(zinc)~1, meuse) # calculates sample variogram values
lzn.fit <- fit.variogram(lzn.vgm, model=vgm(1, "Sph", 900, 1)) # fit model
El variograma muestra y ajuste se pueden graficar en un lienzo a la vez para ver que tan bien fue el ajuste:
plot(lzn.vgm, lzn.fit) # plot the sample values, along with the fit model
3) Realizando Kriging
Setup
Por definicion, la interpolacion es la estimacion de valores en puntos para los que no tenemos mediciones en funcion de los puntos para los que si tenemos mediciones. Entonces, necesitamos dos dominios espaciales: uno que tenga valores asociados con los puntos y otro para el que deseamos estimaciones. En este ejemplo, los dominios espaciales que utilizamos son los de “meuse” y “meuse.grid”:
# load spatial domain to interpolate over
data("meuse.grid")
# to compare, recall the bubble plot above; those points were what there were values for. this is much more sparse
plot1 <- meuse %>% as.data.frame %>%
ggplot(aes(x, y)) + geom_point(size=1) + coord_equal() +
ggtitle("Points with measurements")
# this is clearly gridded over the region of interest
plot2 <- meuse.grid %>% as.data.frame %>%
ggplot(aes(x, y)) + geom_point(size=1) + coord_equal() +
ggtitle("Points at which to estimate")
library(gridExtra)
grid.arrange(plot1, plot2, ncol = 2)
Calculos Preliminares
Una vez preparado lo anterior, estamos listos para realizar kriging a los datos. Esto se puede hacer con la funcion gstat:: krige, que generalmente toma cuatro argumentos:
-
El modelo formula.
-
Un SPDF del dominio espacial que tiene las mediciones.
-
Un SPDF del dominio espacial sobre el cual se aplicara kriging.
-
Un modelo variograma ajustado a los datos.
Tenga en cuenta que los argumentos dos y tres tienen que ser SPDF, es decir, no pueden ser data frame.
Ahora, el paso de kriging se puede realizar en una sola llamada de funcion:
coordinates(meuse.grid) <- ~ x + y # step 3 above
lzn.kriged <- krige(log(zinc) ~ 1, meuse, meuse.grid, model=lzn.fit)
Utilizando Kriging Ordinario
Estos resultados podrian visualizarse como un mapa de calor:
lzn.kriged %>% as.data.frame %>%
ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
scale_fill_gradient(low = "yellow", high="red") +
scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
theme_bw()
Contacto
- Envienos sus comentarios al correo
henrytorrespo@yahoo.com