9. Fields

9.1 Introduction

This handout accompanies Chapter 9 in O’Sullivan and Unwin (2010).

Download the data for this lesson here: download here: counties, precipitation. Or with the script below.

dir.create('data', showWarnings = FALSE)
files <- c('precipitation.csv', 'counties.rds')
for (filename in files) {
    localfile <- file.path('data', filename)
    if (!file.exists(localfile)) {
        download.file(paste0('https://biogeo.ucdavis.edu/data/rspatial/', filename), dest=localfile)

Here is how you can set up and use the continuous function on page 246.

z <- function(x, y) { -12 * x^3 + 10 * x^2 * y - 14 * x * y^2 + 25 * y^3 + 50 }
z(.5, 0.8)
## [1] 58.82

Function zf adds some complexity to make it usable in the ‘interpolate’ function below.

zf <- function(model, xy) {
   x <- xy[,1]
   y <- xy[,2]
   z(x, y)

Now use it

r <- raster(xmn=0.5, xmx=1.4, ymn=0.6, ymx=1.5, ncol=9, nrow=9, crs=NA)
z <- interpolate(r, model=NULL, fun=zf)
names(z) <- 'z'

vt <- persp(z, theta=30, phi=30, ticktype='detailed', expand=.8)


Note that persp returned something invisibly (it won’t be printed when not captured as a variable, vt, in this case), the 3D transformation matrix that we use later. This is not uncommon in R. For example hist and barplot have similar behaviour.

pts <- rasterToPoints(z)
pt <- trans3d(pts[,1], pts[,2], pts[,3], vt)
points(pt, col=rainbow(9, .75, start=.2)[round(pts[,3]/10)-2], pch=20, cex=2)
## Error in plot.xy(xy.coords(x, y), type = type, ...): plot.new has not been called yet

For a more interactive experience, try:

# this opens a new window
plot3D(z, zfac=5)

We will be working with temperature data for California. You can download the climate data used in the examples.

d <- read.csv(file.path(datapath,"precipitation.csv"))
## Error in file.path(datapath, "precipitation.csv"): object 'datapath' not found
## Error in head(d): object 'd' not found
d$prec <- rowSums(d[, c(6:17)])
## Error in rowSums(d[, c(6:17)]): object 'd' not found
plot(sort(d$prec), ylab='Annual precipitation (mm)', las=1, xlab='Stations')
## Error in sort(d$prec): object 'd' not found

A quick map (get the counties here).

dsp <- SpatialPoints(d[,4:3], proj4string=CRS("+proj=longlat +datum=NAD83"))
## Error in coordinates(coords): object 'd' not found
dsp <- SpatialPointsDataFrame(dsp, d)
## Error in is(coords, "SpatialPoints"): object 'dsp' not found
CA <- readRDS("data/counties.rds")
## Warning in readRDS("data/counties.rds"): invalid or incomplete compressed
## data
## Error in readRDS("data/counties.rds"): error reading from connection

cuts <- c(0,200,300,500,1000,3000)
pols <- list("sp.polygons", CA, fill = "lightgray")
## Error in eval(expr, envir, enclos): object 'CA' not found
# set up a palette of interpolated colors
blues <- colorRampPalette(c('yellow', 'orange', 'blue', 'dark blue'))
spplot(dsp, 'prec', cuts=cuts, col.regions=blues(5), sp.layout=pols, pch=20, cex=2)
## Error in spplot(dsp, "prec", cuts = cuts, col.regions = blues(5), sp.layout = pols, : object 'dsp' not found

Transform longitude/latitude to planar coordinates

TA <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +ellps=GRS80 +towgs84=0,0,0")
dta <- spTransform(dsp, TA)
## Error in spTransform(dsp, TA): object 'dsp' not found
cata <- spTransform(CA, TA)
## Error in spTransform(CA, TA): object 'CA' not found

Continue here