Downscalling
2022-09-05
Downscaling of Multiple GCM (Global Circulation Model) at LIBZ
Library
Read the Observation Data (Do not run!)
df_manual_libz %>%
plot.facet()## Warning: Removed 3 row(s) containing missing values (geom_path).
df_observe_pr_libz %>%
plot.facet()
Plot AWS data (Predictand)
df_observe_pr_libz %>%
ggplot(aes(x = rain_mm.sum)) +
geom_histogram(bins = 100, col= "white")
df_observe_pr_libz %>%
ggplot(aes(x = rain_mm.sum)) +
geom_histogram(bins = 100, col= "white") +
scale_x_log10()
df_observe_pr_libz %>%
ggplot(aes(x = rain_mm.sum)) +
geom_histogram(bins = 100, col= "white") +
scale_x_sqrt()
Data Preprocessing
This is an experimental package with a function for downloading and subsetting daily Global Climate Models (GCM) from NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP-CMIP6). Readers can find more data set details in: https://www.nccs.nasa.gov/services/data-collections/land-based-products/nex-gddp-cmip6
Future implementations will consider additional functions to post-process and bias-correct this data set.
NEX-GDDP-CMIP6: https://ds.nccs.nasa.gov/thredds/catalog/AMES/NEX/GDDP-CMIP6/catalog.html
Available variables
hurs: Near-Surface Relative Humidity in [%]. huss: Near-Surface Relative Humidity in [kg/kg]. pr: Precipitation (mean of the daily precipitation rate) in [kg m-2 s-1] or [mm/d]. rlds: Surface Downwelling Longwave Radiation in [W m-2]. rsds: Surface Downwelling Shortwave Radiation in [W m-2]. sfcWind: Daily-Mean Near-Surface Wind Speed in [m s-1] tas: Daily Near-Surface Air Temperature in [K] or [ºC]. tasmax: Daily Maximum Near-Surface Air Temperature in [K] or [ºC]. tasmin: Daily Minimum Near-Surface Air Temperature in [K] or [ºC]. Available scenarios
ssp126 ssp245 ssp370 ssp585
connect and Plot Netcdf file
nc_path <- "C:/Users/DidiA/Documents/download gcm/pr/ssp245/ACCESS-CM2/pr_day_ACCESS-CM2_ssp245_r1i1p1f1_gn_2015.nc"plot geospatial data
gcm.plot(nc_path, day.band = 150)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 53
## projected point(s) not finite
## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 53
## projected point(s) not finite## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the samegcm.levelplot(nc_path, day.band = 150)
Extract Merge all ncdf file on very spesific lat and lon
# lat_range of the weather station
# Precipitation Example
#nc_path <- "C:/Users/DidiA/Documents/download gcm/pr/ssp585/ACCESS-CM2/pr_day_ACCESS-CM2_ssp585_r1i1p1f1_gn_2033.nc"
#nc_path <- "C:/Users/DidiA/Documents/download gcm/pr/ssp370/GFDL-ESM4/pr_day_GFDL-ESM4_ssp370_r1i1p1f1_gr1_2026.nc"
nc_path = "C:/Users/DidiA/Documents/download gcm/pr/ssp245/MIROC6/pr_day_MIROC6_ssp245_r1i1p1f1_gn_2046.nc"
nc <- nc_open(nc_path)
lat <- ncvar_get(nc, "lat") %>%
as.numeric()
lon <- ncvar_get(nc, "lon")
lat_rng <- 0.9265431554788971
lon_rng <- 101.19326699480717
lat_ind = which.min(abs(lat - lat_rng))
lon_ind = which.min(abs(lon - lon_rng))
pr_ts <- ncvar_get(nc, "pr", start = c(lon_ind, lat_ind, 1), count=c(1,1,-1))
#lat_ind <- which(abs(lat_rng-lat)==min(abs(lat_rng-lat)))
#lon_ind <- which(abs(lon_rng-lon)==min(abs(lon_rng-lon)))
df_pr_ts <- as.data.frame(pr_ts) %>%
mutate(rain_mm.sum = pr_ts*86400) %>%
dplyr::select(-pr_ts) %>%
mutate(date =seq.Date(as.Date("2015-01-01"), as.Date("2015-12-31"),by='day'))
ggplot(df_pr_ts) +
aes(x = date, y = rain_mm.sum) +
geom_line(size = 0.5, colour = "#112446") +
theme_bw()
ggplot(df_pr_ts) +
aes(x = rain_mm.sum) +
geom_density(adjust = 1L, fill = "#112446") +
theme_minimal()
#Coordinate Libz
lat_lon=c(0.9265431554788971, 101.19326699480717)
df_cmip6_pr_ssp245_cnrm <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp245/CNRM-CM6-1",
lat_lon=lat_lon)
df_cmip6_pr_ssp370_cnrm <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp370/CNRM-CM6-1",
lat_lon=lat_lon)
df_cmip6_pr_ssp585_cnrm <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp585/CNRM-CM6-1",
lat_lon=lat_lon)
df_cmip6_pr_historical_cnrm <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/historical/CNRM-CM6-1",
lat_lon=lat_lon)
df_cmip6_pr_ssp245_miroc6 <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp245/MIROC6",
lat_lon=lat_lon)
df_cmip6_pr_ssp370_miroc6 <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp370/MIROC6",
lat_lon=lat_lon)
df_cmip6_pr_ssp585_miroc6 <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp585/MIROC6",
lat_lon=lat_lon)
df_cmip6_pr_historical_miroc6 <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/historical/MIROC6",
lat_lon=lat_lon)
df_cmip6_pr_ssp245_gfdl <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp245/GFDL-ESM4",
lat_lon=lat_lon)
df_cmip6_pr_ssp370_gfdl <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp370/GFDL-ESM4",
lat_lon=lat_lon)
df_cmip6_pr_ssp585_gfdl <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp585/GFDL-ESM4",
lat_lon=lat_lon)
df_cmip6_pr_historical_gfdl <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/historical/GFDL-ESM4",
lat_lon=lat_lon)
df_cmip6_pr_ssp245_access <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp245/ACCESS-CM2",
lat_lon=lat_lon)
df_cmip6_pr_ssp370_access <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp370/ACCESS-CM2",
lat_lon=lat_lon)
df_cmip6_pr_ssp585_access <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/ssp585/ACCESS-CM2",
lat_lon=lat_lon)
df_cmip6_pr_historical_access <- gcm.combine.ts(path = "C:/Users/DidiA/Documents/download gcm/pr/historical/ACCESS-CM2",
lat_lon=lat_lon)combine extratcted netcdf data
# scenario (2006-2014) and model (2015-2050) are separated (Dataset option 1 )
df_cmip6_pr <- list(pr_ssp245_cnrm = df_cmip6_pr_ssp245_cnrm,
pr_ssp370_cnrm = df_cmip6_pr_ssp370_cnrm,
pr_ssp585_cnrm = df_cmip6_pr_ssp585_cnrm,
pr_historical_cnrm = df_cmip6_pr_historical_cnrm,
pr_ssp245_miroc6 = df_cmip6_pr_ssp245_miroc6,
pr_ssp370_miroc6 = df_cmip6_pr_ssp370_miroc6,
pr_ssp585_miroc6 = df_cmip6_pr_ssp585_miroc6,
pr_historical_miroc6 = df_cmip6_pr_historical_miroc6,
pr_ssp245_gfdl = df_cmip6_pr_ssp245_gfdl,
pr_ssp370_gfdl = df_cmip6_pr_ssp370_gfdl,
pr_ssp585_gfdl = df_cmip6_pr_ssp585_gfdl,
pr_historical_gfdl = df_cmip6_pr_historical_gfdl,
pr_ssp245_access = df_cmip6_pr_ssp245_access,
pr_ssp370_access = df_cmip6_pr_ssp370_access,
pr_ssp585_access = df_cmip6_pr_ssp585_access,
pr_historical_access = df_cmip6_pr_historical_access
)
df_cmip6_pr_daily <- rbindlist(df_cmip6_pr, use.names=TRUE, idcol=TRUE) %>%
dplyr::rename(id=.id) %>%
bind_rows(df_observe_pr_libz ) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
ungroup() %>%
pivot_wider(names_from = id, values_from = rain_mm.sum) %>%
filter_by_time(.start_date = "2000-01-01", .end_date = "2050-12-31") %>%
mutate(pr_ssp245_cnrm = coalesce(pr_historical_cnrm, pr_ssp245_cnrm ),
pr_ssp370_cnrm = coalesce(pr_historical_cnrm, pr_ssp370_cnrm ),
pr_ssp585_cnrm = coalesce(pr_historical_cnrm, pr_ssp585_cnrm ),
pr_ssp245_miroc6 = coalesce(pr_historical_miroc6, pr_ssp245_miroc6),
pr_ssp370_miroc6 = coalesce(pr_historical_miroc6, pr_ssp370_miroc6),
pr_ssp585_miroc6 = coalesce(pr_historical_miroc6, pr_ssp585_miroc6),
pr_ssp245_access = coalesce(pr_historical_access, pr_ssp245_access ),
pr_ssp370_access = coalesce(pr_historical_access, pr_ssp370_access ),
pr_ssp585_access = coalesce(pr_historical_access, pr_ssp585_access ),
pr_ssp245_gfdl = coalesce(pr_historical_gfdl, pr_ssp245_gfdl),
pr_ssp370_gfdl = coalesce(pr_historical_gfdl, pr_ssp370_gfdl),
pr_ssp585_gfdl = coalesce(pr_historical_gfdl, pr_ssp585_gfdl)) %>%
dplyr::select(-contains("historical")) %>%
pivot_longer(!date, names_to = "id", values_to = "rain_mm.sum") %>%
mutate(raindays = case_when(rain_mm.sum>0.2 ~ 1,
rain_mm.sum<=0.2 ~0)) %>%
group_by(id) %>%
mutate(scenario = as.factor(strsplit(id, "[_]")[[2]][2])) %>% # to use strsplit, should be grouped first
mutate(var = as.factor(strsplit(id, "[_]")[[1]][1] )) %>%
mutate(mod = as.factor(strsplit(id, "[_]")[[3]][3] )) %>%
mutate(id=as.factor(id)) %>%
ungroup()## .date_var is missing. Using: date# Aggregating
df_cmip6_pr_weekly <- df_cmip6_pr_daily %>%
group_by(id) %>%
summarise_by_time(.date_var = date, .by="week",
scenario = scenario, var=var, mod=mod,
rain_mm.sum=sum(rain_mm.sum),
raindays=sum(raindays))
df_cmip6_pr_monthly <- df_cmip6_pr_daily %>%
group_by(id) %>%
summarise_by_time(.date_var = date, .by="month",
scenario = scenario, var=var, mod=mod,
rain_mm.sum=sum(rain_mm.sum),
raindays=sum(raindays))
df_cmip6_pr_yearly <- df_cmip6_pr_daily %>%
group_by(id) %>%
summarise_by_time(.date_var = date,
scenario = scenario, var=var, mod=mod,
rain_mm.sum=sum(rain_mm.sum),
raindays=sum(raindays),
.by="year")
df_cmip6_pr_daily_spread_mm <- df_cmip6_pr_daily %>%
as.data.frame() %>%
dplyr::select(date, id, rain_mm.sum) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
#pad_by_time(date, .start_date = "2006-01-01", .end_date = "2050-12-31") %>%
#filter_by_time(date, .start_date = "2006-01-01", .end_date = "2050-12-31") %>%
ungroup() %>%
spread(id, rain_mm.sum)
list.df.pr <- list(df_cmip6_pr_daily,
df_cmip6_pr_weekly,
df_cmip6_pr_monthly,
df_cmip6_pr_yearly)
list.names <- c("Daily", "Weekly", "monthly", "Yearly")
rain_mm_days <- c("rain_mm.sum", "raindays")Time Series Plot
# Rainfall and Raindays timeseries facet plot
for(i in 1:length(list.df.pr)) {
for(j in 1:length(rain_mm_days)) {
plot.fun <- function(df, labs, mm_days){
df %>%
as.data.frame() %>%
ggplot() +
aes(x = date, y = (!! sym(mm_days)), colour = scenario) +
geom_line(size = 0.5) +
scale_color_hue(direction = 1) +
theme_bw() +
facet_wrap(vars(mod)) +
labs(title=paste0(mm_days," ", labs))
}
df <- list.df.pr[i]
lab <- list.names[i]
mm_days <- rain_mm_days [j]
print(plot.fun(df, labs=lab, mm_days=mm_days))
}}## Warning: Removed 10309 row(s) containing missing values (geom_path).
## Warning: Removed 10309 row(s) containing missing values (geom_path).
## Warning: Removed 10311 row(s) containing missing values (geom_path).
## Warning: Removed 10311 row(s) containing missing values (geom_path).
## Warning: Removed 10319 row(s) containing missing values (geom_path).
## Warning: Removed 10319 row(s) containing missing values (geom_path).
## Warning: Removed 10592 row(s) containing missing values (geom_path).
## Warning: Removed 10592 row(s) containing missing values (geom_path).
for(i in 1:length(list.df.pr)) {
for(j in 1:length(rain_mm_days)) {
plot.fun <- function(df, labs, mm_days){
df %>%
as.data.frame() %>%
ggplot() +
aes(x = date, y = (!! sym(mm_days)), colour = scenario) +
geom_line(size = 0.5) +
scale_color_hue(direction = 1) +
theme_bw() +
#facet_wrap(vars(mod)) +
labs(title=paste0(mm_days," ", labs))
}
df <- list.df.pr[i]
lab <- list.names[i]
mm_days <- rain_mm_days [j]
print(plot.fun(df, labs=lab, mm_days=mm_days))
}}## Warning: Removed 10309 row(s) containing missing values (geom_path).
## Warning: Removed 10309 row(s) containing missing values (geom_path).
## Warning: Removed 10311 row(s) containing missing values (geom_path).
## Warning: Removed 10311 row(s) containing missing values (geom_path).
## Warning: Removed 10319 row(s) containing missing values (geom_path).
## Warning: Removed 10319 row(s) containing missing values (geom_path).
## Warning: Removed 10592 row(s) containing missing values (geom_path).
## Warning: Removed 10592 row(s) containing missing values (geom_path).
Distribution Plot
for(i in 1:length(list.df.pr)) {
plot.dist <- function(df, labs){
df %>%
as.data.frame() %>%
ggplot() +
aes(x = rain_mm.sum, fill = scenario, group = scenario) +
geom_density(adjust = 1L, alpha=.5) +
scale_fill_hue(direction = 1) +
labs(title = labs) +
theme_bw() +
facet_wrap(vars(mod))
}
df <- list.df.pr[i]
lab <- list.names[i]
print(plot.dist(df, labs=lab))
}## Warning: Removed 10456 rows containing non-finite values (stat_density).
## Warning: Removed 10458 rows containing non-finite values (stat_density).
## Warning: Removed 10502 rows containing non-finite values (stat_density).
## Warning: Removed 11687 rows containing non-finite values (stat_density).
## QQ PLOT
library("car")## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:dplyr':
##
## recodedf_cmip6_pr_daily %>%
as.data.frame() %>%
dplyr::select(date, id, rain_mm.sum) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
ungroup() %>%
pivot_wider(names_from = id, values_from = rain_mm.sum) %>%
pull(pr_observed_libz) %>%
qqPlot()
## [1] 14673 909df_cmip6_pr_daily %>%
as.data.frame() %>%
dplyr::select(date, id, rain_mm.sum) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
ungroup() %>%
pivot_wider(names_from = id, values_from = rain_mm.sum) %>%
pull(pr_observed_libz) %>%
hist()
df_cmip6_pr_monthly %>%
as.data.frame() %>%
dplyr::select(date, id, rain_mm.sum) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
ungroup() %>%
pivot_wider(names_from = id, values_from = rain_mm.sum) %>%
pull(pr_observed_libz) %>%
qqPlot()
## [1] 191 167df_cmip6_pr_monthly %>%
as.data.frame() %>%
dplyr::select(date, id, rain_mm.sum) %>%
group_by(id) %>%
distinct(date, .keep_all = TRUE) %>%
ungroup() %>%
pivot_wider(names_from = id, values_from = rain_mm.sum) %>%
pull(pr_observed_libz) %>%
hist()
rnorm(500, mean = 50, sd=5) %>%
qqPlot()
## [1] 230 333hist(rnorm(500, mean = 50, sd=5))
Boxplot
ggplot(df_cmip6_pr_daily) +
aes(x = "", y = rain_mm.sum, fill = scenario) +
geom_boxplot() +
scale_fill_hue(direction = 1) +
labs(title = "daily") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold")) +
facet_wrap(vars(mod))## Warning: Removed 10456 rows containing non-finite values (stat_boxplot).
ggplot(df_cmip6_pr_monthly) +
aes(x = "", y = rain_mm.sum, fill = scenario) +
geom_boxplot() +
scale_fill_hue(direction = 1) +
labs(title = "Monthly") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold")) +
facet_wrap(vars(mod))## Warning: Removed 10502 rows containing non-finite values (stat_boxplot).
ggplot(df_cmip6_pr_yearly) +
aes(x = "", y = rain_mm.sum, fill = scenario) +
geom_boxplot() +
scale_fill_hue(direction = 1) +
labs(title = "yearly") +
theme_bw() +
theme(plot.title = element_text(size = 14L, face = "bold")) +
facet_wrap(vars(mod))## Warning: Removed 11687 rows containing non-finite values (stat_boxplot).
ggplot(df_cmip6_pr_monthly %>%
mutate(month = as.factor(month(date)))) +
aes(
x = "",
y = rain_mm.sum,
fill = month,
colour = scenario
) +
geom_boxplot() +
scale_fill_viridis_d(option = "inferno", direction = 1) +
scale_color_viridis_d(option = "inferno", direction = 1) +
theme_bw() +
facet_wrap(vars(mod))## Warning: Removed 10502 rows containing non-finite values (stat_boxplot).
ggplot(df_cmip6_pr_monthly %>%
mutate(month = as.factor(month(date)))) +
aes(x = "", y = rain_mm.sum, fill = month) +
geom_boxplot() +
scale_fill_viridis_d(option = "inferno", direction = 1) +
theme_bw() +
facet_wrap(vars(scenario))## Warning: Removed 10502 rows containing non-finite values (stat_boxplot).
Correlation Heatmap
## Warning in ggcorr(., palette = "RdBu", label = TRUE, hjust = 1, size = 3, : data
## in column(s) 'date' are not numeric and were ignored
## Warning in ggcorr(., palette = "RdBu", label = TRUE, hjust = 1, size = 3, : data
## in column(s) 'date' are not numeric and were ignored
## Warning in ggcorr(., palette = "RdBu", label = TRUE, hjust = 1, size = 3, : data
## in column(s) 'date' are not numeric and were ignored
## Warning in ggcorr(., palette = "RdBu", label = TRUE, hjust = 1, size = 3, : data
## in column(s) 'date' are not numeric and were ignored
EDA (Explaotory Data Analysis)
Density Xy Plot only Access
##PCA
##
## Attaching package: 'PCAtools'## The following objects are masked from 'package:stats':
##
## biplot, screeplot
## -- removing the lower 10% of variables based on variance
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Comparison between extreme variables
Downscaling
Conventional methods (Using MBC Library)
Data Preparation
# Works for mrs and other data -------------------------------------
oc <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2000-01-01", .end_date = "2017-12-31") %>%
replace(is.na(.), 0) %>%
dplyr::select(pr_observed_libz) %>% # vector of observed samples during the calibration period.
dplyr::rename(pr=pr_observed_libz) %>%
mutate(tas=0, dtr= 0, sfcWind=0, ps =0, huss=0) %>%
data.matrix() ## .date_var is missing. Using: datemc <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2000-01-01", .end_date = "2017-12-31") %>%
replace(is.na(.), 0) %>%
dplyr::select(pr_ssp245_cnrm) %>% # vector of observed samples during the calibration period.
dplyr::rename(pr=pr_ssp245_cnrm) %>%
mutate(tas=0, dtr= 0, sfcWind=0, ps =0, huss=0) %>%
data.matrix() ## .date_var is missing. Using: datemp <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2018-01-01", .end_date = "2050-12-31") %>%
replace(is.na(.), 0) %>%
dplyr::select(pr_ssp245_cnrm) %>% # vector of observed samples during the calibration period.
dplyr::rename(pr=pr_ssp245_cnrm) %>%
mutate(tas=0, dtr= 0, sfcWind=0, ps =0, huss=0) %>%
data.matrix() ## .date_var is missing. Using: date# Works for mrs and QDM --------------------------------------------------------
oc <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2000-01-01", .end_date = "2017-12-31") %>%
distinct(date, .keep_all = TRUE) %>%
pad_by_time() %>%
replace(is.na(.), 0) %>%
pull(pr_observed_libz) %>% # vector of observed samples during the calibration period.
as.matrix() %>%
magrittr::set_colnames(c("pr")) ## .date_var is missing. Using: date## .date_var is missing. Using: date## pad applied on the interval: daymc <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2000-01-01", .end_date = "2017-12-31") %>%
distinct(date, .keep_all = TRUE) %>%
pad_by_time() %>%
replace(is.na(.), 0) %>%
pull(pr_ssp245_cnrm) %>% # vector of model outputs during the calibration period.
as.matrix() %>%
magrittr::set_colnames(c("pr"))## .date_var is missing. Using: date## .date_var is missing. Using: date## pad applied on the interval: daymp <- df_cmip6_pr_daily_spread_mm %>%
filter_by_time(.start_date = "2018-01-01", .end_date = "2050-12-31") %>%
distinct(date, .keep_all = TRUE) %>%
pad_by_time() %>%
replace(is.na(.), 0) %>%
pull(pr_ssp245_cnrm) %>% # vector of model outputs during the projected period.
as.matrix() %>%
magrittr::set_colnames(c("pr")) ## .date_var is missing. Using: date## .date_var is missing. Using: date## pad applied on the interval: day# Data Processing--------------------------------------------
#mbcn <- MBC::MBCn(o.c=oc, m.c=mc, m.p = mp) # Multivariate bias correction (N-pdft)
#test <- mbcn[[1]]
#test <- rbindlist(test)
#mbcp <- MBC::MBCp(o.c=oc, m.c=mc, m.p = mp) # Multivariate bias correction (Pearson correlation)
#mbcr <- MBC::MBCr(o.c=oc, m.c=mc, m.p = mp) # Multivariate bias correction (Spearman rank correlation)
#r2d2 <- MBC::R2D2(o.c=oc, m.c=mc, m.p = mp) # Multivariate bias correction (R2D2)
#mrs <- MBC::MRS(o.c=oc, m.c=mc, m.p = mp) # Multivariate linear rescaling using Cholesky decomposition
qdm <- MBC::QDM(o.c=oc, m.c=mc, m.p = mp) # Univariate bias correction via quantile delta mapping## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values# Okay
# mrs1 <- mrs[1] %>% as.data.frame()
# mrs2 <- mrs[2] %>% as.data.frame()
qdm1 <- qdm[1] %>%
as.data.frame()
qdm2 <- qdm[2] %>%
as.data.frame()
df <- as.data.frame(qdm1) %>%
bind_rows(as.data.frame(qdm2)) %>%
mutate(date=seq(as.Date("2000/1/1"), as.Date("2050/12/31"), "days")) %>%
left_join((df_cmip6_pr_daily_spread_mm %>%
dplyr::select(date, pr_observed_libz)), by="date") %>%
left_join((df_cmip6_pr_daily_spread_mm %>%
dplyr::select(date, pr_ssp245_cnrm) %>%
filter_by_time(.start_date = "2018-01-01" , .end_date ="2050-12-31" )), by="date")## .date_var is missing. Using: datevis_miss(df)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
df_monthly <- df %>%
summarise_by_time(mhat.c=sum(mhat.c),
pr=sum(pr),
pr_observed_libz=sum(pr_observed_libz),
pr_ssp245_cnrm=sum(pr_ssp245_cnrm),
.by="month")## .date_var is missing. Using: datedf_longer <- df %>%
pivot_longer(!date, names_to = "id")
df_monthly_longer <- df_monthly %>%
pivot_longer(!date, names_to = "id")
ggplot(df_longer) +
aes(x = date, y = value, colour = id) +
geom_line(size = 0.5) +
scale_color_hue(direction = 1) +
theme_minimal()## Warning: Removed 35512 row(s) containing missing values (geom_path).
ggplot(df_monthly_longer) +
aes(x = date, y = value, colour = id) +
geom_line(size = 0.5) +
scale_color_hue(direction = 1) +
theme_minimal()## Warning: Removed 1167 row(s) containing missing values (geom_path).
ggpairs(df)## Warning: Removed 12053 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 12053 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 12053 rows containing missing values## Warning: Removed 6575 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6575 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 16884 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6575 rows containing missing values## Warning: Removed 12053 rows containing missing values (geom_point).## Warning: Removed 6575 rows containing missing values (geom_point).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 10309 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6575 rows containing missing values## Warning: Removed 12053 rows containing missing values (geom_point).## Warning: Removed 16884 rows containing missing values (geom_point).## Warning: Removed 10309 rows containing missing values (geom_point).## Warning: Removed 10309 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 16884 rows containing missing values## Warning: Removed 6575 rows containing missing values (geom_point).
## Removed 6575 rows containing missing values (geom_point).## Warning: Removed 16884 rows containing missing values (geom_point).## Warning: Removed 6575 rows containing non-finite values (stat_density).
ggpairs(df_monthly)## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 396 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 216 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 339 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 216 rows containing missing values## Warning: Removed 396 rows containing missing values (geom_point).## Warning: Removed 396 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 396 rows containing missing values## Warning: Removed 216 rows containing missing values (geom_point).## Warning: Removed 216 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 555 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 216 rows containing missing values## Warning: Removed 339 rows containing missing values (geom_point).## Warning: Removed 396 rows containing missing values (geom_point).## Warning: Removed 555 rows containing missing values (geom_point).## Warning: Removed 339 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 555 rows containing missing values## Warning: Removed 216 rows containing missing values (geom_point).
## Removed 216 rows containing missing values (geom_point).## Warning: Removed 555 rows containing missing values (geom_point).## Warning: Removed 216 rows containing non-finite values (stat_density).
QDM
## [1] "mhat.c" "mhat.p" "date" "pr_observed_libz"
## [5] "pr_ssp245_cnrm"
