5f_GTD_Prj03v4_Analysis_TablesSummary.R

################### FORMAT TABLES FOR MANUSCRIPT ##################
#
# Daniel Schlaepfer, 2015-2016
# 
# Format and save tables for manuscript for responses in temperate dryland areas with
# all study_areas_and_shifts x RCP combinations
#	- requests 'get_precalculations' from '5a1_GTD_*.R'
#	- requests 'get_studyareaextents' from '5a1_GTD_*.R'

# Note: aggregating changes per region may not deliver desired/expected values
#	- current value = the median over current extent
#	- future change values = medians over future extent (i.e., stable + leading zones)
#	 => median(future) != median(current) + median(change)
# ==> Use absolute values instead or aggregate by region x shifts
###################################################################


#-------------------------------
#---GLOBAL SETTINGS
comp <- "dropbox"
do_tables <- TRUE
redo_tables <- FALSE

#-------------------------------
#---R packages
pkg_reqd <- c("base")
has_loaded <- sapply(pkg_reqd,
	function(lib) require(lib, character.only = TRUE, quietly = FALSE))
stopifnot(has_loaded)

#---Load data and misc. functions
if (comp %in% c("err", "eleos")) {
	dir.gtd <- "/PATH_TO_PROJECT/Product_PowellCenter/6_Projects_Year1"
} else if (comp == "dropbox") {
	dir.gtd <- "/PATH_TO_PROJECT/Product_PowellCenter/6_Projects_Year1"
}
dir.prj <- file.path(dir.gtd, "Prj03_GlobalVulnerability", "4_Analysis", "4_Analysis_v4")


get_precalculations <- TRUE
get_studyareaextents <- TRUE
source(file.path(dir.prj, "5a1_GTD_Prj03v4_Helper.R"))

stopifnot(done_precalculations, done_studyareaextents)


#---Directories
dir.create(dir.fig_SAM <- file.path(dir.prj, "6_Results", "5_Tables_v4"), showWarnings = FALSE)


#-------------------------------
#------TABLES AND FIGURES

if (do_tables) {
	res_names <- sets_names_extracted2
	res_names_cor <- sets_names_cor2
	type_list_mode <- change_types
	
	do_tabulate_regional_response <- function(dat, dat_current = NULL, shift_used = FALSE, val_type = change_types, varnames, var_sign = -1, probs_current = c(0.025, 0.975), dir_out, ftag = varname) {
#str(dat_current): num [1:6, 1:3, 1:20021, 1:nvars] 1002 NA NA 987 NA ...
# - attr(*, "dimnames") = List of 3
# ..$ : chr [1:3] "Simulation" "MetDef_Any33Cond" "MetDef_ThisCond"
# ..$ : chr [1:2] "Current", "RCP45" "RCP85"
# ..$ : NULL
# ..$ : chr [1:nvars] vars
#str(data): num [1:6, 1:2, 1:16, 1:20021, 1:nvars] 1074 NA NA 1134 NA ...
# - attr(*, "dimnames") = List of 5
# ..$ : chr [1:3] "Simulation" "MetDef_Any33Cond" "MetDef_ThisCond"
# ..$ : chr [1:2] "RCP45" "RCP85"
# ..$ : chr [1:16] "CanESM2" "CESM1-CAM5" "CSIRO-Mk3-6-0" "EC-EARTH" ...
# ..$ : NULL
# ..$ : chr [1:nvars] vars
		inc <- function(var) {
			assign(var, get(var) + 1, inherits = TRUE)
			return(get(var))
		}

		val_type <- match.arg(val_type)
		ftag <- clean_name(ftag)
		
		if (shift_used) stopifnot(exists("dShift2_rankpercell"))
		stopifnot(exists("dShift2"))
		
		ftemp <- file.path(dir_out, fname <- paste0("Table_", ftag, "_", val_type, if (shift_used) "_shifts", ".csv"))
		if (redo_tables || !file.exists(ftemp)) {
			#---Prepare table
			n_resp <- 2 # 2 = length(list(median, range of ranks))
			ncol_per_rcp <- n_resp * if (shift_used) length(shifts) else 1
			id_region <- dLoc$Region
			i_region <- (id_region %in% regions_n)
			if (!(length(var_sign) == length(varnames))) var_sign <- rep(var_sign[1], times = length(varnames))
			if (any(is.na(var_sign))) var_sign[is.na(var_sign)] <- 1
			
			# Data
			has_current_values <- !is.null(dat_current["MetDef_ThisCond", currentSc, , ]) && !(sd(dat_current["MetDef_ThisCond", currentSc, , ], na.rm = TRUE) == 0)
			var_used <- varnames

			
			#---Table
			if (shift_used) {
				ncol_per_var <- n_resp * length(currentSc) + ncol_per_rcp * (length(currentSc) + length(reqRCPs))
				dimn2 <- paste(
							rep(clean_name(var_used), each = ncol_per_var),
							
							c(
								paste0(
									rep(c(currentSc, paste0(currentSc, "_", round(100 * diff(range(probs_current))), "perc_amongCells")), times = 1 + length(shifts)),
									c(rep("", each = n_resp), rep(paste0("_RCP85rank8_", shifts), each = n_resp))
								),
								paste(
									rep(reqRCPs, each = ncol_per_rcp),	
									paste(
										rep(c(paste0("rank", ranks[2]), paste0("rank", ranks[1], "TOrank", ranks[3])), times = length(shifts)),
										rep(shifts, each = n_resp),
									sep = "_"),
								sep = "_")
							),
						sep = "_")
			} else {
				ncol_per_var <- n_resp * length(currentSc) + ncol_per_rcp * (length(reqRCPs))
				dimn2 <- paste(rep(clean_name(var_used), each = ncol_per_var),
									c(currentSc,
									paste0(currentSc, "_", round(100 * diff(range(probs_current))), "perc_amongCells"),
									paste(rep(reqRCPs, each = n_resp),
									c(paste0("rank", ranks[2]), paste0("rank", ranks[1], "TOrank", ranks[3])),
								sep = "_")), sep = "_")
			}
			
			#Table per GCM
			#resSc <- array(NA, dim = c(1 + length(regions_n), length(var_used), length(reqRCPs), length(reqGCMs)), dimnames = list(c("Total", label.regions[regions_n]), var_used, reqRCPs, reqGCMs))
			
			#Final table
			resT <- data.frame(matrix(NA, nrow = 1 + length(regions_n), ncol = ncol_per_var * length(var_used),
					dimnames = list(c("Total", label.regions[regions_n]), dimn2)))
			
			
			iv_add <- 1
			irow <- 0
			for (iv in seq_along(var_used)) {
				#Prepare data
				is_circular <- grepl("_doy", var_used[iv])
				
				if (is_circular) {
					f_quantile <- function(x, ...) circ.quantile(x, int = 365, modulo = "2pi", ...)
					f_median <- function(x, ...) circ.median(x, int = 365, ...)
				} else {
					f_quantile <- quantile
					f_median <- median
				}
				
				if (has_current_values) {
					vdat_current <- dat_current["MetDef_ThisCond", currentSc, , var_used[iv]]
					# Sweep with current data at 'simulation' extent because otherwise 'leading' zone will not get sweeped
					vdat <- dat_sweep(dat[, , , , var_used[iv]], dat_current["Simulation", currentSc, , var_used[iv]], val_type, circular = is_circular)
				} else {
					vdat <- NULL
				}
				
				#Shifts
				seq_shifts <- if (shift_used) seq_along(shifts) else 1
				
				# Current
				if (abs(max(vdat_current, na.rm = TRUE)) >= 50) {
					round_d <- 0
					round_fun <- round
				} else if (abs(max(vdat_current, na.rm = TRUE)) >= 10) {
					round_d <- 1
					round_fun <- round
				} else {
					round_d <- 2
					round_fun <- signif
				}
				
				#	- overall current
				temp <- round_fun(do.call(cbind, tapply(vdat_current, INDEX = id_region, function(x)
										c(f_median(x, na.rm = TRUE), f_quantile(x, probs = probs_current, na.rm = TRUE))))[,regions_n], round_d)
				temp_tot <- round_fun(c(f_median(vdat_current[i_region], na.rm = TRUE), 
										f_quantile(vdat_current[i_region], probs = probs_current, na.rm = TRUE)), round_d)
				
				resT[, inc("irow")] <- c(temp_tot[1], temp[1,])
				resT[, inc("irow")] <- c(paste0("(", temp_tot[2], ", ", temp_tot[3], ")"), paste0("(", temp[2,], ", ", temp[3,], ")"))


				#	- current by ranked shifting zones
				if (shift_used) {
					for (iss in seq_shifts) {
						ids <- !is.na(dShift2_rankpercell[shifts[iss], "RCP85", 2, ])
					
						temp1 <- do.call(cbind, tapply(dat_current["Simulation", currentSc, ids, var_used[iv]], INDEX = id_region[ids], function(x)
												c(f_median(x, na.rm = TRUE), f_quantile(x, probs = probs_current, na.rm = TRUE))))
					
						temp <- round_fun(temp1[, match(regions_n, colnames(temp1))], round_d)
						temp_tot <- round_fun(c(f_median(dat_current["Simulation", currentSc, ids & i_region, var_used[iv]], na.rm = TRUE), 
												f_quantile(dat_current["Simulation", currentSc, ids & i_region, var_used[iv]], probs = probs_current, na.rm = TRUE)), round_d)
				
						resT[, inc("irow")] <- c(temp_tot[1], temp[1,])
						resT[, inc("irow")] <- c(paste0("(", temp_tot[2], ", ", temp_tot[3], ")"), paste0("(", temp[2,], ", ", temp[3,], ")"))
						
					}
				}
				
			
				for (ircp in seq_along(reqRCPs)) {
					# RCPs
					if (abs(f_quantile(vdat[shifts, reqRCPs[ircp], , ], prob = 0.75, na.rm = TRUE)) >= 100) {
						round_d <- 0
						round_fun <- round
					} else if (abs(f_quantile(vdat[shifts, reqRCPs[ircp], , ], prob = 0.75, na.rm = TRUE)) >= 10) {
						round_d <- 1
						round_fun <- round
					} else {
						round_d <- 2
						round_fun <- signif
					}				

					for (iss in seq_shifts) {	
						i_area <- if (shift_used) shifts[iss] else "MetDef_ThisCond"
					
						# get ranks: first aggregate over GCMs, then rank them
						rtemp1 <- matrix(NA, nrow = 1 + length(regions_n), ncol = length(reqGCMs))					
						for (igcm in seq_along(reqGCMs)) {
							#ids <- ZeroToNA(as.integer(apply(dShift2[ircp, igcm, , ][, i_area, drop = FALSE], 1, function(x) any(!is.na(x)))))
							
							rtemp1[-1, igcm] <- tapply(vdat[i_area, ircp, igcm, ], INDEX = id_region, f_median, na.rm = TRUE)[regions_n]
							rtemp1[1, igcm] <- f_median(vdat[i_area, ircp, igcm, ], na.rm = TRUE)
						}
						
						rtemp2 <- round_fun(t(apply(rtemp1, 1, calc_ensemble4)), round_d)

						
						# rank 8
						resT[, inc("irow")] <- rtemp2[, 2]
				
						# rank spread
						resT[, inc("irow")] <- paste0("(", rtemp2[, 1], ", ", rtemp2[, 3], ")")
					}
				}
			}

			write.csv(t(resT), file = ftemp)
		}

		invisible(NULL)
	}



	print("DrylandResponse: regional tables")
	for (id in seq_along(res_names)) {
		dir.create(dir.temp <- file.path(dir.fig_SAM, res_names[id]), showWarnings = FALSE)
	
		# get data objects
		temp <- fix_var_dim(get(res_names[id]))
		data <- temp$data[, -1, -1, , ]
		vars <- temp$vars
		data_current <- temp$data[study_areas_all, , currentSc, , ]
		if (length(dim(data_current)) <= 3) {
			dimns <- dimnames(data_current)
			dim(data_current) <- c(dim(data_current), 1)
			dimnames(data_current) <- modifyList(dimns, list(vars = vars))
		}
		rm(temp)
	
		# obtain correlation with AI
		corAIs <- if (is.character(res_names_cor[id][[1]])) get(res_names_cor[id][[1]])$median_tau["median", , "MetDef_ThisCond"] else 1

		for (it in seq_along(type_list_mode)) {
			for (iss in c(FALSE, TRUE)) { #shift_used
				print(paste(Sys.time(), res_names[id], type_list_mode[it], "shift_used", iss))
			
				do_tabulate_regional_response(dat = data,
							dat_current = data_current,
							shift_used = iss,
							val_type = type_list_mode[it],
							varnames = vars,
							var_sign = corAIs,
							dir_out = dir.temp, ftag = res_names[id])
			}
		}

	}
}