## ---- include=FALSE----------------------------------------------------------- library(knitr) opts_chunk$set( fig.width = 5, fig.height = 5, fig.cap = '', collapse = TRUE ) library(rspatial) library(rasterVis) library(rgdal) library(deldir) ## ----------------------------------------------------------------------------- z <- function(x, y) { -12 * x^3 + 10 * x^2 * y - 14 * x * y^2 + 25 * y^3 + 50 } z(.5, 0.8) ## ----------------------------------------------------------------------------- zf <- function(model, xy) { x <- xy[,1] y <- xy[,2] z(x, y) } ## ----fields5------------------------------------------------------------------ library(raster) 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) ## ----fields6------------------------------------------------------------------ pts <- rasterToPoints(z) pt <- trans3d(pts[,1], pts[,2], pts[,3], vt) plot(pt, col=rainbow(9, .75, start=.2)[round(pts[,3]/10)-2], pch=20, cex=2) ## ----fields10----------------------------------------------------------------- d <- sp_data("precipitation.csv") head(d) d$prec <- rowSums(d[, c(6:17)]) plot(sort(d$prec), ylab='Annual precipitation (mm)', las=1, xlab='Stations') ## ----fields15----------------------------------------------------------------- if (!require("rspatial")) remotes::install_github('rspatial/rspatial') library(rspatial) dsp <- SpatialPoints(d[,4:3], proj4string=CRS("+proj=longlat +datum=NAD83")) dsp <- SpatialPointsDataFrame(dsp, d) CA <- sp_data("counties.rds") cuts <- c(0,200,300,500,1000,3000) pols <- list("sp.polygons", CA, fill = "lightgray") # 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) ## ----------------------------------------------------------------------------- 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") library(rgdal) dta <- spTransform(dsp, TA) cata <- spTransform(CA, TA)