Spatial04
Sara Cendejas-Zarelli
May 10, 2016
library(sp)
library(gstat)
library(ggplot2)
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(raster)Inverse distance weighting
load("wUSPrec_Bnd.Rdata")
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.76sp_map <- fortify(wUSBoundary)## Regions defined for each PolygonswUSPrec.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="orange")
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 weights
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="red")
Here is a grid map of the points overlaying the entire area. Now I can interpolate points over the rest of the grid.
Now I’ll create a spatial points data frame of the IDW with a power of 2.5.
precIDW1 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=2.5)## [inverse distance weighted interpolation]head(precIDW1@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 NAspplot(precIDW1,"var1.pred",main="p=2.5")
precIDW3.5 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=3.5)## [inverse distance weighted interpolation]spplot(precIDW3.5,"var1.pred",main="p=3.5")
precIDW2 <- idw(formula=prec~1,locations=wUSPrec,newdata=wUSPrec.grid,idp=1)## [inverse distance weighted interpolation]spplot(precIDW2,"var1.pred", main="p=1")
As the value of p shrinks, more distant points have more of an effect on the interpolated points. A highter p means that closer points have a greater weight. Thus IDW with a low p tend to produce more accurate plots with less error. I guess the “right” value may be somewhat objective…but should be somewhere around the default that R uses of 2.
Kriging
I couldn’t figure out the shiny app in time. But, I did create a plot and highlighted some outliers. I seems to me that outliers have the potential to effect the accuracy of the interpolation. One could choose to take them out if it were deemed necessary.
show.vgms(max=1500, models = c("Exp", "Mat", "Gau", "Sph"), sill=15, nugget = 5, range=1000)
This looks like the spherical model is best. I’ll stick with that.
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=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Summer precipitation (mm)")
The variogram shows spatial autocorrelation out to about 2.5 e05m (250km). Now, let’s fit the variogram model. I’m using a sill of 2500, a range of 2.5e+05, and a nugget of 600 to start things off.
sph.model <- vgm(psill=2500, model="Sph", range=2.5e+05, nugget=600)
sph.fit <- fit.variogram(object = precvar, model = sph.model)
plot(precvar,pch=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Summer precipitation (mm)",model=sph.fit)
sph.fit## model psill range
## 1 Nug 564.7154 0.0
## 2 Sph 1479.0279 251499.7This shows the nugget at 564, the sill at 1,479, and the range at 251,500
Here’s an exponential model just for kicks…this model is more of a gradual curve and doesn’t really highlight the sill as well as the spherical model.
sph.model <- vgm(psill=2500, model="Exp", range=2.5e+05, nugget=600)
sph.fit <- fit.variogram(object = precvar, model = sph.model)
plot(precvar,pch=20,cex=1.5,col="black",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Summer precipitation (mm)",model=sph.fit)
Let’s get the weights now from the spherial model…
precsph <- krige(prec~1, wUSPrec, wUSPrec.grid, model = sph.fit)## [using ordinary kriging]spplot(precsph,"var1.pred")
I’m going to try the automatic kriging for kicks…
library(automap)
precvar <- autofitVariogram(prec~1,wUSPrec)plot(precvar)
The automatic kriging shows the sill at 2025, the range at 1.13e+05, and the nugget at 538. These numbers are definitely different than the by hand version of kriging.
Now I’ll krige a surface automatically…
precKrig <- autoKrige(prec~1, wUSPrec, wUSPrec.grid)## [using ordinary kriging]plot(precKrig)
I see high error in those areas with few data points. But, I’m not sure how to go about solving this issue. Maybe adjust the neighborhood of points?