Kriging Ordinario en R

Esta entrada está basada en los trabajos del módulo, Estadística Espacial de mi maestría UNIGIS.

En esta estrada desarrollaremos un Kriging Ordinario utilizando el conjunto de datos airquality.csv, el cual contiene datos acerca de contaminación del aire y variables relacionadas para el estado de California, Estados Unidos. El conjunto de datos está constituido por 37 variables y 452 registros medidos desde estaciones ubicadas en toda la extensión de dicho estado.Para facilitar la lectura del código desarrollado, este se presenta a modo de bloques considerando:

Los comentarios en verde,
El codigo de color azul y
El resultado en gris.

# carga de paquetes necesarios

> library(sp)

> library(rgdal)

> library(gstat)

# cargamos el conjunto de datos

> air.quality <- read.csv2("./datos/airquality.csv",

header=T,

dec=",",

na.string=NA)

# comprobamos la dimensión del conjunto de datos

> dim(air.quality)

[1] 452 37

# escalamos la variable de distancia a km.

> air.quality$OZDLYAV <- air.quality$OZDLYAV*1000

Luego, A partir de air.quiality se crea un data.frame de objetos espaciales y proyectamos sus coordenadas. Para esto usamos las funciones:

coordinates: con la cual indicamos las variables del data.frame que contienen las coordenadas geográficas
proj4string: con esta función fijamos el Sistema de Referencia de Coordenadas (SRC, CRS en inglés) a longlat con datum NAD83.
CRS: introducimos todos los parámetros necesarios para trabajar con el SRC correcto.
spTransform: crea el objeto aq a partir de los atrobutos en air.quality y los información del SRC proporcionado.

> coordinates(air.quality) <- ~LONGITUDE + LATITUDE

> proj4string(air.quality) <- CRS('+proj=longlat +datum=NAD83')

TA <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-

4000000 +datum=NAD83 +units=km +ellps=GRS80")

aq <- spTransform(air.quality, TA)

Luego, fijamos la extensión de análisis mediante el archivo counties.rds, el cual contiene el perfil del estado de California.

> cageo <- readRDS('./datos/counties.rds')

> ca <- spTransform(cageo, TA)

r <- raster(ca)

res(r) <- 10

g <- as(r, 'SpatialGrid')

A continuación creamos el variograma empírico usando como argumento formula OZDLYAV~1 y lo graficamos mediante la función plot.

# creamos el variograma empírico

> gs <- gstat(formula=OZDLYAV~1, locations=aq)

> v <- variogram(gs, width=20)

> plot(v)

A partir del variograma empírico, ajustamos el modelo del variograma mediante una aproximación exponencial utilizando la función fit.variogram sobre el variograma v. A continuación graficamos el modelo del variograma v y su ajuste fve

# ajustamos el modelo del variograma

> fve <- fit.variogram(v, vgm(85, "Exp", 75, 20))

> fve

model psill range

1 Nug 21.96600 0.00000

2 Exp 85.52957 72.31404

# graficamos el modelo ajustado

> plot(variogramLine(fve, 400), type="l", ylim=c(0,120))

> points(v[,2:3], pch=20, col="red")

Luego de realizado el ajuste, estamos en la posición de calcular el Kriging Ordinario para este conjunto de datos. Lo visualizaremos junto a su gráfico de varianza.

# Kriging Ordinario

> k <- gstat(formula=OZDLYAV~1, locations=aq, model=fve)

> kp <- predict(k, g)

> spplot(kp)