-
Notifications
You must be signed in to change notification settings - Fork 6
/
DownscalePT_full_CMI6_V1.0.r
247 lines (189 loc) · 8.07 KB
/
DownscalePT_full_CMI6_V1.0.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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
################################################################################
# STATISTICAL DOWNSCALING OF DAILY CLIMATE DATA USING #
# QUANTILE MAPPING TECHNIQUE FOR CMIP6 DATASETS #
################################################################################
#Author: Julio Montenegro Gambini, P.E. ASCE, M.Sc.,
#PhD fellow - Technische Universiteit Delft (TU Delft), Netherlands.
#Current version: 1.0
#©Copyright 2013 2021 Julio Montenegro.
#This script is strictly under license GPLv3
#License details: https://www.gnu.org/licenses/gpl-3.0.en.html
# Please, when using this script, cite as: "Montenegro, J. (2021). Statistical
#downscaling of daily climate data using quantile mapping technique for CMIP6
#datasets"
# Installing or loading the required packages ==================================
library(tidyverse)
library(lubridate)
library(qmap)
library(zoo)
library(latticeExtra)
library(readxl)
#Directories and folders =======================================================
# Setting working directory
setwd('C:/Users/Julio/Downloads/QMAP_CMIP6/qmap2/')
# Name of the folder which contains csv and xlsx excel files
dir_estation <- '585 M3'
# Parameter configuration ======================================================
# DAILY MINIMUM TEMPERATURE PARAMETERS
file_tmin_obs <- 'TASMIN_OBS.xlsx'
file_tmin_mod <- 'TASMIN.csv'
fecha_in_his_min <- as.Date('1981-01-01')
fecha_fin_his_min <- as.Date('2019-12-31')
fecha_in_min <- as.Date('2015-01-01')
fecha_fin_min <- as.Date('2099-12-31')
# DAILY MAXIMUM TEMPERATURE PARAMETERS
file_tmax_obs <- 'TASMAX_OBS.xlsx'
file_tmax_mod <- 'TASMAX.csv'
fecha_in_his_max <- as.Date('1981-01-01')
fecha_fin_his_max <- as.Date('2019-12-31')
fecha_in_max <- as.Date('2015-01-01')
fecha_fin_max <- as.Date('2099-12-31')
# DAILY PRECIPITATION PARAMETERS
file_pr_obs <- 'PR_OBS.xlsx'
file_pr_mod <- 'PR.csv'
fecha_in_his_pr <- as.Date('1981-01-01')
fecha_fin_his_pr <- as.Date('2016-12-31')
fecha_in_pr <- as.Date('2015-01-01')
fecha_fin_pr <- as.Date('2099-12-31')
#Quantile mapping function =====================================================
# Downscaling function for each variable
qp_function <- function(time_ini_his, time_fin_his,
time_ini, time_fin,
file_his, file_mod,
var){
time_ini <- fecha_in_min
time_fin <- fecha_fin_min
time_ini_his <- fecha_in_his_min
time_fin_his <- fecha_fin_his_min
file_his <- file_tmin_obs
file_mod <- file_tmin_mod
var = 'MIN'
j <- 2
## Reading historical baseline and GCM/RCM future time series
data_historica <- read_excel(paste0(dir_estation,'/',file_his))
data_modelada <- read.csv(paste0(dir_estation,'/',file_mod))
names_stations <- names(data_historica)[2:length(data_historica)]
df_reg <- c() # empty dataframe for filling with downscaled data
## Downscaling loop for different time series (station data)
for (j in 2:ncol(data_historica)) {
## Historical baseline configuration
data_his <- data_historica[,c(1,j)]
colnames(data_his) <- c('FECHA','STATION')
data_his <- data_his %>% mutate(FECHA= as.Date(FECHA))
## GCM/RCM data configuration
if (var %in% c('MIN','MAX')) {
correccion <- function(x, na.rm=FALSE) (x-273.15)
}
if (var %in% c('PR')) {
correccion <- function(x, na.rm=FALSE) (x*86400)
}
data_mod <- data_modelada[,c(1,j)] %>%
mutate_if(is.numeric, correccion, na.rm=FALSE)
colnames(data_mod) <- c('FECHA','STATION')
data_model <- data_mod %>% mutate(FECHA= as.Date(FECHA))
## Date filtering
### Historical
var_hist_2 <- data_his %>%
filter(FECHA >= time_ini_his & FECHA <= time_fin_his)
colnames(var_hist_2) <- c('isodate','hist')
var_hist_2 <- var_hist_2 %>%
full_join(
data.frame(isodate = seq(from=time_ini_his, to=time_fin_his, by ='day')),
by = 'isodate')
### Setting the variable of GCM/RCM data
var_model_2 <- data_model %>%
filter(FECHA >= time_ini & FECHA <= time_fin)
colnames(var_model_2) <- c('isodate','mode')
var_model_2 <- var_model_2 %>%
full_join(
data.frame(isodate = seq(from=time_ini, to=time_fin, by ='day')),
by = 'isodate')
## Creating variables to be used
OBS_hist <- var_hist_2 %>%
rename(OBS_hist = hist)
GCM_model <- var_model_2
# Filling missing historical data
if (var=='PR') {
OBS_hist <- OBS_hist %>%
mutate(OBS_hist = ifelse(is.na(OBS_hist),0.1,OBS_hist))
}
if (var=='MAX' | var=='MIN') {
OBS_hist <- OBS_hist %>%
mutate(OBS_hist = ifelse(is.na(OBS_hist),16,OBS_hist))
}
# Filling missing GCM/RCM future data
if (var=='PR') {
GCM_model <- GCM_model %>%
mutate(mode = ifelse(is.na(mode),0.1,mode))
}
if (var=='MAX' | var=='MIN') {
GCM_model <- GCM_model %>%
mutate(mode = ifelse(is.na(mode),15,mode))
}
GCM_hist <- GCM_model
## Conversion of datasets to "ts" objects
data_hist <- OBS_hist %>%
read.zoo()
data_mod <- GCM_model %>%
read.zoo()
data_wt <- GCM_model %>%
read.zoo()
#============================================================================
# QUANTILE MAPPING APPLICATION (EMPIRICAL AS DEFAULT)
seasons_by_year <- list(c("December"), c("January"), c("February"),
c("March"), c("April"), c("May"),
c("June"), c("July"), c("August"),
c("September"), c("October"), c("November"))
#According to system region, month names should be changed (e.g. spanish):
#seasons_by_year <- list(c("Diciembre"), c("Enero"), c("Febrero"),
#c("Marzo"), c("Abril"), c("Mayo"),
#c("Junio"), c("Julio"), c("Agosto"),
#c("Septiembre"), c("Octubre"), c("Noviembre"))
for(i in 1:12) {
obs_sl <- data_hist[months(time(data_hist)) %in% seasons_by_year[[i]]]
mod_sl <- data_mod[months(time(data_mod)) %in% seasons_by_year[[i]]]
#GCM/RCM data, read!: L. Gudmundsson et al. (2012)
if (sum(mod_sl, na.rm = T)==0) {
mod_sl[1]<- 0.001
}
qm_fit <- fitQmapQUANT(obs = coredata(obs_sl),
mod = coredata(mod_sl),
qstep = 0.001,
nboot = 1,
wet.day = 0, # To be changed for temperatures
type = "linear")
mod_sl_qmapped <- doQmapQUANT(coredata(mod_sl), qm_fit, type = "linear")
data_wt[ months(time(data_wt)) %in% seasons_by_year[[i]]] <- mod_sl_qmapped
}
t <- as.data.frame(data_wt)
if (j == 2) {
df_reg <- t
} else{
df_reg <- cbind(df_reg, t)
}
}
# Adding column names
colnames(df_reg) <- names_stations
# Assign the date column (daily)
df_empty<- data.frame(isodate = seq(from=time_ini, to=time_fin, by ='day'))
df_out <- df_empty %>%
cbind(df_reg)
return(df_out)
}
#Applying downscaling function and exporting results ===========================
# Dowscaling for minimum temperature
tmin_reg <- qp_function(fecha_in_min, fecha_fin_min,
fecha_in_his_min, fecha_fin_his_min,
file_tmin_obs, file_tmin_mod, 'MIN')
# Dowscaling for maximum temperature
tmax_reg <- qp_function(fecha_in_max, fecha_fin_max,
fecha_in_his_max, fecha_fin_his_max,
file_tmax_obs, file_tmax_mod, 'MAX')
# Dowscaling for precipitation
pr_reg <- qp_function(fecha_in_pr, fecha_fin_pr,
fecha_in_his_pr, fecha_fin_his_pr,
file_pr_obs, file_pr_mod, 'PR')
# Exporting downscaled data in 3 csv files =====================================
tmin_reg %>% write.csv(file = paste0(dir_estation,'/tmin_reg.csv'),row.names = F)
tmax_reg %>% write.csv(file = paste0(dir_estation,'/tmax_reg.csv'),row.names = F)
tmin_reg %>% write.csv(file = paste0(dir_estation,'/pr_reg.csv'),row.names = F)