Week 7 Work
setwd("/Users/shannoncall/Desktop/ESCI597A/Week07")
load("wUSPrec_Bnd.Rdata")
library(sp)
library(gstat)
library(raster)
library(maptools)
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()
library(ggplot2)
library(automap)
Precipitation Data Map
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
sp_map <-fortify(wUSBoundary)
## Regions defined for each Polygons
wUSPrec.df <- data.frame(wUSPrec)
gg <- ggplot()
gg <- gg + geom_map(data=sp_map,map=sp_map,
aes(x=long, y=lat, map_id=id),
fill= "white",color="black")
gg <- gg + geom_point(data=wUSPrec.df, aes(x=long, y=lat, size=prec),
color="dodgerblue3")
gg <- gg + scale_size(breaks=seq(0,300,by=50))
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + theme(panel.grid=element_blank(), panel.border=element_blank())
gg <- gg + theme(axis.ticks=element_blank(), axis.text=element_blank())
gg <- gg + theme(legend.position="top")
gg
Inverse Distance Weighted Interpolation
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="dodgerblue3")
Spatial points data frame with the inverse distance weight of 0.5
precIDW.5 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=0.5)
## [inverse distance weighted interpolation]
head(precIDW.5@data)
## var1.pred var1.var
## 1 81.88799 NA
## 2 82.10247 NA
## 3 82.30794 NA
## 4 82.46192 NA
## 5 82.55687 NA
## 6 82.63981 NA
spplot(precIDW.5,"var1.pred",main="p=0.5")
Spatial points data frame with the inverse distance weight of 1
precIDW1 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=1)
## [inverse distance weighted interpolation]
head(precIDW1@data)
## var1.pred var1.var
## 1 88.76503 NA
## 2 89.15089 NA
## 3 89.47240 NA
## 4 89.61433 NA
## 5 89.63052 NA
## 6 89.76907 NA
spplot(precIDW1,"var1.pred",main="p=1")
Spatial points data frame with the inverse distance weight of 1.5
precIDW1.5 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=1.5)
## [inverse distance weighted interpolation]
head(precIDW1.5@data)
## var1.pred var1.var
## 1 93.39316 NA
## 2 93.54956 NA
## 3 93.54806 NA
## 4 93.25921 NA
## 5 92.88160 NA
## 6 93.01948 NA
spplot(precIDW1.5,"var1.pred",main="p=1.5")
Spatial points data frame with the inverse distance weight of 2
precIDW2 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=2)
## [inverse distance weighted interpolation]
head(precIDW2@data)
## var1.pred var1.var
## 1 94.95555 NA
## 2 94.77933 NA
## 3 94.42209 NA
## 4 93.65805 NA
## 5 92.70404 NA
## 6 92.56109 NA
spplot(precIDW2,"var1.pred",main="p=2")
Spatial points data frame with the inverse distance weight of 2.5
precIDW2.5 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=2.5)
## [inverse distance weighted interpolation]
head(precIDW2.5@data)
## var1.pred var1.var
## 1 94.97661 NA
## 2 94.77424 NA
## 3 94.34023 NA
## 4 93.26673 NA
## 5 91.69603 NA
## 6 90.97853 NA
spplot(precIDW2.5,"var1.pred",main="p=2.5")
Spatial points data frame with the inverse distance weight of 16
precIDW16 <- idw(formula=prec~1,locations=wUSPrec, newdata=wUSPrec.grid,idp=16)
## [inverse distance weighted interpolation]
head(precIDW16@data)
## var1.pred var1.var
## 1 98.19172 NA
## 2 98.19732 NA
## 3 97.84969 NA
## 4 95.92313 NA
## 5 88.80136 NA
## 6 82.01754 NA
spplot(precIDW16,"var1.pred",main="p=16")
The “right” value seems to be an objective, with a higher p value, the points closer seem to have more weight.
Kriging
precVar <- autofitVariogram(prec~1,wUSPrec)
summary(precVar)
## 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.0008848358
plot(precVar)
Because I’m a monster…
precKrig <- autoKrige(prec~1, wUSPrec, wUSPrec.grid)
## [using ordinary kriging]
plot(precKrig)
Max at 1500
show.vgms(max=1500, models = c("Exp", "Mat", "Gau", "Sph"), sill=15, nugget = 5, range=1000)
precVar<- variogram(prec~1, wUSPrec)
summary(precVar)
## np dist gamma dir.hor dir.ver
## Min. : 14.0 Min. : 19100 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=21,cex=1.5,col="dodgerblue4",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Summer precipitation (mm)")
sph.model <- vgm(psill=2000, model="Sph", range=3e05, nugget=5)
sph.fit <- fit.variogram(object = precVar, model = sph.model)
plot(precVar,pch=21,cex=1.5,col="dodgerblue4",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Summer precipitation (mm)",model=sph.fit)
sph.fit
## model psill range
## 1 Nug 564.7484 0.0
## 2 Sph 1479.0539 251526.2
Cross Validation to come… I am a little iffy about the code.
