2_GTD_GCM_Comparison_v1.R

################### COMPARISON OF SOIL MOISTURE ###################
# 
# Daniel Schlaepfer, 2015-2016
# 
# Comparison of mean monthly soil moisture patterns between SOILWAT and GCMs. The variable 
# 'mrso' (total soil moisture content) was extracted for seven GCMs under historical and 
# future (RCP4.5 and RCP8.5) scenarios from non-downscaled data from the ESGF node 
# https://pcmdi.llnl.gov/. We calculated mean monthly values for the periods of 1980-2005 
# and 2070-2099 for each of our simulated raster cells and compared agreement with 
# equivalent soil moisture values from SOILWAT output. We estimated agreement between 
# models with Duveiller's lambda, which was the best performing symmetric agreement index 
# (Duveiller et al. 2016). lambda ranges between 0 and 1 where 0 indicates no agreement and 1 
# is perfect agreement. lambda is proportional to Pearson's correlation index and accounts for 
# systematic and unsystematic bias.
#
# References:
# Duveiller G, Fasbender D, Meroni M (2016) Revisiting the concept of a symmetric index of 
# agreement for continuous datasets. Scientific Reports, 6, 19401.
#
###################################################################


#---Actions
version <- "v1"
do_extract_from_GCMs <- FALSE
do_compare <- FALSE
do_tables <- TRUE
do_maps <- TRUE

#---R packages
pkg_reqd <- c("raster", "ncdf4", "rgdal", "maps", "reshape2")
has_notloaded <- sapply(pkg_reqd,
	function(lib) !require(lib, character.only = TRUE, quietly = FALSE))

if (any(has_notloaded)) {
	sapply(pkg_reqd[has_notloaded],
		function(lib) install.packages(lib))
	has_loaded <- sapply(pkg_reqd[has_notloaded],
		function(lib) require(lib, character.only = TRUE, quietly = FALSE))
	stopifnot(has_loaded)
}


#---Load functions
dir.gtd <- "/PATH_TO_PROJECT/Product_PowellCenter/6_Projects_Year1"
source(file.path(dir.gtd, "Prj03_GlobalVulnerability", "4_Analysis", "4_Analysis_v4", "5a1_GTD_Prj03v4_Helper.R"), verbose = FALSE, chdir = FALSE)


#---Directories
dir.gtd <- "/PATH_TO_PROJECT/Product_PowellCenter"
dir.prj <- file.path(dir.gtd, "1_Comparison_with_GCM_SoilMoisture")
dir.ext <- file.path(dir.prj, "4_ExtractedData")
dir.out <- file.path(dir.prj, "6_Results")

dir.sana <- file.path(dir.gtd, "6_Projects_Year1", "GTD_SharedAnalysis")
dir.gis <- file.path(dir.sana, "1_GISdata")

dir.gcm <- "/PATH_TO_DATA/Weather_Future/ESGF/Downloads"


#---Settings
WGS84 <- raster::crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

regions_N <- 1:6
label.regions <- c("South America", "Southern Africa", "Eastern Asia", "Western & Central Asia", "Western Mediterranean", "North America")

startyr <- 1979 + 1
endyr <- 2010
nyears <- endyr - startyr + 1
deltaFutureToSimStart_yr <- 90
years <- list(historical = temp <- startyr:endyr, future = temp + deltaFutureToSimStart_yr)
mmonths <- 1:12
gcm_soil_N <- 6

currentSc <- "Current"

reqRCPs <- c("RCP45", "RCP85")
scenarios <- c("historical", reqRCPs)

reqGCMs <- c("CanESM2", "CESM1-CAM5", "CSIRO-Mk3-6-0", "EC-EARTH", "FGOALS-g2", "FGOALS-s2", "GFDL-CM3", "GISS-E2-R", "HadGEM2-CC", "HadGEM2-ES", "inmcm4", "IPSL-CM5A-MR", "MIROC-ESM", "MIROC5", "MPI-ESM-MR", "MRI-CGCM3")

#vars <- c("tran", "evspsblsoi", "evspsblveg", "tsl", "mrso", "mrsos", "mrlsl", "mrros", "mrro", "lai", "landCoverFrac", "residualFrac", "baresoilFrac", "shrubFrac", "treeFrac", "c3PftFrac", "c4PftFrac", "grassFrac", "cropFrac")
vars <- c("mrso", "mrsos", "mrlsl")

#http://www-pcmdi.llnl.gov/ipcc/standard_output.html
#Variable		Variable Long Name				Units
#tran			Transpiration					kg m-2 s-1
#evspsblsoi		Water Evaporation from Soil		kg m-2 s-1
#evspsblveg		Evaporation from Canopy			kg m-2 s-1
#
#tsl				Temperature of soil				K
#
#landCoverFrac	Plant Functional Type Grid Fraction	%
#residualFrac	Fraction of Grid Cell that is Land but Neither Vegetation-Covered nor Bare Soil											%
#baresoilFrac	Bare Soil Fraction				%
#shrubFrac		Shrub Fraction					%
#treeFrac		Tree Cover Fraction				%
#c3PftFrac		Total C3 PFT Cover Fraction		%
#c4PftFrac		Total C4 PFT Cover Fraction		%
#grassFrac		Natural Grass Fraction			%
#cropFrac		Crop Fraction					%
#
#mrso			Total Soil Moisture Content		kg m-2
	# mrso = water in all phases summed over all soil layers, and averaged over the land portion of the grid cell (i.e., compute by dividing the total mass of water contained in the soil layer of the grid cell by the land area in the grid cell); report as "missing" or 0.0 where the land fraction is 0.
#mrsos			Moisture in Upper Portion of Soil Column		kg m-2
	# mrsos = water in all phases in the upper 0.1 meters of soil, and averaged over the land portion of the grid cell (i.e., compute by dividing the total mass of water contained in the soil layer of the grid cell by the land area in the grid cell); report as "missing" or 0.0 where the land fraction is 0; the CMOR singleton dimension default value of 0.1 m can be overridden, if absolutely necessary, by redefining axis "depth1".
#mrlsl			Water Content per Soil Layer	kg m-2
#mrros			Surface Runoff					kg m-2 s-1
#mrro			Total Runoff					kg m-2 s-1
#lai			Leaf Area Index					1

#sftlf			land_area_fraction				%
#rootd			root_depth						m
	# rootd = As a function of longitude and latitude, report the maximum soil depth reachable by plant roots, i.e., the maximum soil depth from which they can extract moisture (if defined in model); report as "missing" or 0.0 where the land fraction is 0

 

#---Our raster and cells
mask_Any33Cond <- raster::raster(file.path(dir.gis, "Any33Cond", "StudyAreaMask_Any33Cond.asc"))
raster::crs(mask_Any33Cond) <- WGS84
xy_Any33Cond <- raster::rasterToPoints(mask_Any33Cond)[, c("x", "y")]
rotate_coords <- function(xy, init = 0) {
	xy[, 1] <- xy[, 1] + ifelse(xy[, 1] < init, 360, 0)
	xy
}
	
if (do_extract_from_GCMs) {
	#---Extract GCM data
	for (iv in seq_along(vars)) {
		ftemp <- file.path(dir.ext, paste0("Extract_", vars[iv], "_ts_", version, ".rds"))
		if (file.exists(ftemp)) next

		k <- 1
		# Find available netCDF files for variable
		ftemps1 <- list.files(dir.gcm, pattern = paste0(vars[iv], "_"))
	
		availGCMs <- sort(unique(sapply(strsplit(ftemps1, split = "_", fixed = TRUE), function(x) x[3])))
		usedGCMs <- availGCMs[availGCMs %in% reqGCMs]
	
		if (length(ftemps1) > 0 && length(usedGCMs) == 0) next
	
		sdepths_m <- lapply(ftemps1, function(f) {
						temp <- ncdf4::nc_open(file.path(dir.gcm, f))
						sdepth <- ncdf4::ncvar_get(temp, "depth")
						ncdf4::nc_close(temp)
						sdepth
					})
		names(sdepths_m) <- ftemps1

		var_soil_N <- if (vars[iv] == "mrlsl") {
				temp <- max(lengths(sdepths_m), na.rm = TRUE)
				if (is.finite(temp)) temp else gcm_soil_N
			} else {
				1L
			}
		save(sdepths_m, var_soil_N, file = sub("_ts_", "_sdepth_", ftemp))
	
		#---Result container
		has_new_data <- FALSE
		res <- array(NA, dim = c(nrow(xy_Any33Cond),
								length(mmonths) * nyears,
								length(scenarios),
								length(usedGCMs),
								var_soil_N),
						 dimnames = list(NULL, NULL, scenarios, usedGCMs, NULL))

		print(object.size(res), unit = "Mb") #4.8 Gb


		for (igcm in seq_along(usedGCMs)) {
			# Select among files those for a specific GCM
			ftemps2 <- ftemps1[grepl(paste0("_", usedGCMs[igcm], "_"), ftemps1, ignore.case = TRUE)]
	
			if (length(ftemps2) == 0) {
				print(paste("No file located for", sQuote(usedGCMs[igcm])))
				next
			}
		
			for (isc in seq_along(scenarios)) {
				if (sum(res[, , isc, igcm, ], na.rm = TRUE) > 0) {
					print(paste("Output container 'res' already contains extracted data for", sQuote(usedGCMs[igcm]), sQuote(scenarios[isc])))
					next
				}
			
				# Select among files those for a specific scenario
				ftemps3 <- ftemps2[grepl(paste0("_", scenarios[isc], "_"), ftemps2, ignore.case = TRUE)]
	
				if (length(ftemps3) == 0) {
					print(paste("No file located for", sQuote(usedGCMs[igcm]), sQuote(scenarios[isc])))
					next
				}
				sc_years <- if (scenarios[isc] == "historical") years[["historical"]] else years[["future"]]
			
				for (it in seq_along(ftemps3)) {
					stopIfNotEqualSpaced <- TRUE
					fdat <- try(raster::brick(file.path(dir.gcm, ftemps3[it]), level = 1), silent = TRUE)
					#funit <- fdat@data@unit
				
					if (inherits(fdat, "try-error")) {
						stopIfNotEqualSpaced <- FALSE
						fdat <- try(raster::brick(file.path(dir.gcm, ftemps3[it]), level = 1, stopIfNotEqualSpaced = stopIfNotEqualSpaced), silent = TRUE) # stopIfNotEqualSpaced = FALSE would treat them as 'index-position' rasters and coordinates would be incorrect
						raster::crs(fdat) <- WGS84
						
						if (inherits(fdat, "try-error")) {
							print(paste("Cannot load raster", ftemps3[it], ":", fdat))
							next
						} else {
							print(paste("Raster cells are not equally spaced:", ftemps3[it], "treated as 'index-position' rasters"))
						}
					}
					
					# Extract time dimension and determine which index is the start year
					temp <- as.POSIXlt(fdat@z$Date)
					fdates <- data.frame(year = 1900 + temp$year, month = 1 + temp$mon)
					tint <- intersect(fdates$year, sc_years)
					nyears_avail <- length(tint)
				
					if (nyears_avail == 0) {
						print(paste("No relevant years in raster", ftemps3[it]))
						next
					}
											
					i_firstlayer <- which(fdates$year == tint[1] & fdates$month == 1)
					# Cell coordinates within (-180, +180) vs. (0, 360) longitude coordinates
					fext <- extent(fdat)
					xy <- if (fext@xmax > 180 && fext@xmin < 0) rotate_coords(xy_Any33Cond, init = fext@xmin) else xy_Any33Cond
			
					# Read netCDF data
					print(paste(Sys.time(), "step:", k, vars[iv], usedGCMs[igcm], scenarios[isc], "for years", paste(range(tint), collapse = "-")))
				
					for (ilevel in seq_len(var_soil_N)) {
						if (ilevel > 1)
							fdat <- try(raster::brick(file.path(dir.gcm, ftemps3[it]), level = ilevel, stopIfNotEqualSpaced = stopIfNotEqualSpaced), silent = TRUE)
						
						if (inherits(fdat, "try-error"))
							next
					
						# Extract requested data
						temp <- raster::extract(fdat, y = xy, layer = i_firstlayer, nl = length(mmonths) * nyears_avail, method = "bilinear", df = TRUE, factors = FALSE)
						if (length(itemp <- grep("ID", colnames(temp))) > 0)
							temp <- temp[, -itemp]

						# Determine for which years/months we extracted data successfully
						ctemp <- strsplit(gsub("X", "", colnames(temp)), split = ".", fixed = TRUE)
						ctemp <- sapply(ctemp, function(x) paste(x[1:2], collapse = ""))
						ctarget <- apply(cbind(rep(sc_years, each = 12), rep(mmonths, times = length(nyears))), 1, function(x) paste0(x[1], formatC(x[2], flag = "0", width = 2)))
						icols <- match(ctemp, table = ctarget, nomatch = 0)

						# Copy values to our result container
						res[, icols, isc, igcm, ilevel] <- as.matrix(temp)
						has_new_data <- TRUE
					
					} # end loop along var_soil_N

					# Save temporary copy
					if (has_new_data) {
						temp <- try(saveRDS(res, file = file.path(dir.ext, paste0("Extract_", k, "_", vars[iv], "_", version, ".rds"))), silent = TRUE)
						if (!inherits(temp, "try-error"))
							unlink(file.path(dir.ext, paste0("Extract_", k - 1, "_", vars[iv], "_", version, ".rds")))
						k <- k + 1
						gc()
					}
				} # end loop along ftemps3
			} # end loop along scenarios
		} # end loop along usedGCMs
	
		if (has_new_data)
			saveRDS(res, file = ftemp)
			
			rm(res, fdat)
	} # end loop along vars




	#---Aggregate time series to mean monthly values
	for (iv in seq_along(vars)) {
		ftemp <- file.path(dir.ext, paste0("Extract_", vars[iv], "_meanmonthly_", version, ".rds"))
		if (file.exists(ftemp)) next
	
		ftempin <- file.path(dir.ext, paste0("Extract_", vars[iv], "_ts_", version, ".rds"))
		res_ts <- readRDS(ftempin)

		res_dims <- dim(res_ts)
		res_dims[2] <- length(mmonths)
		res_dims <- c(res_dims, 2)
		res_dimnames <- dimnames(res_ts)
		res_dimnames[[2]] <- month.abb
		res_dimnames[[length(res_dimnames) + 1]] <- c("mean", "sd")
		res_mon <- array(NA, dim = res_dims, dimnames = res_dimnames)
	
		temp <- apply(res_ts,
			MARGIN = seq_along(dim(res_ts))[-2],
			function(x) tapply(x, rep(1:12, nyears), mean, na.rm = TRUE))
		res_mon[, , , , , "mean"] <- aperm(temp, c(2:1, 3:length(dim(res_ts))))
		
		temp <- apply(res_ts,
			MARGIN = seq_along(dim(res_ts))[-2],
			function(x) tapply(x, rep(1:12, nyears), sd, na.rm = TRUE))
		res_mon[, , , , , "sd"] <- aperm(temp, c(2:1, 3:length(dim(res_ts))))

		res_mon[!is.finite(res_mon)] <- NA
		saveRDS(res_mon, file = ftemp)
		
		rm(res_ts, res_mon)
	}
}

# Metrics of model performance
rmse <- function(x, y, na.rm = FALSE) sqrt(mean((x - y) ^ 2, na.rm = na.rm))

mae <- function(x, y, na.rm = FALSE) mean(abs(x - y), na.rm = na.rm)

r2lin <- function(obs, pred, na.action = na.omit) {
	# Pineiro, G., Perelman, S., Guerschman, J.P. & Paruelo, J.M. (2008). How to evaluate models: Observed vs. predicted or predicted vs. observed? Ecol Model, 216, 316-322.
  z <- lm(obs ~ pred, na.action = na.action)
  r <- z$residuals
  f <- z$fitted.values
	
	mss <- if (attr(z$terms, "intercept")) sum((f - mean(f))^2) else sum(f^2)
  rss <- sum(r^2)
	mss / (mss + rss)
}

AC <- function(x, y, na.rm = FALSE) {
	# Riemann, R., Wilson, B.T., Lister, A. & Parks, S. (2010). An effective assessment protocol for continuous geospatial datasets of forest characteristics using USFS Forest Inventory and Analysis (FIA) data. Remote Sensing of Environment, 114, 2337-2352.
	
	x_mean <- mean(x, na.rm = na.rm)
	y_mean <- mean(y, na.rm = na.rm)
	slope <- sqrt(sum((y - y_mean) ^ 2, na.rm = na.rm) / sum((x - x_mean) ^ 2, na.rm = na.rm))
	int <- y_mean - slope * x_mean
	
	x_pred <- (y - int) / slope
	y_pred <- int + slope * x
	
	SSD <- sum((x - y) ^ 2, na.rm = na.rm) # sum of squared difference
	SPDu <- sum((x - x_pred) * (y - y_pred), na.rm = na.rm)	# unsystematic sum of product difference
	SPDs <- SSD - SPDu
	
	SPOD <- sum((abs(x_mean - y_mean) + abs(x - x_mean)) * (abs(x_mean - y_mean) + abs(y - y_mean)), na.rm = na.rm)	# sum of potential difference
	
	# ACsys = 1 if the GMFR line is perfectly in line with the 1:1 line
	# ACuns = 1 if all points fall directly on the GMFR line.
	c(AC = 1 - SSD / SPOD, ACsys = 1 - SPDs / SPOD, ACuns = 1 - SPDu / SPOD)
}

Duveillers_lambda <- function(x, y, na.rm = FALSE) {
	# Duveiller, G., Fasbender, D. & Meroni, M. (2016). Revisiting the concept of a symmetric index of agreement for continuous datasets. Scientific Reports, 6, 19401.
	
	if (all(is.na(x)) || all(is.na(y))) {
		NA
	} else {
		x_mean <- mean(x, na.rm = na.rm)
		y_mean <- mean(y, na.rm = na.rm)
	
		r <- cor(x, y, method = "pearson", use = if (na.rm) "pairwise.complete.obs" else "everything")
		k <- if (r >= 0) 0 else {2 * abs(sum((x - x_mean) * (y - y_mean), na.rm = na.rm))}	# eq. 14
		delta <- mean((x - y) ^ 2, na.rm = na.rm)
		mu <- var(x, na.rm = na.rm) + var(y, na.rm = na.rm) + (x_mean - y_mean) ^ 2 + k
		lambda <- 1 - delta / mu	# eqs 3 and 15
	}
}

agg_for_region <- function(i, data, iuse) {
	res <- rep(NA, 12)
	i_reg <- dLoc$Region[iuse] == i
	
	if (sum(i_reg) > 0) {
		dat <- data[i_reg, , drop = FALSE]
		if (any(!is.na(dat))) {
			res <- scale(apply(dat, 2, mean, na.rm = TRUE), scale = TRUE, center = TRUE)
		}
	}
	
	res
}

if (do_compare) {
	measures <- c("rmse", "Duveillers_lambda")
	
	#---Simulation soil moisture
	load(file.path(dir.gtd, "6_Projects_Year1", "Prj03_GlobalVulnerability", "1_PC_TempDry_Simulations_Prj03_r2", "4_Data_SWOutputAggregated", "DATA_v4", "Extraction", "dLoc_dbTables_exp01.RData"))
	load(file.path(dir.gtd, "6_Projects_Year1", "Prj03_GlobalVulnerability", "1_PC_TempDry_Simulations_Prj03_r2", "4_Data_SWOutputAggregated", "DATA_v4", "PreCalculations", "dStudy2_dbTables_exp01.RData"))
	load(file.path(dir.gtd, "6_Projects_Year1", "Prj03_GlobalVulnerability", "1_PC_TempDry_Simulations_Prj03_r2", "4_Data_SWOutputAggregated", "DATA_v4", "PreCalculations", "resSoils2_dbTables_exp01.RData"))
	sdepth_gtd_m <- resSoils2["Simulation", currentSc, currentSc, , "SWinput_Soil_maxDepth_cm"] / 100
	rm(resSoils2)
	#i_Any33Cond <- apply(dStudy2["MetDef_ThisCond", , , ], 3, function(x) any(x == 1L, na.rm = TRUE))
	
	load(file.path(dir.gtd, "6_Projects_Year1", "Prj03_GlobalVulnerability", "1_PC_TempDry_Simulations_Prj03_r2", "4_Data_SWOutputAggregated", "DATA_v4", "PreCalculations", "resVWC2_dbTables_exp01.RData"))
	temp <- dimnames(resVWC2)
	i_RCP_gtd <- temp[[2]]
	i_GCM_gtd <- temp[[3]]
	temp <- strsplit(temp[[5]], "_")
	i_months_gtd <- sapply(temp, function(x) x[3])
	i_lyr_gtd <- sapply(temp, function(x) x[2])
	
	# Sort along data extracted from GCMs
	id_gcm <- apply(xy_Any33Cond, 1, paste, collapse = "_")
	id_gtd <- apply(dLoc[, c("X_WGS84", "Y_WGS84")], 1, paste, collapse = "_")
	imatch_gcm_to_gtd <- match(id_gtd, id_gcm, nomatch = 0)
	iuse_gtd <- imatch_gcm_to_gtd > 0 # identical(iuse_gtd, i_Any33Cond) == TRUE
	
	
	#---GCM soil moisture
	if (sum(iuse_gtd) > 0) for (iv in seq_along(vars)) {
		fout1 <- file.path(dir.out, paste0("ComparisonTable_", vars[iv], "_meanmonthly_", version, ".rds"))
		fout2 <- file.path(dir.out, paste0("ComparisonGrids_", vars[iv], "_meanmonthly_", version, ".rds"))
		if (file.exists(fout1) && file.exists(fout2)) next

		fin <- file.path(dir.ext, paste0("Extract_", vars[iv], "_meanmonthly_", version, ".rds"))
		if (!file.exists(fin)) next
	
		resGCM_mmon <- readRDS(fin)
		temp <- dimnames(resGCM_mmon)
		i_RCP_gcm <- temp[[3]]
		i_GCM_gcm <- temp[[4]]
		i_months_gcm <- temp[[2]]
		i_lyr_gcm <- seq_len(dim(resGCM_mmon)[5])

		load(file = sub("_meanmonthly_", "_sdepth_", fin)) # sdepths_m, var_soil_N
	
		#---Comparison
		i_RCP_int <- unique(c(currentSc, intersect(i_RCP_gtd, i_RCP_gcm)))
		i_GCM_int <- intersect(i_GCM_gtd, i_GCM_gcm)
		iuse_lyr_gtd <- switch(EXPR = vars[iv],
							mrso = grepl("allLayers", i_lyr_gtd),
							mrsos = grepl("topLayers", i_lyr_gtd),
							mrlsl = rep(TRUE, length(i_lyr_gtd)))
		
		#---Output containers
		grid_comp <- array(list(),
							dim = c(length(i_RCP_int), length(i_GCM_int), length(measures)),
							dimnames = list(i_RCP_int, i_GCM_int, measures))
		res_comp <- list(cells = array(NA,
							dim = c(length(i_RCP_int), length(i_GCM_int), sum(iuse_gtd), length(measures)),
							dimnames = list(i_RCP_int, i_GCM_int, NULL, measures)),
						regions = array(NA,
							dim = c(length(i_RCP_int), length(i_GCM_int), 1 + length(regions_N), length(measures)),
							dimnames = list(i_RCP_int, i_GCM_int, c("gtd", label.regions), measures)))
		grid_template <- mask_Any33Cond
		grid_template[] <- NA		

		for (ir in seq_along(i_RCP_int)) {
			ir_gcm <- ir_gtd <- i_RCP_int[ir]
			if (i_RCP_int[ir] == currentSc) ir_gcm <- "historical"
		
			for (ig in seq_along(i_GCM_int)) {
				ig_gcm <- ig_gtd <- i_GCM_int[ig]
				if (i_RCP_int[ir] == currentSc) ig_gtd <- currentSc

				d_gcm <- resGCM_mmon[imatch_gcm_to_gtd, month.abb, ir_gcm, ig_gcm, i_lyr_gcm, "mean"]
				d_gtd <- resVWC2["Simulation", ir_gtd, ig_gtd, iuse_gtd, iuse_lyr_gtd]
				
				sdepth_m_gcm <- unique(sdepths_m[grepl(ir_gcm, names(sdepths_m), ignore.case = TRUE) & grepl(ig_gcm, names(sdepths_m), ignore.case = TRUE)])
				if (length(sdepth_m_gcm) == 1) {
					sdepth_m_gcm <- sdepth_m_gcm[[1]]
				} else {
					stop(paste("Soil depth not unique for", vars[iv], ir_gcm, ig_gcm))
				}
				
				if (FALSE) {
					tmask <- grid_template
					tmask[raster::cellFromXY(tmask, dLoc[, c("X_WGS84", "Y_WGS84")])] <- d_gtd
					#tmask[raster::cellFromXY(tmask, dLoc[, c("X_WGS84", "Y_WGS84")])] <- 1
					plot(tmask)
					map(add = TRUE)

					tmask <- grid_template
					tmask[raster::cellFromXY(tmask, xy_Any33Cond)] <- d_gcm
					#tmask[raster::cellFromXY(tmask, xy_Any33Cond)] <- 1
					#tmask[raster::cellFromXY(tmask, xy_Any33Cond)] <- res_ts[, 1, "historical", 5, 1]
					plot(tmask)
					map(add = TRUE)
				}
			
				if (sum(!is.na(d_gcm)) > 0 && sum(!is.na(d_gtd)) > 0) {
					# Per Cell
					ds_gcm <- t(apply(d_gcm, 1, scale, scale = TRUE, center = TRUE))
					ds_gtd <- t(apply(d_gtd, 1, scale, scale = TRUE, center = TRUE))
					ds_cell <- cbind(ds_gcm, ds_gtd)
					ds_used <- complete.cases(ds_cell)
					ds_cell <- ds_cell[ds_used, ]

					# Per Region and global
					ds_gcm_reg <- lapply(regions_N, agg_for_region, data = d_gcm, iuse = iuse_gtd)
					ds_gtd_reg <- lapply(regions_N, agg_for_region, data = d_gtd, iuse = iuse_gtd)
					ds_gcm_glo <- scale(apply(d_gcm, 2, mean, na.rm = TRUE), scale = TRUE, center = TRUE)
					ds_gtd_glo <- scale(apply(d_gtd, 2, mean, na.rm = TRUE), scale = TRUE, center = TRUE)
					
					# Calculate along measures
					for (mf in measures) {
						mfun <- match.fun(mf)
						
						# Measure per cell
						res_comp[["cells"]][ir_gtd, ig_gcm, ds_used, mf] <- apply(ds_cell, 1,
							function(xy) mfun(xy[1:12], xy[13:24], na.rm = TRUE))

						perf1 <- grid_template
						grid_cells <- raster::cellFromXY(perf1, dLoc[which(iuse_gtd)[ds_used], c("X_WGS84", "Y_WGS84")])
						perf1[grid_cells] <- res_comp[["cells"]][ir_gtd, ig_gcm, ds_used, mf]
						grid_comp[ir_gtd, ig_gcm, mf] <- list(perf1)

						if (FALSE) {
							plot(perf1, zlim = c(0, 1))
							title(main = paste0("lambda-agreement of mean monthly soil moisture for ", paste(range(years[[if (ir < 2) 1 else 2]]), collapse = "-"), ":\nGCM (", ig_gcm, ") vs. Soilwat (", ig_gtd, ")"))
							map(add = TRUE)
						}

						# Measure per region
						res_comp[["regions"]][ir_gtd, ig_gcm, label.regions, mf] <- sapply(regions_N,
							function(i) mfun(ds_gcm_reg[[i]], ds_gtd_reg[[i]], na.rm = TRUE))
						
						# Measure globally
						res_comp[["regions"]][ir_gtd, ig_gcm, "gtd", mf] <- 
							mfun(ds_gcm_glo, ds_gtd_glo, na.rm = TRUE)
					}
				}
			}
		}
			
		saveRDS(res_comp, file = fout1)
		saveRDS(grid_comp, file = fout2)
	}
}	



if (do_tables) {
	for (iv in seq_along(vars)) {
		fout1 <- file.path(dir.out, paste0("ComparisonTable_", vars[iv], "_meanmonthly_", version, ".rds"))
		if (!file.exists(fout1)) next
		
		reg_comp <- readRDS(fout1)[["regions"]]
		
		reg_comp_table <- dcast(melt(reg_comp), Var4 + Var1 + Var2 ~ Var3)
		
		write.csv(reg_comp_table, file = sub(".rds", ".csv", fout1), row.names = FALSE)
	}
}


if (do_maps) {
	for (iv in seq_along(vars)) {
		fout2 <- file.path(dir.out, paste0("ComparisonGrids_", vars[iv], "_meanmonthly_", version, ".rds"))
		if (!file.exists(fout2)) next

		grid_comp <- readRDS(fout2)
		
		temp <- dimnames(grid_comp)
		i_RCP_int <- temp[[1]]
		i_GCM_int <- temp[[2]]
		measures <- temp[[3]]
		
		xlim <- c(-125, 132)
		ylim <- c(-55, 60)
		cols <- colorRampPalette(c("orange", "cornflowerblue", "darkblue"))(255)
		cover_north <- raster::raster(raster::extent(-180, 180, ylim[2], 90), crs = WGS84)
		cover_north[] <- 1
		cover_south <- raster::raster(raster::extent(-180, 180, -90, ylim[1]), crs = WGS84)
		cover_south[] <- 1
		
		for (im in seq_along(measures)) {
			fmap <- sub(paste0(version, ".rds"), paste0(measures[im], "_", version, ".pdf"), fout2)
			
			zlim <- switch(EXPR = measures[im],
						Duveillers_lambda = c(0, 1),
						AC = c(0, 1),
						r2lin = c(0, 1),
						NULL)
			zlimc <- NULL
			
			h.panel <- 2; w.panel <- h.panel / diff(ylim) * diff(xlim)
	s		w.edgeL <- 0.4; w.edgeR <- 0.1; h.edgeL <- 0.4
			cex <- 1
			nrows <- length(i_GCM_int); ncols <- length(i_RCP_int)
			
			pdf(height = h.edgeL + h.panel * nrows, width = w.edgeL + w.panel * ncols + w.edgeR, file = fmap)
			layout(
				matrix(c(rep(0, 1 + nrows), sapply(seq_len(ncols), function(x) c((x - 1) * nrows + seq_len(nrows), 0)), rep(0, 1 + nrows)),
					nrow = 1 + nrows,
					ncol = 1 + ncols + 1,
					byrow = FALSE),
				heights = c(rep(h.panel, times = nrows), h.edgeL),
				widths = c(w.edgeL, rep(w.panel, times = ncols), w.edgeR))
			par_old <- par(mgp = c(1, 0, 0), mar = rep(0.5, 4), tcl = 0.3, cex = cex)

			i_panel <- 1
			legend_added <- FALSE
			for (ir in seq_along(i_RCP_int)) for (ig in seq_along(i_GCM_int)) {
				x <- grid_comp[i_RCP_int[ir], i_GCM_int[ig], measures[im]][[1]]
				if (is.null(zlimc)) {
					zlimc <- if (is.null(zlim)) {
							c(raster::cellStats(x, "min", na.rm = TRUE), raster::cellStats(x, "max", na.rm = TRUE))
						} else {
							zlim
						}
				}
				
				if (!is.null(x) && raster::ncell(x) > raster::cellStats(x, 'countNA')) {
					 raster::image(x,
						maxpixels = raster::ncell(x),
						col = cols,
						xlim = xlim, ylim = ylim, zlim = zlimc, 
						xlab = "", ylab = "",
						asp = 1, cex = cex, axes = FALSE)
					
					map("world", lwd = 0.5, col = "gray20", add = TRUE)
					if (!is.null(cover_north))
						raster::image(cover_north, col = "white", add = TRUE)
					if (!is.null(cover_south))
						raster::image(cover_south, col = "white", add = TRUE)

					ats <- axTicks(1)
					axis(1, pos = ylim[1], at = ats, labels = if (ig == length(i_GCM_int)) ats else FALSE, xpd = FALSE)
					ats <- c(-90, (ats <- axTicks(2))[ats >= -90 & ats <= 90], 90)
					axis(2, at = ats, labels = if (ir == 1) ats else FALSE, xpd = FALSE)
					if (ir == 1)
						mtext(side = 2, line = 1.2, "Latitude", cex = par("cex") * cex)
					if (ig == length(i_GCM_int))
						mtext(side = 1, line = 1, "Longitude", cex = par("cex") * cex)
					mtext(side = 3, line = -0.75, adj = 0.025, cex = cex, font = 2, text = tolower(letters[i_panel]))
					
					if (!legend_added) {
													
						col_desc <- make_colors(zlim = zlimc, zextremes = zlimc,
							var_sign = 1,
							val_crit = -Inf,
							colors_Pos = NULL,
							colors_Neg = NULL,
							colors_nozero = rev(cols),
							colors_zero = NULL,
							colors_below = NULL,
							colors_above = NULL,
							n_colors = 255)
						
						add_legend(zlimc, zlimc, col_desc, grid = x,
							box = c(-30, 100, ylim[1] + 2, ylim[1] + 10), horiz = TRUE,
							srt = 0, cex = 0.75 * cex)
						legend_added <- TRUE
					}
					
				} else {
					plot.new()
				}
				
				i_panel <- i_panel + 1
			}
						
			par(par_old)
			dev.off()
		}
	}
}