Inverse Distance Weighting and Finding P

library(sp)
library(gstat)
library(raster)
library(ggmap)
## Loading required package: ggplot2
library(ggalt)
library(colorspace)
## 
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
## 
##     RGB
library(gridExtra)
birds <- readRDS("birdRichnessMexico.rds")
gridMex <- readRDS("gridMex.rds")

spplot(birds,"nSpecies",colorkey=TRUE)

Inverse distance weights with p=2.5

mex.idw <- gstat::idw(nSpecies ~ 1, birds, gridMex, idp = 3)
## [inverse distance weighted interpolation]
spplot(mex.idw,"var1.pred")

#Try different value of p and compare the R squared for your models to find the best fit.

mex.idw2<- gstat::idw(nSpecies ~ 1, birds, gridMex, idp = 2.5)
## [inverse distance weighted interpolation]
spplot(mex.idw2,"var1.pred")

Making Variograms

birdsVar <- variogram(birds$nSpecies~1, birds)
summary(birdsVar)
##        np              dist             gamma           dir.hor 
##  Min.   : 155.0   Min.   :  50455   Min.   : 237.5   Min.   :0  
##  1st Qu.: 757.5   1st Qu.: 298845   1st Qu.:1908.6   1st Qu.:0  
##  Median : 893.0   Median : 558141   Median :2440.4   Median :0  
##  Mean   : 810.0   Mean   : 559756   Mean   :2433.7   Mean   :0  
##  3rd Qu.: 956.0   3rd Qu.: 818467   3rd Qu.:3248.7   3rd Qu.:0  
##  Max.   :1083.0   Max.   :1078585   Max.   :4111.7   Max.   :0  
##     dir.ver     id    
##  Min.   :0   var1:15  
##  1st Qu.:0            
##  Median :0            
##  Mean   :0            
##  3rd Qu.:0            
##  Max.   :0
plot(birdsVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Number of Bird Species Mexico")

The sill shows up at about 800 kilometers with a nugget of 100 kilometers

sph.model <- vgm(psill=800000, model="Sph", range=700000, nugget=100000)
sph.fit <- fit.variogram(object = birdsVar, model = sph.model)
plot(birdsVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Birds Species richness Mexico",model=sph.fit)

exp.model <- vgm(psill=800000, model="Exp", range=700000, nugget=100000)
exp.fit <- fit.variogram(object = birdsVar, model = exp.model)
plot(birdsVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Birds Species richness Mexico",model=exp.fit)

gau.model <- vgm(psill=800000, model="Gau", range=700000, nugget=100000)
gau.fit <- fit.variogram(object = birdsVar, model = gau.model)
plot(birdsVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Birds Species richness Mexico",model=gau.fit)

The spherical model did not fit the semivariance as well as an exponential model.

Kriging

 birdsSph <- krige(birds$nSpecies~1, birds, gridMex, model = sph.fit)
## [using ordinary kriging]
 spplot(birdsSph,"var1.pred")

 spplot(birdsSph,"var1.var")

birdsexp <- krige(birds$nSpecies~1, birds, gridMex, model = exp.fit)
## [using ordinary kriging]
 spplot(birdsexp,"var1.pred")

 spplot(birdsexp,"var1.var")

birdsgau <- krige(birds$nSpecies~1, birds, gridMex, model = gau.fit)
## [using ordinary kriging]
 spplot(birdsgau,"var1.pred")

 spplot(birdsgau,"var1.var")

The standard error is smaller in areas where there are actual data points. The Gausian model has the lowest standard error eventhough the variogram model fit looked worse. The species richness is predicted to be higher in the south which is similar to what the idw analysis predicted.

library(automap)
birdsKrig <- autoKrige(nSpecies~1, birds, gridMex)
## [using ordinary kriging]
plot(birdsKrig)