Work 7
Faythe
May 15, 2016
Inverse distance weighting
Below I attempt to cross validate the IDW power values, but was unsure how to compare my predicted values to my actual values similar to what I did for the krige validation. I applied four different power values: 0.25,0.75, 1.25, and 3.33. They all look very different to me.
setwd("E:/spring16")
library(sp)
library(raster)
load("wUSPrec_Bnd.Rdata")
spplot(wUSPrec,"prec",colorkey=F,main="Precipitation")
e <- extent(wUSPrec)
wUSPrec.grid <- raster(e)
proj4string(wUSPrec.grid) <- proj4string(wUSPrec)
res(wUSPrec.grid) <- c(1e4,1e4)
wUSPrec.grid <- as(wUSPrec.grid,"SpatialPointsDataFrame")
plot(wUSPrec.grid,col="grey")
points(wUSPrec,pch=20,col="darkmagenta")
summary(wUSPrec)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## long -302251.8 381330.4
## lat -422252.6 556846.7
## Is projected: TRUE
## proj4string :
## [+proj=aea +lat_1=41 +lat_2=47 +lat_0=44 +lon_0=-120 +x_0=0 +y_0=0
## +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0]
## Number of points: 113
## Data attributes:
## prec
## Min. : 23.69
## 1st Qu.: 40.88
## Median : 60.92
## Mean : 75.00
## 3rd Qu.: 97.62
## Max. :272.76##I subsampled here to try to cross validate but was unsure how to compare predicted and actural values...some more tinkering is needed.
library(gstat)
sam75=sample(1:113,75)
p.model=wUSPrec[sam75,]
p.valid=wUSPrec[-sam75,]
e1 <- extent(p.model)
wUSPrec.grid1 <- raster(e1)
proj4string(wUSPrec.grid1) <- proj4string(p.model)
res(wUSPrec.grid1) <- c(1e4,1e4)
wUSPrec.grid1 <- as(wUSPrec.grid1,"SpatialPointsDataFrame")
prec.idw1 <- idw(formula=prec~1,locations=p.model,newdata=wUSPrec.grid1,idp=0.5)## [inverse distance weighted interpolation]spplot(prec.idw1,"var1.pred",main="p=0.5")
prec.idw1## class : SpatialPointsDataFrame
## features : 6468
## extent : -297251.8, 352748.2, -418992.8, 551007.2 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=41 +lat_2=47 +lat_0=44 +lon_0=-120 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0
## variables : 2
## names : var1.pred, var1.var
## min values : 60.4535980525424, Inf
## max values : 88.6445647566772, -Infresid.idw=p.valid$prec-prec.idw1$var1.pred## Warning in p.valid$prec - prec.idw1$var1.pred: longer object length is not
## a multiple of shorter object lengthsummary(resid.idw)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -62.570 -27.360 -6.273 8.257 39.830 208.700##Here are my different power plots.
prec.idw2 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=0.25)## [inverse distance weighted interpolation]spplot(prec.idw2,"var1.pred",main="p=0.25")
prec.idw3 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=0.75)## [inverse distance weighted interpolation]spplot(prec.idw3,"var1.pred",main="p=0.75")
prec.idw4 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=1.25)## [inverse distance weighted interpolation]spplot(prec.idw2,"var1.pred",main="p=1.25")
prec.idw5 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=3.33)## [inverse distance weighted interpolation]spplot(prec.idw5,"var1.pred",main="p=3.33")
Kriging
First I fit a variogram the good’ol fashion way, then I validate a model. When I used the auto version, my estimates were not too far off those r chose.
precVar <- variogram(prec~1, wUSPrec)
summary(precVar)## np dist gamma dir.hor dir.ver
## Min. : 14.0 Min. : 19101 Min. : 690.9 Min. :0 Min. :0
## 1st Qu.:157.0 1st Qu.:106259 1st Qu.:1380.6 1st Qu.:0 1st Qu.:0
## Median :239.0 Median :199112 Median :1901.7 Median :0 Median :0
## Mean :199.6 Mean :199504 Mean :1691.5 Mean :0 Mean :0
## 3rd Qu.:259.5 3rd Qu.:291831 3rd Qu.:1981.7 3rd Qu.:0 3rd Qu.:0
## Max. :291.0 Max. :384270 Max. :2327.3 Max. :0 Max. :0
## id
## var1:15
##
##
##
##
## plot(precVar,pch=15,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Precipitation")
sph.model <- vgm(psill=2000, model="Sph", range=250000, nugget=500)
sph.fit <- fit.variogram(object = precVar, model = sph.model)
plot(precVar,pch=15,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Precipitation",model=sph.fit)
precSph <- krige(prec~1, wUSPrec, wUSPrec.grid, model = sph.fit)## [using ordinary kriging]spplot(precSph,"var1.pred")
spplot(precSph,"var1.var")
cross validation
I took a subsample of 75 to create a model and then tested it with the remaining 38. An R2 was then calculated.
sam75=sample(1:113,75)
p.model=wUSPrec[sam75,]
p.valid=wUSPrec[-sam75,]
v75.fit=fit.variogram(variogram(prec~1,p.model), vgm(psill=2000, model="Sph", range=250000, nugget=500))
m.valid.pr=krige(prec~1,p.model,p.valid, model=v75.fit)## [using ordinary kriging]resid.kr=p.valid$prec-m.valid.pr$var1.pred
summary(resid.kr)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -63.820 -14.480 -4.942 -3.120 3.782 120.400resid.mean=p.valid$prec-mean(p.valid$prec)
R2=1-sum(resid.kr^2)/sum(resid.mean^2)
R2## [1] 0.3572943bubble(m.valid.pr,"var1.var")
When I ran this the first time my R2 was 0.44.
Auto variogram and Kriging
library(automap)
pprecVar <- autofitVariogram(prec~1,wUSPrec)
summary(pprecVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 13 18557.98 633.7695 0 0 var1
## 2 21 33256.55 667.3642 0 0 var1
## 3 30 44106.34 903.7051 0 0 var1
## 4 59 56966.96 951.1490 0 0 var1
## 5 184 85646.40 1417.2941 0 0 var1
## 6 276 125642.37 1569.0008 0 0 var1
## 7 530 178543.95 1832.8309 0 0 var1
## 8 586 240222.51 2141.5306 0 0 var1
## 9 622 303366.38 1971.5391 0 0 var1
## 10 899 375807.17 1949.7610 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 538.1447 0.0 0
## 2 Ste 1486.5017 113262.1 2
## Sums of squares betw. var. model and sample var.[1] 0.0008848358plot(pprecVar,pch=12)
##My estimates were not too far off.
##Autokriging
precKrig <- autoKrige(prec~1, wUSPrec, wUSPrec.grid)## [using ordinary kriging]plot(precKrig)