Interpolation_answers

Set up

Load the following packages:

library(gstat)
library(tidyverse)
library(sp)
library(raster)

…and the precip.txt file and the counties.RDS file that stores the county shapefile for California:

precip <- read.table("C:/Users/Maggie/Downloads/precip.txt", sep = ",", header=T)
ca <- readRDS("C:/Users/Maggie/Downloads/counties.RDS")

Q1

• Use spplot() to visualize precipitation in April. Note that to do this, you will have to coerce the .txt file to a SpatialPointsDataFrame(). We’ll use this SPDF for the rest of the tutorial.

# create spdf
p1 <- CRS("+proj=longlat +datum=NAD83") # original projection
prec <- SpatialPointsDataFrame(coords=cbind(precip$LONG,precip$LAT), proj4string=p1, data=precip)
# pull out april
apr <- prec[,c(1:5,9)]
# plot april
spplot(apr,"APR")

Q2

Perform Inverse Distance Weighting interpolation on your April precipitation data using a grid with n = 5000 cells. As you build your grid from the counties.RDS file, be sure to reproject the data to the projection used in the tutorial:

• CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

# reproject
proj <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
apr <- spTransform(apr,proj)
ca <- spTransform(ca,apr@proj4string)
plot(ca)

ca@proj4string; apr@proj4string

## CRS arguments:
##  +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0
## +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
## +towgs84=0,0,0

## CRS arguments:
##  +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0
## +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
## +towgs84=0,0,0

# make an empty grid
grd <- SpatialPixels(SpatialPoints(makegrid(ca, n=5000)), proj4string=proj4string(ca))
grd <- grd[ca]
plot(grd)

# perform IDW with standard inputs
p.idw <- idw(formula=APR~1, apr, grd, idp=2, nmax=12)

## [inverse distance weighted interpolation]

plot(p.idw) # use raster extract at known pts to calc RMSE

Q3

Create a sample variogram for the APR attribute of the SPDF you created. Assume you are not worried about outlier observations.

p.vgm <- variogram(APR~1, apr)
plot(p.vgm)

Q4

Fit a model variogram to your data. Justify your choice. * psill: point at which p.vgm levels off on y-axis * range: point at which p.vgm levels off on x-axis * nugget: point at which p.vgm crosses y-axis (x=0) * eyeballed and tried multiple models, exp fit the best

show.vgms()

p.fit <- fit.variogram(p.vgm, model=vgm(psill=850, "Exp", range=400000, nugget = 50))
plot(p.vgm, pch=20, cex=1.5, col="black", ylab="Semivariance", xlab="Distance (m)", model=p.fit)

Q5

Krige your data using your model variogram.

p.krig <- krige(APR~1, apr, grd, model=p.fit)

## [using ordinary kriging]

plot(p.krig)