-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy path5e_GTD_Prj03v2_Analysis_Scatterplots.R
executable file
·191 lines (150 loc) · 9.36 KB
/
5e_GTD_Prj03v2_Analysis_Scatterplots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
################### PLOTS SCATTERS BETWEEN RESPONSE VARIABLES #####
#
# Daniel Schlaepfer, 2015-2016
#
# Plots scatters between response variables
# - requests 'get_precalculations' from '5a1_GTD_*.R'
# - requests 'get_studyareaextents' from '5a1_GTD_*.R'
#
###################################################################
#-------------------------------
#---GLOBAL SETTINGS
comp <- "dropbox"
do_scatter <- TRUE
redo_scatter <- FALSE
#-------------------------------
#---R packages
pkg_reqd <- c("MASS", "digest")
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", "3_Analysis_v3")
get_precalculations <- TRUE
get_studyareaextents <- TRUE
source(file.path(dir.prj, "5a1_GTD_Prj03v2_Helper.R"))
stopifnot(done_precalculations, done_studyareaextents)
#---Directories
dir.create(dir.fig_SAM <- file.path(dir.prj, "6_Results", "4_Scatterplots_v4"), showWarnings = FALSE)
#-------------------------------
#------Functions
add_shifts_to_current_dat_v3 <- function(dat, dat_current) {
res <- array(NA, dim = dim(dat), dimnames = dimnames(dat))
# Copy existing extents
for (i1 in 1:dim(dat)[[1]]) {
if (dimnames(dat)[[1]][i1] %in% dimnames(dat_current)[[1]]) {
# Repeat data for RCPs and GCMs/ranks
for (i2 in 1:dim(dat)[[2]]) for (i3 in 1:dim(dat)[[3]]) res[i1, i2, i3, ] <- dat_current[i1, ]
}
}
# Extract data for shifts
for (iss in seq_along(shifts)) for (i2 in 1:dim(dat)[[2]]) for (i3 in 1:dim(dat)[[3]]) {
ishift <- !is.na(dat[shifts[iss], i2, i3, ]) #if it has values then in shifting zone (except few response variables that have NAs...
res[shifts[iss], i2, i3, ishift] <- dat_current["MetDef_Any33Cond", ishift]
}
return(res)
}
add_shifts_to_current_dat_v4 <- function(dat, dat_current) {
res <- array(NA, dim = dim(dat), dimnames = dimnames(dat))
# Copy existing extents
for (i1 in 1:dim(dat)[[1]]) {
if (dimnames(dat)[[1]][i1] %in% dimnames(dat_current)[[1]]) {
# Repeat data for RCPs and GCMs/ranks
for (i2 in 1:dim(dat)[[2]]) for (i3 in 1:dim(dat)[[3]]) res[i1, i2, i3, ] <- dat_current[i1, ]
}
}
# Extract data for shifts
for (iss in seq_along(shifts)) for (i2 in 1:dim(dat)[[2]]) for (i3 in 1:dim(dat)[[3]]) {
ishift <- !is.na(dat[shifts[iss], i2, i3, ]) #if it has values then in shifting zone (except few response variables that have NAs...
res[shifts[iss], i2, i3, ishift] <- dat_current["MetDef_Any17Cond", ishift]
}
return(res)
}
#-------------------------------
#------TABLES AND FIGURES
if (do_scatter) {
res_names <- sets_names_extracted2
res_names_ranks <- sets_names_ranks2
res_labels <- list(label_climate, label_response, label_definition, label_soils)
type_list_mode <- change_types[1:2]
do_scatter_regional_shifted_response <- function(datX, datY, datX_current, datY_current, axisX_type = c("Current", "Future"), valX_type, valY_type, varnameX, varnameY, labelX = varnameX, labelY = varnameY, val_critX = 0, val_critY = 0, probsX = c(0.005, 0.995), probsY = c(0.005, 0.995), dir_out, ftagX = varnameX, ftagY = varnameY) {
ftagX <- clean_name(ftagX)
ftagY <- clean_name(ftagY)
fname <- paste0("Scatter_", ftagY, "_", valY_type, "_by_", ftagX, "_", valX_type, "_", axisX_type, "Xaxis.pdf")
if (nchar(fname) >= 255) {
tempY <- if (nchar(ftagY) <= 50) ftagY else paste0(substr(ftagY, 1, 25), "-", substr(ftagY, 26, 50))
tempX <- if (nchar(ftagX) <= 50) ftagY else paste0(substr(ftagX, 1, 25), "-", substr(ftagX, 26, 50))
fname <- paste0("Scatter_", tempY, "-", digest(ftagY), "_", valY_type, "_by_", tempX, "-", digest(ftagX), "_", valX_type, "_", axisX_type, "Xaxis.pdf")
}
ftemp <- file.path(dir_out, fname)
if (redo_scatter || !file.exists(ftemp)) {
#---Prepare plot
# Data
if (axisX_type == "Future") {
datX <- dat_sweep(datX, datX_current["MetDef_Any33Cond", ], valX_type, circular = grepl("_doy", varnameX))
} else if (axisX_type == "Current") {
datX <- add_shifts_to_current_dat_v4(dat = datX, dat_current = datX_current)
}
datY <- dat_sweep(datY, datY_current["MetDef_Any33Cond", ], valY_type, circular = grepl("_doy", varnameY))
# Range limits
zlimitsX <- calc_response_limits(dat_cur = datX_current["MetDef_ThisCond", ], dat = datX, get_ircps = NULL, get_ranks = NULL, Any33Cond = FALSE, type = valX_type, probs = probsX)
zlimitsY <- calc_response_limits(dat_cur = datY_current["MetDef_ThisCond", ], dat = datY, get_ircps = NULL, get_ranks = NULL, Any33Cond = FALSE, type = valY_type, probs = probsY)
#---Plot
scatter_regional_response(dat_curX = datX_current["MetDef_ThisCond", ], dat_curY = datY_current["MetDef_ThisCond", ], dat_ensX = datX, dat_ensY = datY,
xlab1 = labelX, xlab2 = valX_type, ylab1 = labelY, ylab2 = valY_type,
xlims = zlimitsX$zlims, xextremes = zlimitsX$zextremes, ylims = zlimitsY$zlims, yextremes = zlimitsY$zextremes,
val_critX = val_critX, val_critY = val_critY,
dir.out = dir_out, fname = fname)
}
invisible(NULL)
}
print("DrylandResponse: regional scatter plots")
# get data objects
scatter_list <- list(list(idx = 1, idy = 1, ivx = 2, ivy = 1, ixas = 2), #MAP vs MAT
list(idx = 1, idy = 1, ivx = 4, ivy = 1, ixas = 2), #MAP vs Seasonality
list(idx = 1, idy = 4, ivx = 1, ivy = 2, ixas = 2), #Sand0 vs MAP
list(idx = 1, idy = 4, ivx = 2, ivy = 2, ixas = 2), #Sand0 vs MAT
list(idx = 1, idy = 2, ivx = 1, ivy = 8, ixas = 2), #TranspirationBottomToTranspirationTotal_fraction vs MAP
list(idx = 1, idy = 2, ivx = 4, ivy = 8, ixas = 2), #TranspirationBottomToTranspirationTotal_fraction vs Seasonality
list(idx = 2, idy = 2, ivx = 1, ivy = 2, ixas = 2), #ThermalSnowfreeDryPeriods_SWPcrit3000kPa_bottomLayers_DrySpellsAllLayers_maxDuration_days vs ThermalSnowfreeDryPeriods_SWPcrit3000kPa_topLayers_DrySpellsAllLayers_maxDuration_days
list(idx = 2, idy = 2, ivx = 8, ivy = 8, ixas = 1:2), #TranspirationBottomToTranspirationTotal_fraction vs TranspirationBottomToTranspirationTotal_fraction
list(idx = 2, idy = 2, ivx = 1, ivy = 27, ixas = 2), #DrySoilPeriods_SWPcrit1500kPa_MissingWater_allLayers_AnnualSum_m vs ThermalSnowfreeDryPeriods_SWPcrit3000kPa_topLayers_DrySpellsAllLayers_maxDuration_days
list(idx = 2, idy = 2, ivx = 1, ivy = 24, ixas = 2), #DrySoilPeriods_SWPcrit1500kPa_MissingWater_topLayers_PerEventPerDay_mmH2O vs ThermalSnowfreeDryPeriods_SWPcrit3000kPa_topLayers_DrySpellsAllLayers_maxDuration_days
list(idx = 2, idy = 2, ivx = 7, ivy = 11, ixas = 2), #TmaxAbovePos34degC_days vs TeeriEtAl1976_NSadj_FreezeFreeGrowingPeriod_days
list(idx = 2, idy = 2, ivx = 7, ivy = 12, ixas = 2), #TmaxAbovePos40degC_days vs TeeriEtAl1976_NSadj_FreezeFreeGrowingPeriod_days
list(idx = 2, idy = 2, ivx = 5, ivy = 26, ixas = 2), #ThermalSnowfreeWetPeriods_SWPcrit1500kPa_allLayers_AvailableWater_m vs WetSoilPeriods_SWPcrit1500kPa_NSadj_topLayers_AllLayersWet_Duration_Total_days
list(idx = 2, idy = 2, ivx = 5, ivy = 30, ixas = 2), #ThermalSnowfreeWetPeriods_SWPcrit1500kPa_allLayers_AvailableWater_mm_/_TeeriEtAl1976_NSadj_FreezeFreeGrowingPeriod_days vs WetSoilPeriods_SWPcrit1500kPa_NSadj_topLayers_AllLayersWet_Duration_Total_days
list(idx = 2, idy = 2, ivx = 8, ivy = 9, ixas = 2), #TtoAET vs TranspirationBottomToTranspirationTotal_fraction
list(idx = 2, idy = 2, ivx = 1, ivy = 9, ixas = 2) #TtoAET vs ThermalSnowfreeDryPeriods_SWPcrit3000kPa_topLayers_DrySpellsAllLayers_maxDuration_days
)
for (ils in scatter_list) {
dataX <- get(res_names_ranks[ils$idx])
dataY <- get(res_names_ranks[ils$idy])
dataX_cur <- get(res_names[ils$idx])[, currentSc, currentSc, , ]
dataY_cur <- get(res_names[ils$idy])[, currentSc, currentSc, , ]
varsX <- dimnames(dataX)[[5]]
varsY <- dimnames(dataY)[[5]]
for (itx in seq_along(type_list_mode)) for (ity in seq_along(type_list_mode)) {
ixas <- if (type_list_mode[itx] == "Mean") c("Current", "Future")[ils$ixas] else "Future"
for (ixa in ixas) {
print(paste(Sys.time(), res_names_ranks[ils$idx], res_names_ranks[ils$idy], varsX[ils$ivx], varsY[ils$ivy], type_list_mode[itx], type_list_mode[ity], ixa))
do_scatter_regional_shifted_response(datX = dataX[, , , , varsX[ils$ivx]], datY = dataY[, , , , varsY[ils$ivy]],
datX_current = dataX_cur[, , varsX[ils$ivx]], datY_current = dataY_cur[, , varsY[ils$ivy]],
axisX_type = ixa,
valX_type = type_list_mode[itx], valY_type = type_list_mode[ity],
varnameX = varsX[ils$ivx], varnameY = varsY[ils$ivy],
labelX = res_labels[[ils$idx]][ils$ivx], labelY = res_labels[[ils$idy]][ils$ivy],
val_critX = switch(EXPR = type_list_mode[itx], Mean = Inf, AbsoluteChange = 0, RelativeChange = 1),
val_critY = switch(EXPR = type_list_mode[ity], Mean = Inf, AbsoluteChange = 0, RelativeChange = 1),
probsX = switch(EXPR = type_list_mode[itx], Mean = c(0, 1), AbsoluteChange = c(0.005, 0.995), RelativeChange = c(0.05, 0.95)),
probsY = switch(EXPR = type_list_mode[ity], Mean = c(0, 1), AbsoluteChange = c(0.005, 0.995), RelativeChange = c(0.05, 0.95)),
dir_out = dir.fig_SAM)
}
}
}
}