Carga de librerías

library(spdep)

## Loading required package: sp

## Loading required package: Matrix

library(akima)
library(geoR)

## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.1 (built on 2015-04-15) is now loaded
## --------------------------------------------------------------

data(oldcol)
attach(COL.OLD)

summary(CRIME)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1783 20.0500 34.0000 35.1300 48.5900 68.8900

Test de hipótesis para la normalidad

par(mfrow = c(1,2))
hist(CRIME)
qqnorm(CRIME)
qqline(CRIME, col = 2)

shapiro.test(CRIME)

## 
##  Shapiro-Wilk normality test
## 
## data:  CRIME
## W = 0.9722, p-value = 0.2956

Test de hipótesis para la normalidad

par(mfrow = c(2,2))
hist(X)
qqnorm(X)
qqline(X, col = 2)
shapiro.test(X)

## 
##  Shapiro-Wilk normality test
## 
## data:  X
## W = 0.96967, p-value = 0.2349

hist(Y)
qqnorm(Y)
qqline(Y, col = 2)

shapiro.test(Y)

## 
##  Shapiro-Wilk normality test
## 
## data:  Y
## W = 0.96468, p-value = 0.1475

Interpolación

par(mfrow=c(1,2))
int.crime <- interp.new(Y, X, CRIME)
image(int.crime, ylab="X (latitud)", xlab="Y (longitud)", 
      main ="Interpolación para CRIME")
contour(int.crime, ylab="X (latitud)", xlab="Y (longitud)", 
        main = "Curvas de Nivel para CRIME")

Descriptivos geoestadísticos

obj <- cbind(X, Y, CRIME)
CRIME.geo <- as.geodata(obj, coords.col = 1:2, data.col = 3)
plot(CRIME.geo)

summary(CRIME.geo)

## Number of data points: 49 
## 
## Coordinates summary
##         X     Y
## min 24.25 24.96
## max 51.24 44.07
## 
## Distance summary
##        min        max 
##  0.7421561 27.0128174 
## 
## Data summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1783 20.0500 34.0000 35.1300 48.5900 68.8900

maxd = summary(CRIME.geo)$distances.summary[[2]] #distancia maxima

Variograma empírico

par(mfrow=c(1,2))
CRIME.var <- variog(CRIME.geo, uvec = seq(0, maxd/2, l = 12))

## variog: computing omnidirectional variogram

plot(CRIME.var)
CRIME.var.c <- variog(CRIME.geo, uvec = seq(0, maxd/2, l = 12), bin.cloud = T)

## variog: computing omnidirectional variogram

plot(CRIME.var.c,bin.cloud=T)

Variograma direccional

par(mfrow = c(1,1))
vario.4 <- variog4(CRIME.geo, max.dist = maxd/2)

## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram

plot(vario.4, omnidirectional = T, legend = F)

Modelos geoespaciales

CRIME.var.fit <- variofit(vario = CRIME.var, #variograma empírico
                          ini.cov.pars = c(375, 13), #valores escogidos al observar el var.empirico
                          cov.model = "matern") #modelo de variograma

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim

#c("matern", "exponential", "gaussian", "cauchy", "spherical", "cubic",
#"wave", "circular", "powered.exponential", "gneiting")
plot(CRIME.var, main="Variograma empírico y ajuste teórico")
lines.variomodel(cov.model = "matern",
                 cov.pars = CRIME.var.fit$cov.pars,
                 nug = CRIME.var.fit$nugget,
                 kappa = CRIME.var.fit$kappa,
                 max.dist = maxd/2,
                 col = 2)

#param_eye <- eyefit(CRIME.var)
#CRIME.var.fit <- variofit(vario = CRIME.var,
#                          ini.cov.pars = param_eye)

Predicción espacial

delta <- 5
min1 <- summary(CRIME.geo)$coords.summary[1,1]-delta
#El 20 es para que las mediciones del borde aparezcan al graficar
min2 <- summary(CRIME.geo)$coords.summary[1,2]-delta
max1 <- summary(CRIME.geo)$coords.summary[2,1]+delta
max2 <- summary(CRIME.geo)$coords.summary[2,2]+delta
loci <- expand.grid(seq(min1, max1, delta), seq(min2, max2, delta))
kconv <- krige.conv(geodata = CRIME.geo,
                    locations = loci,
                    krige = krige.control(obj.m = CRIME.var.fit))

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

contour(kconv,
        filled = TRUE,
        coords.data = CRIME.geo$coords,
        color = terrain.colors,
        main = "Estimación en el dominio usando Kriging")

contour(kconv,
        val=sqrt(kconv$krige.var),
        filled = TRUE,
        coords.data = CRIME.geo$coords,
        color = terrain.colors,
        main = "Estimación de la Varianza en el dominio usando Kriging")