Elias Medeiros
2 de fevereiro de 2017
Carregando os pacotes
require(gstat)
require(rgdal)
require(sp)
require(chron)
require(spacetime)
require(zoo)
require(xts)
require(maptools)
require(googleVis)
library(plotKML)
require(plotGoogleMaps)
library(RColorBrewer)Lendo dados do gstat
## [1] "STSDF"
## attr(,"package")
## [1] "spacetime"Visualizando as coordenadas
Visualização dos dados via gráfico de video
#dd2005 = data.frame(DE_RB_2005)
#dados <- data.frame("time" = as.Date(dd2005$time),
# "posto" = dd2005$station_european_code,
# dd2005[,1:2], "pm10" = dd2005$PM10)
#pgvis <- gvisMotionChart(dados, idvar='posto', timevar = 'time',
# options=list(showSidePanel=T))#op <- options(gvis.plot.tag='chart')
#plot(pgvis)Registro de dias observados em cada uma das 69 estações
barplot(sort(table(DE_RB_2005@index[,1])), col = heat.colors(69),
main="", ylab="número de dias", xaxt="n")
Lendo as bordas da região em estudo
data(air)
DE_NUTS1 <- spTransform(DE_NUTS1, CRS("+init=epsg:32632"))
class(DE_NUTS1)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"plot(DE_NUTS1)
Selecionando alguns dias para serem esboçados graficamente
set.seed(123)
smplDays <- sort(sample(365,8))
smplDays## [1] 17 105 149 190 287 320 340 362stplot(as(DE_RB_2005[,smplDays],"STFDF"),
col.regions=bpy.colors(120)[-(1:20)],
sp.layout = list("sp.polygons", DE_NUTS1), scales=list(draw=F),
key.space="right", colorkey=T, cuts=0:70,
main=NULL)
Variograma espaço-temporal (amostral)
\[\begin{eqnarray} \gamma(h, \mu) = \frac{1}{2}E\left[\varepsilon(s,t) - \varepsilon(s+h,t+\mu) \right]^{2} \label{var1} \end{eqnarray}\]#empVgm <- variogramST(PM10~1, DE_RB_2005, tlags=0:6)empVgm$dist <- empVgm$dist/1000
empVgm$avgDist <- empVgm$avgDist/1000
empVgm$spacelag <- empVgm$spacelag/1000plot(empVgm, wireframe=F, scales=list(arrows=F))
plot(empVgm, wireframe=T, scales=list(arrows=F))
Variograma puramente espacial
spEmpVgm <- empVgm[empVgm$timelag == 0,]
class(spEmpVgm) <- c("gstatVariogram", "data.frame")
spEmpVgm <- spEmpVgm[-1,1:3]
spEmpVgm$dir.hor <- 0
spEmpVgm$dir.ver <- 0spVgmMod <- fit.variogram(spEmpVgm, vgm(psill=80, model="Exp",
range=300, nugget=20))plot(spEmpVgm, spVgmMod)
Variograma puramente temporal
tmpEmpVgm <- empVgm[empVgm$spacelag == 0,]
class(tmpEmpVgm) <- c("gstatVariogram","data.frame")
tmpEmpVgm <- tmpEmpVgm[-1,c("np","timelag","gamma")]
colnames(tmpEmpVgm) <- c("np", "dist", "gamma")
tmpEmpVgm$dir.hor <- 0
tmpEmpVgm$dir.ver <- 0tmpVgmMod <- fit.variogram(tmpEmpVgm, vgm(psill=120, model="Sph",
range=6, nugget=30))plot(tmpEmpVgm, tmpVgmMod)
Estimando a anisotropia
spVgmMod$range## [1] 0.0000 224.3012linStAni <- estiStAni(empVgm, interval=spVgmMod$range)
linStAni ## [1] 117.3222Ajustes de diferentes modelos de variograma espaço-tempo
Modelo separable
\[ C_{sep}\left(h,\mu\right) = C_{s}\left(h\right)\times C_{t}\left(\mu\right) \] seu variograma é dado por:
\[ \gamma_{sep}\left(h,\mu\right) = sill\times\left(\bar{\gamma}_{s}(h) + \bar{\gamma}_{t}(\mu) - \bar{\gamma}_{s}(h)\bar{\gamma}_{t}(\mu)\right) \]
separableModel <- vgmST("separable",
space=vgm(psill=0.9,"Exp", range=200, nugget=0.1),
time =vgm(0.9,"Sph", range=3.5, nugget=0.1),
sill=124)
fitSepModel <- fit.StVariogram(empVgm, separableModel, fit.method = 7,
stAni = linStAni, method = "L-BFGS-B",
control = list(parscale=c(100,1,10,1,100)),
lower = c(10,0,.1,0,0.1),
upper = c(2000,1,12,1,200))
fitSepModel## space component:
## model psill range
## 1 Nug 0.1417641 0.0000
## 2 Exp 0.8582359 558.1318
## time component:
## model psill range
## 1 Nug 0 0.000000
## 2 Sph 1 5.581331
## sill: 123.940491343563Modelo product-sum
\[ C_{ps}\left(h,\mu\right) = kC_{s}\left(h\right)C_{t}\left(\mu\right) + C_{s}\left(h\right) + C_{t}\left(\mu\right) \] com \(k > 0\). O variograma correspondente pode ser escrito como: \[ \gamma_{ps}\left(h,\mu\right) = \left(k\cdot sill_{t} + 1\right)\gamma_{s}(h) + \left(k\cdot sill_{s} + 1\right)\gamma_{t}(\mu) - k\gamma_{s}(h)\gamma_{t}(\mu) \]
prodSumModel <- vgmST("productSum",
space=vgm(10, "Exp", 200, 1),
time= vgm(10, "Sph", 2, 1),
k=2)
fitProdSumModel <- fit.StVariogram(empVgm, prodSumModel, fit.method = 7,
stAni = linStAni, method = "L-BFGS-B",
control = list(parscale =
c(1,10,1,1,0.1,1,10)),
lower = rep(0.0001,7))
fitProdSumModel## space component:
## model psill range
## 1 Nug 1.199504 0.0000
## 2 Exp 6.851819 543.4864
## time component:
## model psill range
## 1 Nug 0.000100 0.000000
## 2 Sph 8.721038 5.503135
## k: 1.59118539176836Modelo sum-metric
\[ C_{sm}(h,\mu) = C_{s}(h) + C_{t}(\mu) + C_{joint}\left(\sqrt{h^{2} + \left(\mathbf{k}\cdot\mu\right)^2}\right) \] Este modelo permite efeitos pepita espaciais, temporais e conjuntos. Assim, o variograma é dado por \[ \gamma_{sm}(h,\mu) = \gamma_{s}(h) + \gamma_{t}(\mu) + \gamma_{joint}\left(\sqrt{h^{2} + \left(\mathbf{k}\cdot\mu\right)^2}\right) \]
#sumMetricModel <- vgmST("sumMetric",
# space = vgm(20, "Sph", 150, 1),
# time = vgm(10, "Exp", 2, 0.5),
# joint = vgm(80, "Sph", 1500, 2.5),
# stAni = 120)
#fitSumMetricModel <- fit.StVariogram(empVgm, sumMetricModel,
# fit.method = 7, stAni=linStAni,
# method = "L-BFGS-B",
# lower = c(sill.s = 0, range.s = 10, nugget.s = 0,
# sill.t = 0, range.t = 0.1, nugget.t = 0,
# sill.st= 0, range.st = 10, nugget.st = 0,
# anis = 40),
# upper = c(sill.s = 200, range.s = 1E3, nugget.s = 20,
# sill.t = 200, range.t = 75, nugget.t = 20,
# sill.st= 200, range.st = 5E3, nugget.st = 20,
# anis = 500),
# control = list(parscale = c(1,100,1,1,0.5,1,1,100,1,100),
# maxit=1e4))fitSumMetricModel## space component:
## model psill range
## 1 Nug 0.00000 0.0000
## 2 Sph 16.42764 67.3639
## time component:
## model psill range
## 1 Nug 0.00000 0.0000000
## 2 Exp 9.31931 0.9122949
## joint component:
## model psill range
## 1 Nug 7.261825 0.0000
## 2 Sph 91.515174 998.8057
## stAni: 184.58922445157Selecionando o melhor modelo
O melhor modelo que se ajusta pode ser feito baseando em um critério de validação cruzada na krigagem. Aqui iremos selecionar, por conveniência, por meio de uma análise visual:
plot(empVgm, list(fitSepModel, fitProdSumModel, fitSumMetricModel),
wireframe=T, all=T, zlim=c(0,140), ylim=c(0,6.1), xlim=c(0,300),
scales=list(arrows = F, cex=.8,
x=list(at=0:3*100),
y=list(at=0:6, labels=c("0 ","","2 ","","4 ","","6 ")),
z=list(at=0:5*25,
labels=c("0 ","","50 ","","100 ",""))),
at=0:100*1.4,xlab=list("space [km]", rot=27, cex=0.8),
ylab=list("time [days]", rot=-40, cex=0.8),
zlab=list(NULL, rot=94, cex=0.8))
plot(empVgm, list(fitSepModel, fitProdSumModel, fitSumMetricModel),
wireframe=T, all=T, zlim=c(-10,25), ylim=c(0,6.1), xlim=c(0,300),
diff=TRUE,
scales=list(arrows = F, cex=.8,
x=list(at=0:3*100),
y=list(at=0:6, labels=c("0 ","","2 ","","4 ","","6 ")),
z=list(at=-2:5*5,
labels=c("-10 ","","0 ","","10 ","","20 ",""))),
xlab=list("space [km]", rot=27, cex=0.8),
ylab=list("time [days]", rot=-40, cex=0.8),
zlab=list(NULL, rot=94, cex=0.8))
Krigagem
data(DE_RB_2005)
gride <- SpatialGrid(GridTopology(DE_RB_2005@sp@bbox[,1]%/%10000*10000,
c(10000,10000),
cells.dim=ceiling(apply(
DE_RB_2005@sp@bbox,1,diff)/10000)))
proj4string(gride) <- CRS("+init=epsg:32632")
fullgrid(gride) <- F
data(air)
DE_NUTS1 <- spTransform(DE_NUTS1, CRS("+init=epsg:32632"))
ind <- over(gride, as(DE_NUTS1,"SpatialPolygons"))
gridDE <- gride[!is.na(ind)]plot(gridDE)
Predição Espaço-Temporal
#sumPred <- krigeST(PM10~1, data=DE_RB_2005[,tIDS],
# newdata=DE_pred, fitSumMetricModel, nmax=50,
# stAni=fitSumMetricModel$stAni/24/3600)stplot(sumPred, col.regions=bpy.colors(120)[-(1:20)],
scales=list(draw=F),
main="spatio-temporal sum-metric model", at=0:70,
sp.layout = list(list("sp.polygons", DE_NUTS1, first=FALSE,
col=gray(0.5)),
list("sp.points", DE_RB_2005@sp, col=gray(0.25), pch=3, cex=.5)))