library("remote") ###### EXAMPLE I ########################################################## ########################################################################### ## library("rworldmap") library("rgdal") library("rgeos") library("gridExtra") library("RColorBrewer") ## load data data("vdendool", package = "remote") data("coastsCoarse", package = "rworldmap") ## calculate 4 leading modes nh_modes <- eot(x = vdendool, y = NULL, n = 4, reduce.both = FALSE, standardised = FALSE, verbose = TRUE) ## create coastal outlines ster <- CRS("+proj=stere +lat_0=90 +lon_0=-45") xmin <- -180 xmax <- 180 ymin <- 20 ymax <- 90 # Coordinates for bounding box bb <- cbind(x = c(xmin, xmin, xmax, xmax, xmin), y = c(ymin, ymax, ymax, ymin, ymin)) # Create bounding box SP <- SpatialPolygons(list(Polygons(list(Polygon(bb)), "1")), proj4string = CRS(proj4string(coastsCoarse))) gI <- gIntersects(coastsCoarse, SP, byid = TRUE) out <- vector(mode = "list", length = length(which(gI))) ii <- 1 for (i in seq(along = gI)) if (gI[i]) { out[[ii]] <- gIntersection(coastsCoarse[i, ], SP) row.names(out[[ii]]) <- row.names(coastsCoarse)[i] ii <- ii + 1 } nh_coasts <- do.call("rbind", out) nh_coasts_ster <- spTransform(nh_coasts, ster) lout <- list("sp.lines", nh_coasts_ster, col = "grey30", grid = TRUE) ## define colours clrs <- colorRampPalette(rev(brewer.pal(9, "RdBu"))) ## project modes to polar stereographic CRS mode <- lapply(seq(nmodes(nh_modes)), function(i) { projectRaster(nh_modes[[i]]@r_predictor, crs = ster) }) ## create title for each panel title <- lapply(seq(mode), function(i) { paste("Mode ", i, " : ", "EV = ", round(if (i > 1) { nh_modes[[i]]@cum_exp_var * 100 - nh_modes[[i - 1]]@cum_exp_var * 100 } else { nh_modes[[i]]@cum_exp_var * 100 }, 1), " : ", "BP = ", round(nh_modes[[i]]@coords_bp[, 1], 1), ", ", round(nh_modes[[i]]@coords_bp[, 2], 1), sep = "") }) ## create panel plots p <- lapply(seq(mode), function(i) { spplot(mode[[i]], sp.layout = lout, main = list(title[[i]], cex = 0.7), col.regions = clrs(1000), at = seq(-1, 1, 0.2), par.settings = list(axis.line = list(col = 0)), colorkey = list(height = 0.75, width = 1)) }) f <- function(...) grid.arrange(..., heights = 1, ncol = 2) do.call(f, p) #Final plot vp1 <- viewport(x = 0.5, y = 0.5, height = 1, width = 1, just = c("centre", "centre"), name = "text.vp") pushViewport(vp1) grid.text("DJF HGT 700 hPa NCEP/NCAR Reanalysis (1948 - 1998)", y = unit(0.52, "npc")) popViewport() ########################################################################### ########################################################################### ###### EXAMPLE II ######################################################### ########################################################################### data("australiaGPCP", package = "remote") data("pacificSST", package = "remote") ### deseason data sst_pred <- deseason(pacificSST, cycle.window = 12) gpcp_resp <- deseason(australiaGPCP, cycle.window = 12) ### denoise data (keeping 90 % of the variance) sst_pred_dns <- denoise(sst_pred, expl.var = 0.9) gpcp_resp_dns <- denoise(gpcp_resp, expl.var = 0.9) ### calculate first 3 leading modes aus_modes <- eot(x = sst_pred_dns, y = gpcp_resp_dns, n = 1, standardised = FALSE, reduce.both = FALSE, verbose = TRUE) ### plot first EOT showing the location of the base point plot(aus_modes, y = 1, show.bp = TRUE, arrange = "long") ########################################################################### ########################################################################### ###### EXAMPLE III ######################################################## ########################################################################### library("ggplot2") library("reshape") library("latticeExtra") library("gridExtra") load("gimmsKiliNDVI.RData") load("modisKiliNDVI.RData") ### define index for training data pred_ind <- 1:24 ### create training (pred) and evaluation (eval) sets mod_stck_pred <- modisKiliNDVI[[pred_ind]] mod_stck_eval <- modisKiliNDVI[[-pred_ind]] gimms_stck_pred <- gimmsKiliNDVI[[pred_ind]] gimms_stck_eval <- gimmsKiliNDVI[[-pred_ind]] ### calculate EOT ndvi_modes <- eot(x = gimms_stck_pred, y = mod_stck_pred, n = 10, standardised = FALSE, reduce.both = FALSE, verbose = TRUE) ### calculate number of modes necessary for explaining 98% variance nm <- nXplain(ndvi_modes, 0.98) ### prediction using claculated intercept, slope and GIMMS NDVI values mod_predicted <- predict(object = ndvi_modes, newdata = gimms_stck_eval, n = nm) mod_observed <- mod_stck_eval pred_vals <- getValues(mod_predicted) obs_vals <- getValues(mod_observed) ### error scores ME <- colMeans(pred_vals - obs_vals, na.rm = TRUE) MAE <- colMeans(abs(pred_vals - obs_vals), na.rm = TRUE) RMSE <- sqrt(colMeans((pred_vals - obs_vals)^2, na.rm = TRUE)) R <- diag(cor(pred_vals, obs_vals, use = "complete.obs")) Rsq <- R * R ### visualise error scores scores <- data.frame(ME, MAE, RMSE, R, Rsq) melt_scores <- melt(scores) p <- ggplot(melt_scores, aes(factor(variable), value)) p <- p + geom_boxplot() + theme_bw() + xlab("") + ylab("") print(p) ## spatial residuals plot clrs_hcl <- function(n) { hcl(h = seq(230, 0, length.out = n), c = 60, l = 70, fixup = TRUE) } clrs_ndvi <- colorRampPalette(brewer.pal(11, "BrBG")) clrs_resids <- colorRampPalette(c(clrs_hcl(2)[2], "grey10", clrs_hcl(2)[1])) plotResid <- function(pred.rst, resp.rst, bg = NULL, txt = "", n = 5000) { if (is.null(bg)) bg <- resp.rst lm1 <- resp.rst - pred.rst highest.pix <- pred.rst highest.pix[] <- NA highest.pix[which(values(lm1) > 0.1)] <- lm1[which(values(lm1) > 0.1)] lowest.pix <- pred.rst lowest.pix[] <- NA lowest.pix[which(values(lm1) < -0.1)] <- lm1[which(values(lm1) < -0.1)] lev <- spplot(bg, col.regions = clrs_ndvi(1200), at = seq(-0.1, 1, 0.01), as.table = TRUE, main = "Predicted NDVI", colorkey = list(space = "top", width = 1, height = 0.75)) lowest.lev <- spplot(lowest.pix, col.regions = clrs_resids(1200), at = seq(-1, 1, 0.01), colorkey = FALSE, as.table = TRUE) highest.lev <- spplot(highest.pix, col.regions = clrs_resids(1200), at = seq(-1, 1, 0.01), colorkey = FALSE, as.table = TRUE, panel = function(...) { panel.levelplot(...) panel.text(labels = txt, col = "black", x = 345000, y = 9682000, cex = 0.7)}) out <- lev + as.layer(lowest.lev) + as.layer(highest.lev) return(out) } resid_plot <- lapply(seq(ncol(pred_vals)), function(i) { panel.name <- substr(strsplit(names(mod_stck_eval)[i], "_")[[1]][5], 2, 7) plotResid(pred.rst = mod_observed[[i]], resp.rst = mod_predicted[[i]], bg = NULL, txt = panel.name, n = 5000) }) outLayout <- function(x, y) { update(c(x, y, layout = c(4, length(resid_plot)/4)), between = list(y = 0.3, x = 0.3)) } out_res <- Reduce(outLayout, resid_plot) ## png(paste("resids_", Sys.Date(), ".png", sep = ""), ## width = 19, height = 27, units = "cm", res = 300) print(out_res) downViewport(trellis.vpname(name = "figure")) vp1 <- viewport(x = 0.5, y = 0, height = 0.07, width = 0.75, just = c("centre", "top"), name = "key.vp") pushViewport(vp1) draw.colorkey(key = list(col = clrs_resids(1200), width = 1, at = seq(-1, 1, 0.01), space = "bottom"), draw = TRUE) upViewport() vp2 <- viewport(x = 0.5, y = 0, height = 0.07, width = 0.075, just = c("centre", "top"), name = "key.vp") pushViewport(vp2) draw.colorkey(key = list(col = "white", width = 1, at = seq(-0.1, 0.1, 0.1), tick.number = 2, space = "bottom"), draw = TRUE) grid.text("Residuals", x = 0.5, y = 0, just = c("centre", "top")) ## dev.off() ## time series comparison set.seed(123) excl_cells <- attributes(na.exclude(mod_predicted[[1]][]))$na.action "%w/o%" <- function(x, y) x[!x %in% y] valid_cells <- seq(mod_predicted[[1]][]) %w/o% excl_cells kili_sites <- sample(valid_cells, 50, replace = FALSE) time <- seq.Date(as.Date("2005-01-01"), as.Date("2006-12-31"), by = "months") p_vals <- as.data.frame(extract(mod_predicted, kili_sites)) names(p_vals) <- time p_vals$pID <- as.character(kili_sites) p_vals_long <- melt(p_vals) o_vals <- as.data.frame(extract(mod_observed, kili_sites)) names(o_vals) <- time o_vals$pID <- as.character(kili_sites) o_vals_long <- melt(o_vals) pred_p <- xyplot(value ~ as.Date(variable) | pID, data = p_vals_long, type = "l", as.table = TRUE, layout = c(5, 10), par.settings = list(layout.heights = list(strip = 0.75)), between = list(y = 0.2, x = 0.2), col = "black", ylab = "NDVI", xlab = "", panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(h = mean(y), lty = 2) }, strip = strip.custom( bg = "grey20", par.strip.text = list(col = "white", font = 2, cex = 0.8)), ylim = c(0, 1), xlim = c(as.Date(p_vals_long$variable[1]), as.Date(p_vals_long$variable[nrow(p_vals_long)]))) obs_p <- xyplot(value ~ as.Date(variable) | pID, data = o_vals_long, type = "l", as.table = TRUE, col = "grey60", panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(h = mean(y), lty = 2, col = "grey60") }) ## png(paste("dnsc_vs_mod_", Sys.Date(), ".png", sep = ""), ## width = 19, height = 27, units = "cm", res = 300) pred_p + as.layer(obs_p) ## dev.off() ########################################################################### ###########################################################################