## ----setup, include=FALSE, cache=FALSE------------------------------ # set global chunk options ## this code is only required to set up knitr for reproducible documentation ## it has nothing to do with the exercise itself, please ignore! opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), out.format='latex', use.highlight=TRUE) opts_chunk$set(fig.path='kgraph/exRKGLS', fig.align='center', fig.show='hold', prompt=FALSE, fig.width=6, fig.height=6, out.width='.7\\linewidth', size='scriptsize') # when fig.cap needs to use objects in the current chunk opts_knit$set(eval.after = 'fig.cap') ## comment=NA, ## this does not comment out the output ## options to be read by formatR options(replace.assign=TRUE,width=70) ## supress warnings about rgdal being archived. options("rgdal_show_exportToProj4_warnings"="none") ## ----echo=FALSE, results='hide'------------------------------------- rm(list=ls()) ## ----child-data, child='MappingRegionalClimate_explore.Rnw'--------- ## ----explore-set-parent, echo=FALSE, cache=FALSE-------------------- set_parent('MappingRegionalClimate.Rnw') ## ----load-packages-explore------------------------------------------ require(sf) ## ----load-saved-dataset, eval=FALSE--------------------------------- ## load(file="./StationsDEM_covariates.RData", verbose=TRUE) ## ----load-saved-dataset-yes, eval=TRUE, echo=FALSE------------------ rm(list=ls()) # make sure this is all the data we need to go forward load(file="./ds/NEweather/StationsDEM_covariates.RData", verbose=TRUE) ## ----summarize-ds--------------------------------------------------- summary(ne.df[,c("E", "N", "ANN_GDD50","ELEVATION_")]) ## ----show-stations, out.width='\\linewidth', w=8, h=8, fig.cap="Location of weather stations"---- unique(st_geometry_type(state.ne.m)) state.ne.m.boundary <- st_geometry(st_cast(state.ne.m, "MULTILINESTRING")) unique(st_geometry_type(state.ne.m.boundary)) plot(st_coordinates(ne.m), asp=1, xlab="E", ylab="N", pch=20) text(st_coordinates(ne.m), labels=substr(ne.m$STATION_NA, 1, 4), cex=0.5, pos=2, col=as.numeric(ne.m$STATE)) plot(state.ne.m.boundary, add=TRUE, col="darkgray", lwd=2, sf_max.plot = 1) grid() ## ----load-ggplot2--------------------------------------------------- require(ggplot2) ## ----display-annual, out.width='\\textwidth', w=11, h=10, fig.cap="Annual Growing Degree days, 50F"---- ggplot(data=ne.df) + aes(x=E, y=N) + geom_point(aes(size=ANN_GDD50, colour=ELEVATION_), shape=20) + xlab("E") + ylab("N") + coord_fixed() ## ----display-annual-sf, out.width='\\textwidth', w=11, h=10, fig.cap="Annual Growing Degree days, base 50F"---- ggplot() + geom_sf(data = ne.m, aes(size=ANN_GDD50, colour=ELEVATION_), shape = 20) + labs(x = "Longitude", y = "Latitude", size = "GDD base 50F", colour = "Elevation, feet a.s.l.") + geom_sf(data = state.ne.m.boundary, col = "darkgray", size = 2) ## ----display-annual-sf-2, out.width='\\textwidth', w=11, h=10, fig.cap="Annual Growing Degree days, base 50F"---- (which.lowest.gdd <- which.min(ne.m$ANN_GDD50)) ggplot() + geom_sf(data = ne.m[-which.lowest.gdd, ], aes(size=ANN_GDD50, colour = ELEVATION_), shape = 10) + geom_sf(data = ne.m[which.lowest.gdd, ], colour = "red") + labs(x = "Longitude", y = "Latitude", size = "GDD base 50F", colour = "Elevation, feet a.s.l.") + geom_sf(data = state.ne.m.boundary, col = "darkgray", size = 2) ## ----export-to-kml-------------------------------------------------- ne.wgs84 <- st_transform(ne.m, st_crs(4326)) # this is the EPSG code for global WGS84 long/lat require(plotKML) shape = "http://maps.google.com/mapfiles/kml/pal2/icon18.png" station.names = substr(ne.wgs84$STATION_NA,1,12) kml_open(file.name='ne.stations.kml') kml_layer(ne.wgs84, colour=ANN_GDD50, colour_scale=SAGA_pal[[1]], shape=shape, points_names=station.names) kml_close('ne.stations.kml') ## ----child-data, child='MappingRegionalClimate_OLS.Rnw'------------- ## ----OLS-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('MappingRegionalClimate.Rnw') ## ----scatter-covars------------------------------------------------- p1 <- ggplot() + geom_point( aes(x=ELEVATION_, y=ANN_GDD50, colour=STATE), data=ne.df ) p2 <- ggplot() + geom_point( aes(x=E, y=ANN_GDD50, colour=STATE), data=ne.df ) + xlab("Easting") p3 <- ggplot() + geom_point( aes(x=N, y=ANN_GDD50, colour=STATE), data=ne.df ) + xlab("Northing") ## ----print-scatter-covars, out.width='0.31\\linewidth'-------------- print(p3) print(p2) print(p1) ## ----scatter-xform, out.width='0.48\\linewidth', fig.cap="GDD50 vs. square root of elevation"---- ggplot() + geom_point( aes(x=sqrt(ELEVATION_), y=ANN_GDD50, colour=STATE), data=ne.df ) ## ----cor-preds------------------------------------------------------ with(ne.df, cor(ANN_GDD50, N)) with(ne.df, cor(ANN_GDD50, E)) with(ne.df, cor(ANN_GDD50, sqrt(ELEVATION_))) ## ----build-naive-elev----------------------------------------------- m.ols.elev <- lm(ANN_GDD50 ~ sqrt(ELEVATION_), data=ne.df) summary(m.ols.elev) ## ----build-naive-two------------------------------------------------ m.ols <- lm(ANN_GDD50 ~ sqrt(ELEVATION_) + N, data=ne.df) summary(m.ols) ## ----plot-lm-fits, fig.cap = "Actual vs. predicted, OLS model"------ plot(ne.m$ANN_GDD50 ~ fitted(m.ols), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by OLS", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----dx-naive, out.width='\\textwidth', w=12, h=4------------------- par(mfrow=c(1,3)) plot(m.ols, which=c(1,2,5)) par(mfrow=c(1,1)) ## ----find-min-max-ols-resid----------------------------------------- (ix <- c(which.max(residuals(m.ols)), which.min(residuals(m.ols)))) ne.df[ix,] ne.df[ix,c("LATITUDE_D", "LONGITUDE_")] ## ----correct-it----------------------------------------------------- ne.m[ix[1],"ELEVATION_"] <- ne.m[ix[1],"ELEV_FT"] <- 1192 ne.df[ix[1],"ELEVATION_"] <- ne.df[ix[1],"ELEV_FT"] <- 1192 ne.m[ix[2],"ELEVATION_"] <- ne.m[ix[2],"ELEV_FT"] <- 1141 ne.df[ix[2],"ELEVATION_"] <- ne.df[ix[2],"ELEV_FT"] <- 1141 ## ----re-run-ols----------------------------------------------------- m.ols <- lm(ANN_GDD50 ~ sqrt(ELEVATION_) + N, data=ne.df) summary(m.ols) ## ----ols-model-dx2, out.width='\\textwidth', w=12, h=4-------------- par(mfrow=c(1,3)) plot(m.ols, which=c(1,2,5)) par(mfrow=c(1,1)) ## ----plot-lm-fits-2, fig.cap = "Actual vs. predicted, OLS model with corrected observations"---- plot(ne.m$ANN_GDD50 ~ fitted(m.ols), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by OLS (corrected data)", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----ols-resid-summary---------------------------------------------- summary(residuals(m.ols)) ## ----add-east, eval=FALSE, echo=FALSE, results="hide", fig.keep='none'---- ## ## this code would be useful to answer the challenge ## m.ols.ts1 <- lm(ANN_GDD50 ~ N + E, data=ne.df) ## m.ols.ts1.z <- lm(ANN_GDD50 ~ N + E + sqrt(ELEVATION_), data=ne.df) ## m.ols.ts1.z.i <- lm(ANN_GDD50 ~ (N + E)*sqrt(ELEVATION_), data=ne.df) ## summary(m.ols.ts1) ## summary(m.ols) ## summary(m.ols.ts1.z) ## summary(m.ols.ts1.z.i) ## anova(m.ols.ts1.z.i, m.ols.ts1.z, m.ols.ts1) ## anova(m.ols.ts1.z, m.ols) ## AIC(m.ols.ts1.z.i, m.ols.ts1.z, m.ols.ts1, m.ols) ## par(mfrow=c(1,3)) ## plot(m.ols.ts1.z.i, which=c(1,2,5)) ## par(mfrow=c(1,1)) ## ix <- which(abs(residuals(m.ols.ts1.z)) > 400) ## ne.m$ols.ts1.z.resid <- ne.df$ols.ts1.z.resid <- residuals(m.ols.ts1.z) ## ne.df[ix,c("STATION_NA","STATE","ELEVATION_","N","E", "ols.ts1.z.resid")] ## ne.df$diff.resid.ols.ts1.z <- (ne.m$ols.ts1.z.resid - ne.m$ols.resid) ## ----add-OLS-resid-df----------------------------------------------- ne.m$ols.resid <- residuals(m.ols) ## ----bubble-function, out.width='0.75\\linewidth'------------------- bubble.sf <- function(.point.obj.name, .field.name, .field.label, .title="") { # make a plus/minus indicator eval(parse( text = paste0("pm <- factor(", .point.obj.name, "$", .field.name, "> 0)") )) # rename them levels(pm) <- c("-, overprediction", "+, underprediction") # plot eval(parse( text = paste0("ggplot(", .point.obj.name, ") + geom_sf(aes(colour = pm, size = abs(", .field.name, ")), shape = 1) + labs(size = paste('+/-', .field.label), colour = '', title = '", .title, "') + scale_colour_manual(values = c('red', 'green'))") )) } ## ----bubble-naive--------------------------------------------------- bubble.sf("ne.m", "ols.resid", "residuals, GDD50F", "OLS") ## ----load-gstat----------------------------------------------------- require(gstat) ## ----vgm-ols-resid, out.width='0.5\\textwidth'---------------------- v.r.ols <- variogram(ols.resid ~ 1, locations=ne.m, cutoff=100000, width=16000) plot(v.r.ols, pl=T) ## ----vgm-ols-resid-model, out.width='0.5\\textwidth'---------------- (vmf.r.ols <- fit.variogram(v.r.ols, vgm(15000, "Exp", 20000, 20000))) plot(v.r.ols, pl=T, model=vmf.r.ols) ## ----vgm-ols-cloud, out.width='0.5\\textwidth'---------------------- vc <- variogram(ols.resid ~ 1, locations=ne.m, cutoff=12000, cloud=T) plot(vc, pch=20, cex=2) ## ----why-ols-resid-diff--------------------------------------------- vc.df <- as.data.frame(vc) vc.close <- subset(vc.df, vc.df$dist < 8000) ## sort by separation, look for anomalies. vc.close[order(vc.close$dist),c("dist","gamma","left","right")] ## ------------------------------------------------------------------- print(ne.m[c(106,107),c("STATE","STATION_NA","LATITUDE_D","LONGITUDE_", "ELEVATION_", "ANN_GDD50","ols.resid")]) ## ----child-data, child='MappingRegionalClimate_GLS.Rnw'------------- ## ----GLS-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('MappingRegionalClimate.Rnw') ## ----fit-reml------------------------------------------------------- require(nlme) vmf.r.ols[1:2,] m.gls <- gls(model=ANN_GDD50 ~ sqrt(ELEVATION_) + N, data=ne.df, correlation=corExp( value=c(vmf.r.ols[2,"range"]), form=~E + N, nugget=FALSE)) ## ----summary-reml--------------------------------------------------- summary(m.gls) ## ----intervals-reml------------------------------------------------- intervals(m.gls, level=0.95)$coef ## ----summary-intervals-coef----------------------------------------- intervals(m.gls, level=0.95)$corStruct ## ----compare-gls-ols-corr------------------------------------------- intervals(m.gls, level=0.95)$corStruct vmf.r.ols[2,"range"] ## ----plot-gls-fits-------------------------------------------------- plot(ne.m$ANN_GDD50 ~ predict(m.gls), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by GLS", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----fit-reml-resid------------------------------------------------- ne.m$gls.resid <- residuals(m.gls) summary(ne.m$gls.resid) bubble.sf("ne.m", "gls.resid", "residuals, GDD50F", "GLS") ## ----gls-resid-vgm, out.width='0.5\\textwidth'---------------------- v.r.gls <- variogram(gls.resid ~ 1, ne.m, cutoff=120000, width=12000) (vmf.r.gls <- vgm(psill=var(ne.m$gls.resid), model="Exp", range=intervals(m.gls)$corStruct["range","est."], nugget=TRUE)) plot(v.r.gls, pl=T, model=vmf.r.gls) ## ----show-spherical------------------------------------------------- plot(v.r.ols, pl=T) ## ----fit-spherical-------------------------------------------------- (vmf.r.ols.sph <- fit.variogram(v.r.ols, vgm(psill=20000, model="Sph", range=40000, nugget=15000))) plot(v.r.ols, pl=T, model=vmf.r.ols.sph) ## ----compute-prop-nugget-------------------------------------------- (prop.nugget <- vmf.r.ols.sph[1,"psill"]/sum(vmf.r.ols.sph[,"psill"])) ## ----fit-reml-with-nugget------------------------------------------- m.gls.2 <- gls(model=ANN_GDD50 ~ sqrt(ELEVATION_) + N, data=ne.df, correlation=corSpher( value=c(vmf.r.ols.sph[2,"range"], prop.nugget), form=~E + N, nugget=TRUE)) ## ----compare-corr-structs------------------------------------------- intervals(m.gls, level=0.95)$corStruct intervals(m.gls.2, level=0.95)$corStruct ## ----compare-ranges------------------------------------------------- intervals(m.gls)$corStruct["range", 2]*3 intervals(m.gls.2)$corStruct["range", 2] ## ----compare-regr-coeff--------------------------------------------- intervals(m.gls, level=0.95)$coef[,2] intervals(m.gls.2, level=0.95)$coef[,2] ## ----compare-coeff-------------------------------------------------- coefficients(m.gls) coefficients(m.ols) round(100*(coefficients(m.gls) - coefficients(m.ols))/coefficients(m.ols),2) ## ----bubble-diffs--------------------------------------------------- ne.m$diff.gls.ols.resid <- (ne.m$gls.resid - ne.m$ols.resid) summary(ne.m$diff.gls.ols.resid) bubble.sf("ne.m", "diff.gls.ols.resid", "residual difference", "GLS - OLS") ## ----scatter-diff, out.width='\\linewidth', w=10, h=6--------------- ggplot() + geom_point(aes(x=sqrt(ne.m$ELEVATION_), y=ne.m$diff.gls.ols.resid, size=ne.m$ANN_GDD50, colour=ne.m$STATE)) + xlab("Sqrt(Elevation)") + ylab("GLS - OLS residual") + geom_hline(yintercept=0, linewidth=1.2) + labs(colour = "State", size = "Annual GDD50F") ## ----big-diff-resid------------------------------------------------- ix <- which.max(ne.m$diff.gls.ols.resid) ne.df[ix,c("STATION_NA","STATE","ELEVATION_","N","E")] ix <- which.min(ne.m$diff.gls.ols.resid) ne.df[ix,c("STATION_NA","STATE","ELEVATION_","N","E")] ## ----child-data, child='MappingRegionalClimate_pred_OLS_GLS.Rnw'---- ## ----pred-OLS-GLS-set-parent, echo=FALSE, cache=FALSE--------------- set_parent('MappingRegionalClimate.Rnw') ## ----predict-4km---------------------------------------------------- dem.ne.m.df$pred.ols <- predict(m.ols, newdata=dem.ne.m.df) dem.ne.m.df$pred.gls <- predict(m.gls, newdata=dem.ne.m.df) summary(dem.ne.m.df[,-(1:3)]) ## ----define-display-map-function------------------------------------ display.prediction.map <- function(.prediction, .plot.title, .legend.title, .plot.limits=NULL, .palette="YlGnBu") { ggplot(data=dem.ne.m.df) + # note must find the indirect name in the dataframe geom_point(aes(x=E, y=N, color=.data[[as.name(.prediction)]])) + xlab("E") + ylab("N") + coord_fixed() + geom_point(aes(x=E, y=N, color=ANN_GDD50), data=ne.df, shape=I(20)) + ggtitle(.plot.title) + scale_colour_distiller(name="GDD50", space="Lab", palette=.palette, limits=.plot.limits) } ## ----plot-4km, out.width='.8\\textwidth'---------------------------- ols.ix <- which(names(dem.ne.m.df)=="pred.ols") gls.ix <- which(names(dem.ne.m.df)=="pred.gls") # set up limits of the scale for the target variable # use the extremes +/-10 GDD, to make sure that all values are shown (gdd.pred.lim <- round( c(min(dem.ne.m.df[,c(ols.ix, gls.ix)])-10, max(dem.ne.m.df[,ols.ix, gls.ix])+10))) display.prediction.map("pred.ols", "Annual GDD, base 50F, OLS prediction", "GDD50", gdd.pred.lim) display.prediction.map("pred.gls", "Annual GDD, base 50F, GLS prediction", "GDD50", gdd.pred.lim) ## ----diffs-4km-gls-ols---------------------------------------------- summary(dem.ne.m.df$diff.gls.ols <- dem.ne.m.df$pred.gls - dem.ne.m.df$pred.ols) ## ----define-difference-map-function--------------------------------- display.difference.map <- function(.diff.map.name, .diff.map.title, .legend.name, .palette="Spectral") { ggplot(data=dem.ne.m.df) + geom_point(aes(x = .data[["E"]], y = .data[["N"]], color=.data[[as.name(.diff.map.name)]])) + xlab("E") + ylab("N") + coord_fixed() + ggtitle(.diff.map.title) + scale_colour_distiller(name=.legend.name, space="Lab", palette=.palette) } ## ----display-diffs-4km-gls-ols, out.width='.8\\textwidth'----------- display.difference.map("diff.gls.ols", "Annual GDD, base 50F, GLS - OLS predictions", "+/- GDD50") ## ----child-data, child='MappingRegionalClimate_RK.Rnw'-------------- ## ----RKGLS-set-parent, echo=FALSE, cache=FALSE---------------------- set_parent('MappingRegionalClimate.Rnw') ## ----print-gls-vmf-------------------------------------------------- print(vmf.r.gls) ## ----gls-resid-ok-grid, out.width='.6\\textwidth'------------------- summary(ne.m$gls.resid) hist(ne.m$gls.resid, freq=F, xlab="GDD50", main="GLS residuals") rug(ne.m$gls.resid) ## ----krige-ok-gls-resid-1------------------------------------------- class(dem.ne.m.sf) system.time( ok.gls.resid <- krige(gls.resid ~ 1, loc=ne.m, newdata=dem.ne.m.sf, model=vmf.r.gls) ) summary(ok.gls.resid) hist(ok.gls.resid$var1.pred, main = "OK deviations from GLS trend surface", xlab = "GDD50") abline(v=0, col="red") rug(ok.gls.resid$var1.pred) ## ----krige-ok-gls-resid-2------------------------------------------- dem.ne.m.df$ok.gls.resid <- ok.gls.resid$var1.pred dem.ne.m.df$ok.gls.resid.var <- ok.gls.resid$var1.var ## ----krige-ok-gls-resid, out.width='.9\\textwidth'------------------ ggplot() + geom_point(aes(x=E, y=N, colour=ok.gls.resid), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Residuals from GLS trend surface, GDD base 50F") + scale_colour_distiller(name="GDD50", space="Lab", palette="RdBu") ## ------------------------------------------------------------------- mean(dem.ne.m.df$ok.gls.resid) # spatial mean over the grid mean(ne.m$gls.resid) # arithmetic mean at observation points ## ----show-ok-gls-resid-var, out.width='.9\\textwidth'--------------- ggplot() + geom_point(aes(x=E, y=N, colour=sqrt(ok.gls.resid.var)), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + geom_point(aes(x=E, y=N), data=ne.df, size=0.5, colour="black", shape=I(20)) + ggtitle("Residuals from GLS trend, kriging prediction standard deviation") + scale_colour_distiller(name="GDD50", space="Lab", palette="BuPu", trans="reverse") ## ----rk-gls, out.width='.9\\textwidth'------------------------------ summary(dem.ne.m.df$pred.rkgls <- dem.ne.m.df$pred.gls + dem.ne.m.df$ok.gls.resid) display.prediction.map("pred.rkgls", "Annual GDD, base 50F, GLS-RK prediction", "GDD50") ## ----mask-predictions, out.width='.9\\textwidth', fig.caption="Limited to the study area"---- require(terra) str(dem.ne.m.df) # by default the first two fields are taken as the coordinates # get the CRS from the plygon object that will be the mask pred.ne.m.rast <- rast(dem.ne.m.df, crs=st_crs(state.ne.m)$proj4string) print(pred.ne.m.rast) pred.ne.m.rast <- mask(pred.ne.m.rast, state.ne.m) terra::plot(pred.ne.m.rast) ## ----mask-predictions-1, out.width='.9\\textwidth', fig.caption="RK-GLS prediction"---- require(RColorBrewer) terra::plot(pred.ne.m.rast, y = "pred.rkgls", col = colorRampPalette( brewer.pal(n = 9, name = "YlGnBu"))(64), main = "RK-GLS prediction, annual GDD50F") ## ----child-data, child='MappingRegionalClimate_KED.Rnw'------------- ## ----KED-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('MappingRegionalClimate.Rnw') ## ----add-N-to-sf-objects-------------------------------------------- ne.m$N <- st_coordinates(ne.m)[, 2] dem.ne.m.sf$N <- st_coordinates(dem.ne.m.sf)[, 2] ## ----ked-vgm, out.width='0.6\\textwidth', w=8, h=6------------------ v.ked <- variogram(ANN_GDD50 ~ sqrt(ELEVATION_) + N, locations=ne.m, cutoff=100000, width=16000) plot(v.ked, pl=T) ## ----ked-vgm-fit, out.width='0.6\\textwidth', w=8, h=6-------------- (vmf.ked <- fit.variogram(v.ked, vgm(15000, "Exp", 20000, 20000))) plot(v.ked, plot.numbers=TRUE, model=vmf.ked) ## ----ked-pred-grid-------------------------------------------------- k.ked <- krige(ANN_GDD50 ~ sqrt(ELEVATION_)+ N, locations=ne.m, newdata=dem.ne.m.sf, model=vmf.ked) summary(k.ked) ## ----ked-show, out.width='.9\\textwidth'---------------------------- dem.ne.m.df$pred.ked <- k.ked$var1.pred display.prediction.map("pred.ked", "Annual GDD, base 50F, KED prediction", "GDD50") ## ----diffs-4km-rkgls-ked-------------------------------------------- summary(dem.ne.m.df$diff.rkgls.ked <- dem.ne.m.df$pred.rkgls - dem.ne.m.df$pred.ked) ## ----display-diffs-4km-rkgls-ked, out.width='.8\\textwidth'--------- display.difference.map("diff.rkgls.ked", "Difference annual GDD base 50F, RK-GLS - KED", "+/- GDD50") ## ----mask-predictions-ked------------------------------------------- # by default the first two fields are taken as the coordinates pred.ked.rast <- rast(dem.ne.m.df, crs=st_crs(state.ne.m)$proj4string)["pred.ked"] pred.ked.rast <- mask(pred.ked.rast, state.ne.m) terra::plot(pred.ked.rast, y = "pred.ked", col = colorRampPalette( brewer.pal(n = 9, name = "YlGnBu"))(64), main = "KED prediction [N, sqrt(ELEVATION)], annual GDD50F") ## ----krige-loocv-ked, results='hide'-------------------------------- kcv.ked <- krige.cv(ANN_GDD50 ~ sqrt(ELEVATION_)+ N, locations=ne.m, model=vmf.ked) ## ----krige-loocv-ked-2---------------------------------------------- summary(kcv.ked$residual) ## ------------------------------------------------------------------- (loocv.ked.rmse <- sqrt(sum(kcv.ked$residual^2)/length(kcv.ked$residual))) ## ----ok-loocv-bubble-ked-------------------------------------------- bubble.sf("kcv.ked", "residual", "GDD50", .title="LOOCV KED residuals") ## ----ked-local------------------------------------------------------ ked.n.max <- 61 # change this to try other numbers of neighbours k.ked.nn <- krige(ANN_GDD50 ~ sqrt(ELEVATION_)+ N, locations=ne.m, newdata=dem.ne.m.sf, model=vmf.ked, nmax=ked.n.max) summary(k.ked.nn) ## ----ked-local-display-map------------------------------------------ dem.ne.m.df$pred.ked.nn <- k.ked.nn$var1.pred dem.ne.m.df$pred.ked.nn.sd <- sqrt(k.ked.nn$var1.var) display.prediction.map("pred.ked.nn", paste("Annual GDD, base 50F, KED prediction,", ked.n.max, "nearest neighbours"), "GDD50") ## ----ked-compare-global-local--------------------------------------- summary(dem.ne.m.df$diff.ked <- dem.ne.m.df$pred.ked - dem.ne.m.df$pred.ked.nn) display.difference.map("diff.ked", paste("Difference annual GDD base 50F, KED global - KED", ked.n.max,"nearest neighbours"), "+/- GDD50") ## ----ked-local-cv, results='hide'----------------------------------- kcv.ked.nn <- krige.cv(ANN_GDD50 ~ sqrt(ELEVATION_)+ N, locations=ne.m, model=vmf.ked, nmax=ked.n.max) ## ----ked-local-cv-2------------------------------------------------- summary(kcv.ked.nn$residual) summary(kcv.ked$residual) (loocv.ked.nn.rmse <- sqrt(sum(kcv.ked.nn$residual^2)/length(kcv.ked.nn$residual))) (loocv.ked.rmse <- sqrt(sum(kcv.ked$residual^2)/length(kcv.ked$residual))) ## ----ked-global-local-cv-bubble------------------------------------- summary(kcv.ked$diff <- kcv.ked$residual - kcv.ked.nn$residual) bubble.sf("kcv.ked", "diff", "delta GDD50", .title=paste("KED - KED", ked.n.max, "nearest neighbour residuals")) ## ------------------------------------------------------------------- k.ok <- krige(ANN_GDD50 ~ sqrt(ELEVATION_)+ N, locations=ne.m, newdata=dem.ne.m.sf, model=NULL) ## ------------------------------------------------------------------- k.okr <- krige(ols.resid ~ 1, locations=ne.m, newdata=dem.ne.m.sf, model=vmf.r.ols) ## ------------------------------------------------------------------- k.ok$rk.pred <- k.ok$var1.pred + k.okr$var1.pred ## ------------------------------------------------------------------- k.ok$diff.pred <- k.ok$rk.pred - k.ked$var1.pred summary(k.ok$rk.pred); summary(k.ked$var1.pred) summary(k.ok$diff.pred) ## ----display-okrk-ked-diff, out.width='.8\\textwidth'--------------- dem.ne.m.df$diff.okrk.ked.pred <- k.ok$diff.pred display.difference.map("diff.okrk.ked.pred", "Naive RK vs.\ KED surface, GDD base 50F", "+/- GDD50") ## ----child-data, child='MappingRegionalClimate_GAM.Rnw'------------- ## ----GAM-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('MappingRegionalClimate.Rnw') ## ----plot-marginal-loess, out.width='0.9\\linewidth', fig.caption = "Relation of target variable with predictors"---- g1 <- ggplot(ne.df, aes(x=E, y=ANN_GDD50)) + geom_point() + geom_smooth(method="loess") g2 <- ggplot(ne.df, aes(x=N, y=ANN_GDD50)) + geom_point() + geom_smooth(method="loess") g3 <- ggplot(ne.df, aes(x=ELEVATION_, y=ANN_GDD50)) + geom_point() + geom_smooth(method="loess") g4 <- ggplot(ne.df, aes(x=sqrt(ELEVATION_), y=ANN_GDD50)) + geom_point() + geom_smooth(method="loess") require(gridExtra) grid.arrange(g1, g2, g3, g4, ncol = 2) ## ----load-mgcv------------------------------------------------------ require(mgcv) ## ----fit-gam-------------------------------------------------------- m.g.xy <- gam(ANN_GDD50 ~ s(E, N) + s(ELEVATION_), data=ne.df) summary(m.g.xy) summary(residuals(m.g.xy)) ## ----gam-resid-bubble----------------------------------------------- ne.m$resid.m.g.xy <- residuals(m.g.xy) bubble.sf("ne.m", "resid.m.g.xy", "GDD50", "Residuals from GAM") ## ----gam-resid-vgm, out.width='.55\\linewidth'---------------------- vr <- variogram(resid.m.g.xy ~ 1, loc=ne.m, cutoff=50000, width=5000) plot(vr, pl=T) ## ----plot-gam-2D, out.width='\\linewidth'--------------------------- plot.gam(m.g.xy, rug=T, se=T, select=1, scheme=1, theta=30+130, phi=30) ## ----plot-gam-2D-se, out.width='\\linewidth'------------------------ vis.gam(m.g.xy, plot.type="persp", color="terrain", theta=160, zlab="Annual GDD50", se=1.96) ## ----plot-gam-1D-elev----------------------------------------------- plot.gam(m.g.xy, select=2, rug=T, se=T, residuals=T, pch=20, shade=T, seWithMean=T, shade.col="lightblue") ## ----plot-gam-fits-------------------------------------------------- (rmse.gam <- sqrt(sum(residuals(m.g.xy)^2)/length(residuals(m.g.xy)))) plot(ne.m$ANN_GDD50 ~ predict(m.g.xy, newdata=ne.df), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.df$STATE), pch=20, col=1:4) grid(); abline(0,1) ## ----predict-gam, out.width='0.8\\textwidth'------------------------ tmp <- predict.gam(object=m.g.xy, newdata=dem.ne.m.df, se.fit=TRUE) summary(tmp$fit) summary(tmp$se.fit) dem.ne.m.df$pred.gam <- tmp$fit dem.ne.m.df$pred.gam.se <- tmp$se.fit display.prediction.map("pred.gam", "Annual GDD, base 50F, GAM prediction", "GDD50") ## ----predict-gam-se, out.width='0.8\\textwidth'--------------------- ggplot() + geom_point(aes(x=E, y=N, colour=pred.gam.se), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Annual GDD base 50F, Standard error of GAM prediction") + scale_colour_distiller(name="GDD50 s.e.", space="Lab", palette="RdYlGn", direction=-1) ## ----diffs-4km-gls-gam---------------------------------------------- summary(dem.ne.m.df$diff.gls.gam <- dem.ne.m.df$pred.gls - dem.ne.m.df$pred.gam) ## ----display-diffs-4km-gls-gam, out.width='.8\\textwidth'----------- display.difference.map("diff.gls.gam", "Difference annual GDD base 50F, GLS - GAM", "+/- GDD50") ## ----diffs-4km-rkgls-gam-------------------------------------------- dem.ne.m.df$diff.rkgls.gam <- dem.ne.m.df$pred.rkgls - dem.ne.m.df$pred.gam summary(dem.ne.m.df$diff.rkgls.gam) ## ----display-diffs-4km-rkgls-gam, out.width='.8\\textwidth'--------- ggplot() + geom_point(aes(x=E, y=N, colour=diff.rkgls.gam), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Difference annual GDD base 50F, GLS/RK - GAM") + scale_colour_distiller(name="GDD50", space="Lab", palette="Spectral") ## ----child-data, child='MappingRegionalClimate_RF.Rnw'-------------- ## ----RF-set-parent, echo=FALSE, cache=FALSE------------------------- set_parent('MappingRegionalClimate.Rnw') ## ----load-rpart----------------------------------------------------- require(rpart) ## ----set-seed-rt, echo=F, results='hide'---------------------------- # set seed for the random number generator to get a consistent result for the document # this should *not* be done when following the exercise set.seed(318316) ## ----trees-comp-rpart----------------------------------------------- m.rt <- rpart(ANN_GDD50 ~ N + E + ELEVATION_, data=ne.df, minsplit=2, cp=0.003) ## ----summary-rt----------------------------------------------------- print(m.rt) ## ----trees-plot-unpruned, out.width='\\textwidth'------------------- require(rpart.plot) rpart.plot(m.rt, digits=3, type=4, extra=1) ## ----trees-import-rpart--------------------------------------------- x <- m.rt$variable.importance data.frame(variableImportance = 100 * x / sum(x)) ## ----trees-plotcp-rp, out.width='\\textwidth'----------------------- printcp(m.rt) plotcp(m.rt) ## ----trees-prune-rpart---------------------------------------------- (m.rt.p <- prune(m.rt, cp=0.0045)) ## ----trees-plot-pruned, out.width='\\textwidth'--------------------- rpart.plot(m.rt.p, digits=3, type=4, extra=1) ## ----trees-predict-cal-pts------------------------------------------ summary(p.rt.p <- predict(m.rt.p, newdata=ne.df)) ## ----trees-count-unique-pred---------------------------------------- length(unique(p.rt.p)) ## ----trees-residuals-histo------------------------------------------ summary(residuals.rt.p <- ne.df$ANN_GDD50 - p.rt.p) hist(residuals.rt.p, main="Residuals from regression tree fit", xlab="ANN_GDD50") rug(residuals.rt.p) sqrt(mean(residuals.rt.p^2)/length(residuals.rt.p)) ## ----trees-predict-rpart, out.width='0.75\\textwidth'--------------- plot(ne.df$ANN_GDD50 ~ p.rt.p, asp=1, pch=20, xlab="fitted by regression tree", ylab="actual", xlim=c(500,4200), ylim=c(500,4200), col=ne.df$STATE, main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----pred-rt-area--------------------------------------------------- dem.ne.m.df$pred.rt <- predict(m.rt.p, newdata=dem.ne.m.df) summary(dem.ne.m.df$pred.rt) ## ----plot-rt-area--------------------------------------------------- display.prediction.map("pred.rt", "Annual GDD, base 50F, regression tree prediction", "GDD50") ## ------------------------------------------------------------------- require(ranger) ## ----set-seed-rf, echo=F, results='hide'---------------------------- # set seed for the random number generator to get a consistent result for the document # this should *not* be done when following the exercise set.seed(318316) ## ----random-forest-------------------------------------------------- m.rf <- ranger(ANN_GDD50 ~ ELEVATION_ + N + E, data=ne.df, num.trees=1200, importance="permutation") # proportional importance ranger::importance(m.rf)/sum(ranger::importance(m.rf))*100 ranger::importance(m.rf)/dim(ne.df)[1] ## ----plot-rf-fits--------------------------------------------------- summary(rf.fits <- predict(m.rf, data = ne.df)$predictions) plot(ne.m$ANN_GDD50 ~ rf.fits, col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----rf-fits-errors------------------------------------------------- summary(rf.resid <- ne.m$ANN_GDD50 - rf.fits) (rf.me <- mean(rf.resid)) (rf.rmse <- sqrt(mean(rf.resid^2)/length(rf.resid))) ## ----plot-rf-oob---------------------------------------------------- summary(rf.oob <- m.rf$predictions) plot(ne.m$ANN_GDD50 ~ rf.oob, col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest (OOB)", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid() abline(0,1) ## ----rf-oob-errors-------------------------------------------------- summary(rf.oob.resid <- ne.m$ANN_GDD50 - rf.oob) (rf.oob.me <- mean(rf.oob.resid)) (rf.oob.rmse <- sqrt(mean(rf.oob.resid^2)/length(rf.oob.resid))) ## ----rf-compare-errors---------------------------------------------- (rf.oob.me/rf.me) (rf.oob.rmse/rf.rmse) ## ----rf-resid------------------------------------------------------- ne.m$rf.resid <- (ne.df$ANN_GDD50 - rf.fits) summary(ne.m$rf.resid) bubble.sf("ne.m", "rf.resid", "GDD50", "Random Forest fitted residuals, actual-predicted") ## ----rf-resid-oob--------------------------------------------------- ne.m$rf.resid.oob <- (ne.df$ANN_GDD50 - rf.oob) summary(ne.m$rf.resid.oob) bubble.sf("ne.m", "rf.resid.oob", "GDD50", "Random Forest OOB residuals, actual-predicted") ## ------------------------------------------------------------------- summary(ne.m$rf.resid.oob/ne.m$rf.resid) ## ----which-rf-worst------------------------------------------------- (ix <- order(abs(ne.m$rf.resid.oob), decreasing=TRUE)[1:8]) ne.m[ix,c("STATE","STATION_NA","ELEVATION_", "gls.resid", "rf.resid")] names(ne.m) ## ----vgm-rf-resid, out.width='0.5\\textwidth'----------------------- v.rf <- variogram(rf.resid ~ 1, locations=ne.m, cutoff=100000, width=16000) plot(v.rf, pl=T) ## ------------------------------------------------------------------- summary(ne.m$ols.resid) summary(ne.m$rf.resid) summary(ne.m$gls.resid) sd(ne.m$ols.resid) sd(ne.m$rf.resid) sd(ne.m$gls.resid) ## ----predict-4km-rf------------------------------------------------- dem.ne.m.df$pred.rf <- predict(m.rf, data=dem.ne.m.df)$prediction summary(dem.ne.m.df$pred.rf) ## ----display-4km-rf, out.width='.8\\textwidth'---------------------- display.prediction.map("pred.rf", "Annual GDD, base 50F, random forest prediction", "GDD50") ## ----diffs-4km-gls-rf----------------------------------------------- summary(dem.ne.m.df$diff.gls.rf <- dem.ne.m.df$pred.gls - dem.ne.m.df$pred.rf) ## ----display-diffs-4km-gls-rf, out.width='.8\\textwidth'------------ display.difference.map("diff.gls.rf", "Annual GDD, base 50F, GLS - RF predictions", "+/- GDD50") ## ----load-caret----------------------------------------------------- require(caret) ## ----list-models---------------------------------------------------- length(names(getModelInfo())) head(names(getModelInfo()),24) ## ----caret-model-info----------------------------------------------- getModelInfo("ranger")$ranger$parameters ## ----set-seed-caret, echo=F, results='hide'------------------------- # set seed for the random number generator to get a consistent result for the document # this should *not* be done when following the exercise set.seed(318316) ## ----caret---------------------------------------------------------- dim(preds <- ne.df[, c("E", "N", "ELEVATION_")]) length(response <- ne.df[, "ANN_GDD50"]) system.time( ranger.tune <- train(x = preds, y = response, method="ranger", tuneGrid = expand.grid(.mtry = 1:3, .splitrule = "variance", .min.node.size = 1:10), trControl = trainControl(method = 'cv')) ) ## ----caret-show-results, out.width='.65\\textwidth', fig.width=8, fig.height=5---- print(ranger.tune) names(ranger.tune$result) ix <- which.min(ranger.tune$result$RMSE) ranger.tune$result[ix, c(1,3,4)] ix <- which.max(ranger.tune$result$Rsquared) ranger.tune$result[ix, c(1,3,5)] ix <- which.min(ranger.tune$result$MAE) ranger.tune$result[ix, c(1,3,6)] plot.train(ranger.tune, metric="RMSE") plot.train(ranger.tune, metric="Rsquared") plot.train(ranger.tune, metric="MAE") ## ----fit-ranger----------------------------------------------------- (ranger.rf <- ranger(ANN_GDD50 ~ N + E + ELEVATION_, data=ne.df, mtry=2, min.node.size=5)) ## ----plot-ranger-fits----------------------------------------------- summary(ranger.fits <- predict(ranger.rf, data=ne.df)$predictions) plot(ne.df$ANN_GDD50 ~ ranger.fits, col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by ranger", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) ## ----predict-4km-ranger--------------------------------------------- summary(dem.ne.m.df$pred.ranger <- predict(ranger.rf, data=dem.ne.m.df)$predictions) ## ----display-4km-ranger, out.width='.8\\textwidth'------------------ display.prediction.map("pred.ranger", "Annual GDD, base 50F, ranger RF prediction", "GDD50") ## ----load-cubist---------------------------------------------------- library(Cubist) ## ----vignette-cubist, eval=FALSE------------------------------------ ## vignette("cubist") ## ----set-seed-cubist, echo=F, results='hide'------------------------ # set seed for the random number generator to get a consistent result for the document # this should *not* be done when following the exercise set.seed(316318) ## ----cubist-train--------------------------------------------------- all.preds <- ne.df[, c("N", "E", "ELEVATION_")] all.resp <- ne.df[ , "ANN_GDD50"] system.time( cubist.tune <- train(x = all.preds, y = all.resp, "cubist", tuneGrid = expand.grid(.committees = 1:12, .neighbors = 0:8), trControl = trainControl(method = 'cv')) ) ## ----cubist-train-results, fig.width=8, fig.height=5---------------- print(cubist.tune) plot(cubist.tune, metric="RMSE") plot(cubist.tune, metric="Rsquared") plot(cubist.tune, metric="MAE") ## ----cubist-fit-best------------------------------------------------ require(Cubist) c.model <- cubist(x = all.preds, y = all.resp, committees=6) summary(c.model) ## ----cubist-model-eval-fit------------------------------------------ # predictive accuracy cubist.fits <- predict(c.model, newdata=all.preds, neighbors=cubist.tune$bestTune$neighbors) ## Test set RMSE sqrt(mean((cubist.fits - all.resp)^2)) cor(cubist.fits, all.resp)^2 # R^2 plot(ne.df$ANN_GDD50 ~ cubist.fits, col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by cubist", ylab="Actual", main="Annual GDD50") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) ## ----cubist-model-predict------------------------------------------- summary(dem.ne.m.df$pred.cubist <- predict(c.model, newdata=dem.ne.m.df, neighbours=cubist.tune$bestTune$neighbors)) ## ----display-4km-cubist, out.width='.8\\textwidth'------------------ display.prediction.map("pred.cubist", "Annual GDD, base 50F, Cubist prediction", "GDD50") ## ----build-extended-matrix------------------------------------------ pred.field.names <- c("ELEVATION_", "E", "N", "dist.lakes", "dist.coast", "mrvbf", "tri3", "pop15", "pop2pt5") pred.field.names.base <- c("ELEVATION_", "E", "N") (model.formula <- paste("ANN_GDD50 ~", paste0(pred.field.names, collapse=" + "))) (model.formula.base <- paste("ANN_GDD50 ~", paste0(pred.field.names.base, collapse=" + "))) names(dem.ne.m.df) ## ----pairs,fig.height=10, fig.width=10------------------------------ cor(ne.df[, pred.field.names]) corrplot::corrplot(cor(ne.df[, pred.field.names]), diag = FALSE, type = "upper", method = "ellipse", addCoef.col = "black") ## ----pca,fig.height=7, fig.width=7---------------------------------- pc <- prcomp(ne.df[, pred.field.names], scale. = TRUE, retx=TRUE) summary(pc) pc$rotation biplot(pc) biplot(pc, choices=3:4) ## ----rf-extended---------------------------------------------------- dim(preds <- ne.df[, pred.field.names]) length(response <- ne.df[, "ANN_GDD50"]) system.time( ranger.tune <- train(x = preds, y = response, method="ranger", tuneGrid = expand.grid(.mtry = 2:7, .splitrule = "variance", .min.node.size = 1:10), trControl = trainControl(method = 'cv')) ) ## ----ranger.tune.results, fig.width=6, fig.height=4----------------- ix <- which.min(ranger.tune$result$RMSE) ranger.tune$result[ix, c(1,3,4)] ix <- which.max(ranger.tune$result$Rsquared) ranger.tune$result[ix, c(1,3,5)] ix <- which.min(ranger.tune$result$MAE) ranger.tune$result[ix, c(1,3,6)] plot.train(ranger.tune, metric="RMSE") plot.train(ranger.tune, metric="Rsquared") ## ----rf.full-------------------------------------------------------- rf.ext <- ranger(model.formula, data=ne.df, importance="impurity", mtry=4, min.node.size=3, oob.error=TRUE, num.trees=1024) print(rf.ext) str(rf.ext, max.level=1) round(rf.ext$prediction.error/rf.ext$num.samples) plot(rf.ext$predictions, ne.df$ANN_GDD50, asp=1, col=ne.df$STATE, pch=20, main="ANN_GDD50", ylab="Measured", xlab="Ranger RF fit") abline(0,1); grid() round(ranger::importance(rf.ext)/sum(ranger::importance(rf.ext)),3) ## ----rf.map, eval=TRUE, fig.width=6, fig.height=6------------------- dem.ne.m.df$pred.rf.ext <- predict(rf.ext, data=dem.ne.m.df)$predictions summary(dem.ne.m.df$pred.rf.ext) display.prediction.map("pred.rf.ext", "Annual GDD, base 50F, ranger prediction, 9 predictors", "GDD50") ## ----rf.ext.min.depth.distribution, fig.width=8, fig.height=10, out.width="0.9\\linewidth"---- require(randomForestExplainer) tmp <- min_depth_distribution(rf.ext) plot_min_depth_distribution(tmp) ## ----Shapley-predictor---------------------------------------------- require(iml) # matrix of predictors to be evaluated X <- ne.df[, pred.field.names] # the predictor <- Predictor$new(model = rf.ext, data = X, y = ne.df[, "ANN_GDD50"]) ## ----Shapley-values-Mansfield--------------------------------------- ix <- which.min(ne.df[, "ANN_GDD50"]) ne.df[ix, 2:3] X[ix,] shapley <- iml::Shapley$new(predictor, x.interest = X[ix, ]) shapley$plot() ## ----Shapley-values-Ithaca------------------------------------------ ix <- which(ne.df[, "STATION_NA"] == "PITTSBURGH INTL AP") X[ix,] shapley <- iml::Shapley$new(predictor, x.interest = X[ix, ]) shapley$plot() ## ----SHAP----------------------------------------------------------- require(fastshap) # a prediction function pfun <- function(object, newdata) { predict(object, data = newdata)$predictions } # matrix of predictors to be evaluated X <- ne.df[ , pred.field.names] fshap <- fastshap::explain(object = rf.ext, X = X, shap_only = FALSE, # also return feature and baseline values pred_wrapper = pfun, nsim = 24) names(fshap) head(fshap$shapley_values) ne.df[1, 2:3] ne.df[1, pred.field.names] ## ----shapviz-importance--------------------------------------------- require(shapviz) sv.fshap <- shapviz(fshap, X = ne.df[ , pred.field.names], bg_X = ne.df) # small dataset, can see all of them class(sv.fshap) sv_importance(sv.fshap, kind = "bee") ## ----shapviz-dependence--------------------------------------------- sv_dependence(sv.fshap, v = "dist.lakes", color_var = "N") ix <- which(fshap$shapley_values[, "dist.lakes"] > 350) ne.df[ix, 2:5] ## ----display.rf.diff, fig.width=6, fig.height=6--------------------- summary(dem.ne.m.df$diff.rf.9.3 <- dem.ne.m.df$pred.rf.ext - dem.ne.m.df$pred.rf) display.difference.map("diff.rf.9.3", "Difference RF 9 and RF 3 predictors", "+/- GDD50") ## ----child-data, child='MappingRegionalClimate_TPS.Rnw'------------- ## ----TPS-set-parent, echo=FALSE, cache=FALSE------------------------ set_parent('MappingRegionalClimate.Rnw') ## ----TPS------------------------------------------------------------ require(fields) ne.tps <- ne.df ne.tps$coords <- matrix(c(ne.df$E, ne.df$N), byrow=F, ncol=2) surf.1 <-Tps(ne.tps$coords, ne.tps$ANN_GDD50) class(surf.1) summary(surf.1) ## ----TPS-grid------------------------------------------------------- resolution <- 9 st_bbox(state.ne.m) (ne.sq.km <- diff(st_bbox(state.ne.m)[c(1,3)]) * diff(st_bbox(state.ne.m)[c(2,4)])/10^6) (approx.n.grid.cells <- ceiling(ne.sq.km/(resolution^2))) states.grid <- st_sample(state.ne.m, size = approx.n.grid.cells, type="regular") class(states.grid) states.grid.df <- as.data.frame(states.grid) states.grid.df$coords <- as.matrix(st_coordinates(states.grid)) str(states.grid.df) ## ----TPS-predict---------------------------------------------------- surf.1.pred <- predict.Krig(surf.1, states.grid.df$coords) class(surf.1.pred) dim(surf.1.pred) summary(as.vector(surf.1.pred)) ## ----TPS-map, out.width='\\textwidth'------------------------------- plot(states.grid.df$coords, pch=20, asp=1, cex=.6, col=sp::bpy.colors(256)[cut(surf.1.pred, 256)], xlab="E", ylab="N", main="Annual GDD50") plot(state.ne.m.boundary, add = TRUE, col = "darkgray") ## ----get-value-tps-------------------------------------------------- pt <- st_sfc(st_point(x=c(-76.402175, 42.453271))) st_crs(pt) <- 4326 pt <- st_transform(pt, st_crs(states.grid)) st_coordinates(pt) (pt.pred <- surf.1.pred <- predict.Krig(surf.1, st_coordinates(pt))) ## ----child-data, child='MappingRegionalClimate_OK.Rnw'-------------- ## ----OK-set-parent, echo=FALSE, cache=FALSE------------------------- set_parent('MappingRegionalClimate.Rnw') ## ----compute-cutoff, echo=FALSE------------------------------------- cutoff.default <- (sqrt(diff(st_bbox(ne.m)[c(1,3)])^2 + diff(st_bbox(ne.m)[c(2,4)])^2))/3 binwidth.default <- cutoff.default/15 ## ----display-ordinary-v, out.width='0.6\\textwidth', w=8, h=6------- library(gstat) st_bbox(ne.m) (cutoff.default <- (sqrt(diff(st_bbox(ne.m)[c(1,3)])^2 + diff(st_bbox(ne.m)[c(2,4)])^2))/3) (binwidth.default <- cutoff.default/15) plot(v.o <- variogram(ANN_GDD50 ~ 1, locations=ne.m), plot.numbers=TRUE) ## ----display-ordinary-v-long, out.width='0.6\\textwidth', w=8, h=6---- plot(v.o <- variogram(ANN_GDD50 ~ 1, locations=ne.m, cutoff=500000, width=20000), plot.numbers=TRUE) ## ----show-ordinary-v-short, out.width='0.6\\textwidth', w=8, h=6---- v.o <- variogram(ANN_GDD50 ~ 1, locations=ne.m, cutoff=220000) plot(v.o, plot.numbers=TRUE) ## ----model-ordinary-2, out.width='0.6\\textwidth', w=8, h=6--------- (vmf.o <- fit.variogram(v.o, vgm(psill=210000, model="Sph", range=220000, nugget=40000))) plot(v.o, plot.numbers=TRUE, model=vmf.o) ## ----krige-ok------------------------------------------------------- # dem.ne.m.sp was set up in \S6.2 "Adjusting the grid for prediction" k.ok <- krige(ANN_GDD50 ~ 1, locations=ne.m, newdata=dem.ne.m.sf, model=vmf.o, nmax=24, block=c(1000,1000)) summary(k.ok) ## ----close-sta-eith------------------------------------------------- e.ith.pt <- subset(ne.m, (ne.m$STATION_NA=="ITHACA CORNELL UNIV")) dist.pt <- st_distance(ne.m, e.ith.pt) (ix <- order(dist.pt)[1:24]) print(cbind(ne.df[ix,c('STN_NAME','STATE','ELEVATION_','ANN_GDD50')], dist=round(dist.pt[ix,]/1000,1))) ## ----krige-plot-ok-pred, out.width='0.65\\textwidth'---------------- ggplot(k.ok) + geom_sf(aes(col=var1.pred)) ## ----krige-plot-ok-pred-df, out.width='0.8\\textwidth'-------------- dim(dem.ne.m.df) dem.ne.m.df$pred.ok <- k.ok$var1.pred display.prediction.map("pred.ok", "Annual GDD base 50F, OK prediction", "GDD50") ## ----krige-plot-ok-pred-sd-df, out.width='0.8\\textwidth'----------- dem.ne.m.df$sd.ok <- sqrt(k.ok$var1.var) ggplot() + geom_point(aes(x=E, y=N, colour=sd.ok), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Annual GDD base 50F, Standard error of OK prediction") + scale_colour_distiller(name="GDD50 s.d.", space="Lab", palette="RdYlGn", direction=-1) ## ----krige-min-sd--------------------------------------------------- summary(dem.ne.m.df$sd.ok) ix <- which.min(dem.ne.m.df$sd.ok) k.ok[ix,"var1.pred"] round(100*dem.ne.m.df[ix,"sd.ok"]/k.ok$var1.pred[ix],1) ## ----krige-loocv-ok, results='hide'--------------------------------- kcv.ok <- krige.cv(ANN_GDD50 ~ 1, locations=ne.m, model=vmf.o) summary(kcv.ok$residual) ## ------------------------------------------------------------------- (loocv.ok.rmse <- sqrt(sum(kcv.ok$residual^2)/length(kcv.ok$residual))) ## ----ok-loocv-bubble-ok--------------------------------------------- bubble.sf("kcv.ok", "residual", "GDD50", "LOOCV OK residuals") ## ------------------------------------------------------------------- ne.m[which.min(kcv.ok$residual),2:6] ne.m[which.max(kcv.ok$residual),2:6] ## ----idw------------------------------------------------------------ k.idw <- idw(ANN_GDD50 ~ 1, locations=ne.m, newdata=dem.ne.m.sf, idp=2, nmax=24) summary(k.idw$var1.pred) summary(k.ok$var1.pred) ## ----idw-plot, out.width='.8\\textwidth'---------------------------- dem.ne.m.df$pred.idw <- k.idw$var1.pred display.prediction.map("pred.idw", "Annual GDD base 50F, IDW2 prediction", "GDD50") ## ----idw-loocv-1, results='hide'------------------------------------ kcv.idw <- krige.cv(ANN_GDD50 ~ 1, locations=ne.m, model=NULL) ## ----idw-loocv-2---------------------------------------------------- summary(kcv.idw$residual) summary(kcv.ok$residual) (loocv.idw.rmse <- sqrt(sum(kcv.idw$residual^2)/length(kcv.idw$residual))) (loocv.ok.rmse) ## ----idw-loocv-bubble----------------------------------------------- bubble.sf("kcv.idw", "residual", "GDD50", "IDW^2 LOOCV residuals") ## ----child-data, child='MappingRegionalClimate_Thiessen.Rnw'-------- ## ----Thiessen-set-parent, echo=FALSE, cache=FALSE------------------- set_parent('MappingRegionalClimate.Rnw') ## ----voronoi, out.width='0.9\\linewidth'---------------------------- v <- terra::voronoi(vect(ne.m), bnd=vect(state.ne.m)) class(v) plot(v) ## ----voronoi-mask, out.width='0.9\\linewidth'----------------------- v <- crop(v, vect(state.ne.m)) plot(v) ## ----predict-voroni-over, out.width='\\linewidth'------------------- class(v); dim(v) # convert from SpatVector to sf v.sf <- st_as_sf(v) summary(v.sf$geometry) plot(v.sf["STATE"], main="Thiessen polygons assigned to each State") # spatial join v3 <- st_join(v.sf, ne.m) class(v3) ## ----plot-pred-voronoi, out.width='\\linewidth'--------------------- ggplot(v3) + geom_sf(aes(bg = ANN_GDD50.y)) + labs(bg = "Annual GDD50", title = "Thiessen polygon prediction") ## ----get-value-thiessen-pt------------------------------------------ # an arbitrary point of interest require(tmaptools) (query.pt <- geocode_OSM("1256 Poplar Ridge Rd, Aurora, NY 13026")) require(sf) (pt <- st_sfc(st_point(query.pt$coords))) class(pt) # pt <- st_sfc(st_point(x=c(-76.65689, 42.73498))) st_crs(pt) <- 4326 pt <- st_transform(pt, st_crs(states.grid)) st_coordinates(pt) pt <- st_as_sf(pt) class(pt) ## ----get-value-thiessen--------------------------------------------- class(pt) class(v3) pt.pred <- st_join(pt, v3) print(pt.pred["ANN_GDD50.y"]) ## ----inter-station-dist--------------------------------------------- dm <- st_distance(ne.m) str(dm) ## ----find-nn-------------------------------------------------------- diag(dm) <- max(dm)*1.1 nn <- apply(dm, 1, which.min) str(nn) ## ----predict-nn----------------------------------------------------- nn.gdd <- ne.m[nn,"ANN_GDD50"] str(nn.gdd) ## ----compare-nn-gdd------------------------------------------------- obs <- ne.m[,"ANN_GDD50"] summary(diff <- st_drop_geometry(obs - nn.gdd)) hist(diff$ANN_GDD50, xlab = "Annual GDD50", main="Cross-validation errors, Thiessen polygons") rug(diff$ANN_GDD50) (me <- mean(diff$ANN_GDD50)) (rmse <- sqrt(sum(diff$ANN_GDD50^2)/length(diff$ANN_GDD50))) ## ------------------------------------------------------------------- rf.oob.me; rf.oob.rmse ## ----plot-loocv-voronoi, out.width='\\linewidth'-------------------- v3$resid.thiessen <- diff$ANN_GDD50 ggplot(v3) + geom_sf(aes(bg = resid.thiessen)) + scale_fill_gradient2() + labs(bg = "Annual GDD50", title = "X-validation error, Thiessen polygon prediction") ## ----child-data, child='MappingRegionalClimate_compare.Rnw'--------- ## ----compare-set-parent, echo=FALSE, cache=FALSE-------------------- set_parent('MappingRegionalClimate.Rnw') ## ----list-files-shp-no, eval=F-------------------------------------- ## list.files(".", pattern=".shp$") ## ----list-files-shp-yes, echo=F, eval=T----------------------------- ## this is the file specification on the author's system; he organizes datasets into a subdirectory named `ds', and then ## each data source into its own sub-subdirectory list.files("./ds/NEweather", pattern=".shp$") ## ----list-files-xml-no, eval=F-------------------------------------- ## list.files(".", pattern=".xml$") ## ----list-files-xml-yes, echo=F, eval=T----------------------------- list.files("./ds/NEweather", pattern=".xml$") ## ----read-july-no, eval=FALSE--------------------------------------- ## varname <- "maat7100" ## tmp <- st_read(dsn=".", layer=varname, ## int64_as_string = FALSE, quiet = TRUE) ## head(tmp) ## ----read-july-usa-yes, echo=FALSE---------------------------------- varname <- "maat7100" tmp <- st_read(dsn="./ds/NEweather", layer=varname, int64_as_string = FALSE, quiet = TRUE) head(tmp) ## ----july-4states--------------------------------------------------- ix <- (tmp$STATE %in% c("NY","NJ","PA","VT")) tmp <- tmp[ix,] names(tmp) ## ----july-crs------------------------------------------------------- st_crs(tmp)$proj4string tmp.m <- st_transform(tmp, ne.crs) st_crs(tmp.m)$proj4string ## ----add-july------------------------------------------------------- names(tmp.m) ne.m$ANN_TNORM_ <- tmp.m$ANN_TNORM_ rm(tmp, tmp.m) ## ------------------------------------------------------------------- ne.coords <- st_coordinates(ne.m) ne.m <- cbind(ne.m, E = ne.coords[,1], N = ne.coords[ , 2]) rm(ne.coords) ## ----standardized--------------------------------------------------- summary(ne.m$ANN_GDD50) summary(ne.m$ANN_TNORM_) gdd50.std <- (ne.m$ANN_GDD50 - mean(ne.m$ANN_GDD50))/sd(ne.m$ANN_GDD50) t.ann.std <- (ne.m$ANN_TNORM_ - mean(ne.m$ANN_TNORM_))/sd(ne.m$ANN_TNORM_) summary(gdd50.std) summary(t.ann.std) ## ----check-corr-t-gdd----------------------------------------------- cor(ne.m$ANN_GDD50, ne.m$ANN_TNORM_) cor(gdd50.std, t.ann.std) plot(gdd50.std, t.ann.std, asp=1, col = ne.m$STATE, pch = 20, xlab="Standardized annual GDD50", ylab="Standardized mean annual T") abline(0,1); grid() ## ----build-two-models----------------------------------------------- summary(m.ols.t <- lm(t.ann.std~ sqrt(ELEVATION_)+N+E, data=ne.m)) summary(m.ols.gdd <- lm(gdd50.std ~ sqrt(ELEVATION_)+N+E, data=ne.m)) coef(m.ols.t) coef(m.ols.gdd) coef(m.ols.t)/coef(m.ols.gdd) ## ----show-ols-fits-two, out.width='\\textwidth', w=12, h=6---------- summary(ne.m$ols.resid.t <- residuals(m.ols.t)) summary(ne.m$ols.resid.gdd <- residuals(m.ols.gdd)) p1 <- bubble.sf("ne.m", "ols.resid.t", "OLS residuals", "residuals, standardized annual T") p2 <- bubble.sf("ne.m", "ols.resid.gdd", "xx", "residuals, standardized annual GDD50") require(gridExtra) grid.arrange(p1, p2, ncol = 2) ## ----compare-resids-two--------------------------------------------- summary(ne.m$ols.resid.diff <- ne.m$ols.resid.t - ne.m$ols.resid.gdd) bubble.sf("ne.m", "ols.resid.diff", "difference", "Difference of residuals, standardized T - standardized GDD") ## ----model-resid-vgms----------------------------------------------- require(gstat) v.r.ols.t <- variogram(ols.resid.t ~ 1, locations=ne.m, cutoff=120000, width=12000) (vmf.r.ols.t <- fit.variogram(v.r.ols.t, vgm(psill=0.1, model="Exp", range=20000, nugget=0.02))) # v.r.ols.gdd <- variogram(ols.resid.gdd ~ 1, locations=ne.m, cutoff=120000, width=12000) (vmf.r.ols.gdd <- fit.variogram(v.r.ols.gdd, vgm(psill=0.12, model="Exp", range=20000, nugget=0.02))) # # a common y-axis scale ymax <- max(sum(vmf.r.ols.t[,"psill"]), sum(vmf.r.ols.gdd[,"psill"]))*1.1 p1 <- plot(v.r.ols.t, pl=T, model=vmf.r.ols.t, ylim = c(0, ymax)) p2 <- plot(v.r.ols.gdd, pl=T, model=vmf.r.ols.gdd, ylim = c(0, ymax)) print(p1, split=c(1,1,1,2), more=T) print(p2, split=c(1,2,1,2), more=F) ## ----fit-gls-two---------------------------------------------------- # require(nlme) require(nlme) (p.nugget <- vmf.r.ols.t[1,"psill"]/sum(vmf.r.ols.t[,"psill"]) + 0.001) m.gls.t <- gls(model=t.ann.std ~ sqrt(ELEVATION_) + N + E, data=ne.m, correlation=corExp( value=c(vmf.r.ols.t[2,"range"], p.nugget), form=~E + N, nugget=TRUE)) # (p.nugget <- vmf.r.ols.gdd[1,"psill"]/sum(vmf.r.ols.gdd[,"psill"]) + 0.001) m.gls.gdd <- gls(model=gdd50.std ~ sqrt(ELEVATION_) + N + E, data=ne.m, correlation=corExp( value=c(vmf.r.ols.gdd[2,"range"], p.nugget), form=~E + N, nugget=TRUE)) summary(m.gls.t) summary(m.gls.gdd) ## ----compare-coeff-gls-two------------------------------------------ coefficients(m.gls.t) coefficients(m.gls.gdd) coefficients(m.gls.t)/coefficients(m.gls.gdd) ## ----compare-corr-struct-gls-two------------------------------------ intervals(m.gls.t)$corStruct intervals(m.gls.gdd)$corStruct intervals(m.gls.t)$corStruct/intervals(m.gls.gdd)$corStruct ## ----gls-resid-two-------------------------------------------------- ne.m$gls.resid.t <- residuals(m.gls.t) ne.m$gls.resid.gdd <- residuals(m.gls.gdd) ## ----show-vgm-gls-corr-std------------------------------------------ (p.nugget <- intervals(m.gls.t)$corStruct["nugget","est."]) (t.sill <- var(ne.m$gls.resid.t)) (vmf.r.gls.t <- vgm(psill=t.sill*(1-p.nugget), model="Exp", range=intervals(m.gls.t)$corStruct["range","est."], nugget=t.sill*p.nugget)) v.r.gls.t <- variogram(gls.resid.t ~ 1, locations=ne.m, cutoff=120000, width=12000) # (p.nugget <- intervals(m.gls.gdd)$corStruct["nugget","est."]) (t.sill <- var(ne.m$gls.resid.gdd)) (vmf.r.gls.gdd <- vgm(psill=t.sill*(1-p.nugget), model="Exp", range=intervals(m.gls.gdd)$corStruct["range","est."], nugget=t.sill*p.nugget)) v.r.gls.gdd <- variogram(gls.resid.gdd ~ 1, locations=ne.m, cutoff=120000, width=12000) ymax <- max(sum(vmf.r.gls.t[,"psill"]), sum(vmf.r.gls.gdd[,"psill"]))*1.1 p1 <- plot(v.r.gls.t, model=vmf.r.gls.t, pl=T, ylim = c(0, ymax)) p2 <- plot(v.r.gls.gdd, model=vmf.r.gls.gdd, pl=T, ylim = c(0, ymax)) print(p1, split=c(1,1,1,2), more=T) print(p2, split=c(1,2,1,2), more=F) ## ----gls-predict-two------------------------------------------------ dem.ne.m.df$pred.gls.t <- predict(m.gls.t, newdata=dem.ne.m.df) dem.ne.m.df$pred.gls.gdd <- predict(m.gls.gdd, newdata=dem.ne.m.df) dem.ne.m.df$diff.gls.t.gdd <- dem.ne.m.df$pred.gls.t - dem.ne.m.df$pred.gls.gdd ## ----plot-4km-gdd-t, out.width='.8\\textwidth'---------------------- (std.pred.lim <- c(min(dem.ne.m.df[,c("pred.gls.t","pred.gls.gdd")]), max(dem.ne.m.df[,c("pred.gls.t","pred.gls.gdd")]))) display.prediction.map("pred.gls.t", "Mean Annual Temperature, standardized, GLS prediction", "GDD50", std.pred.lim, .palette="RdGy") display.prediction.map("pred.gls.gdd", "Annual GDD, base 50F, standardized, GLS prediction", "GDD50", std.pred.lim, .palette="RdGy") ## ----diffs-4km-gls-t-gdd-------------------------------------------- summary(dem.ne.m.df$diff.gls.t.gdd <- dem.ne.m.df$pred.gls.t - dem.ne.m.df$pred.gls.gdd) ## ----plot-4km-gdd-t-diff, out.width='.8\\textwidth'----------------- display.difference.map("diff.gls.t.gdd", "Difference, MAT - GDD50, standardized", "+/- s.d.", .palette="BrBG") ## ----ok-gls-resid-t-gdd-std----------------------------------------- ok.gls.resid.t <- krige(gls.resid.t ~ 1, loc=ne.m, newdata=dem.ne.m.sf, model=vmf.r.gls.t) ok.gls.resid.gdd <- krige(gls.resid.gdd ~ 1, loc=ne.m, newdata=dem.ne.m.sf, model=vmf.r.gls.gdd) dem.ne.m.df$ok.gls.resid.t <- ok.gls.resid.t$var1.pred dem.ne.m.df$ok.gls.resid.gdd <- ok.gls.resid.gdd$var1.pred ## ----ok-gls-resid-t-gdd-std-maps, out.width='.8\\textwidth'--------- ggplot() + geom_point(aes(x=E, y=N, colour=ok.gls.resid.t), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Residuals from GLS trend surface, MAT (std)") + scale_colour_distiller(name="T (std)", space="Lab", palette="RdBu") ggplot() + geom_point(aes(x=E, y=N, colour=ok.gls.resid.gdd), data=dem.ne.m.df) + xlab("E") + ylab("N") + coord_fixed() + ggtitle("Residuals from GLS trend surface, GDD base 50F (std)") + scale_colour_distiller(name="GDD50 (std)", space="Lab", palette="RdBu") ## ----ok-gls-resid-t-gdd-diff, out.width='.8\\textwidth'------------- summary(dem.ne.m.df$ok.gls.resid.diff.t.gdd <- dem.ne.m.df$ok.gls.resid.t - dem.ne.m.df$ok.gls.resid.gdd) display.difference.map("ok.gls.resid.diff.t.gdd", "Difference: kriged residuals from GLS trend surface", "+/- s.d.", .palette="BrBG") ## ----final-gls-rk-std, out.width='.8\\textwidth'-------------------- summary(dem.ne.m.df$pred.rkgls.std.t <- dem.ne.m.df$pred.gls.t + dem.ne.m.df$ok.gls.resid.t) summary(dem.ne.m.df$pred.rkgls.std.gdd <- dem.ne.m.df$pred.gls.gdd + dem.ne.m.df$ok.gls.resid.gdd) display.prediction.map("pred.rkgls.std.t", "GLS-RK prediction, Mean Annual Temperature, standardized", "MAAT (std)", std.pred.lim, .palette="YlOrBr") display.prediction.map("pred.rkgls.std.gdd", "GLS-RK prediction, Annual GDD, base 50F, standardized", "GDD50 (std)", std.pred.lim, .palette="YlOrBr") ## ----final-gls-rk-std-diff, out.width='.8\\textwidth'--------------- # summary(dem.ne.m.df$diff.rkgls.std.t.gdd <- (dem.ne.m.df$pred.rkgls.std.t - dem.ne.m.df$pred.rkgls.std.gdd)) display.difference.map("diff.rkgls.std.t.gdd", "GLS-RK predictions, difference, MAT - GDD50, standardized", "+/- s.d.", .palette="BrBG") ## ----build-rf-two--------------------------------------------------- require(randomForest) m.rf.std.t <- randomForest(t.ann.std ~ ELEVATION_ + N + E + dist.lakes + dist.coast, data=ne.df, ntree=1200, importance=TRUE) m.rf.std.gdd <- randomForest(gdd50.std ~ ELEVATION_ + N + E, data=ne.df, ntree=1200, importance=TRUE) ## ----compare-rf-importance------------------------------------------ randomForest::importance(m.rf.std.t) randomForest::importance(m.rf.std.gdd) ## ----compare-rf-11-oob, out.width='0.45\\linewidth'----------------- plot(t.ann.std ~ predict(m.rf.std.t, newdata=ne.m), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest", ylab="Actual", main="Mean Annual T (std)") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) plot(gdd50.std ~ predict(m.rf.std.t, newdata=ne.m), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest", ylab="Actual", main="Annual GDD50 (std)") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) plot(t.ann.std ~ predict(m.rf.std.t), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest (OOB)", ylab="Actual", main="Mean Annual T (std)") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) plot(gdd50.std ~ predict(m.rf.std.t), col=ne.m$STATE, pch=20, asp=1, xlab="Fitted by random forest (OOB)", ylab="Actual", main="Annual GDD50 (std)") legend("topleft", levels(ne.m$STATE), pch=20, col=1:4) grid(); abline(0,1) ## ----predict-rf-two------------------------------------------------- dem.ne.m.df$pred.rf.std.t <- predict(m.rf.std.t, newdata=dem.ne.m.df) dem.ne.m.df$pred.rf.std.gdd <- predict(m.rf.std.gdd, newdata=dem.ne.m.df) (std.pred.lim <- c(min(dem.ne.m.df[,c("pred.rf.std.t","pred.rf.std.gdd")]), max(dem.ne.m.df[,c("pred.rf.std.t","pred.rf.std.gdd")]))) display.prediction.map("pred.rf.std.t", "Mean Annual Temperature, standardized, RF prediction", "degrees C", std.pred.lim, .palette="YlOrBr") display.prediction.map("pred.rf.std.gdd", "Annual GDD, base 50F, standardized, RF prediction", "GDD50", std.pred.lim, .palette="YlOrBr") ## ----predict-rf-two-diff-------------------------------------------- summary(dem.ne.m.df$diff.rf.std.t.gdd <- dem.ne.m.df$pred.rf.std.t - dem.ne.m.df$pred.rf.std.gdd) display.difference.map("diff.rf.std.t.gdd", "RF predictions, difference, MAT - GDD50, standardized", "+/- s.d.", .palette="BrBG") ## ----child-data, child='MappingRegionalClimate_answers.Rnw'--------- ## ----answers-set-parent, echo=FALSE, cache=FALSE-------------------- set_parent('MappingRegionalClimate.Rnw') ## ----get-rt-row, echo=FALSE----------------------------------------- ix.rt.l <- which(row.names(m.rt$frame)=="2") ix.rt.r <- which(row.names(m.rt$frame)=="3") ## ----child-data, child='MappingRegionalClimate_challenge.Rnw'------- ## ----challenge-set-parent, echo=FALSE, cache=FALSE------------------ set_parent('MappingRegionalClimate.Rnw') ## ----child-data, child='MappingRegionalClimate_colours.Rnw'--------- ## ----colours-set-parent, echo=FALSE, cache=FALSE-------------------- set_parent('MappingRegionalClimate.Rnw') ## ----seq-brewers, out.width='\\linewidth', h=6, w=8----------------- library(RColorBrewer) display.brewer.all(type="seq") ## ----display-annual-colored, out.width='\\textwidth'---------------- ggplot(data=ne.df) + aes(x=E, y=N) + geom_point(aes(size=ANN_GDD50, colour=ANN_GDD50), shape=20) + scale_colour_distiller(space="Lab", palette="Greens") + xlab("E") + ylab("N") + coord_fixed() ## ----seq-brewers-div, out.width='\\linewidth', h=6, w=8------------- display.brewer.all(type="div")