rm(list=ls())
-
-knitr::opts_knit$set(root.dir="/Volumes/vmfileshare/ClimateData/")
-
-library(terra)
-library(sp)
-library(exactextractr)
-
-dd <- "/Volumes/vmfileshare/ClimateData/"Bias correction techniques in general require observational data to compare with climate projections in order to appropriately correct the bias.
-The HadUK grid is a 1km x 1km gridded dataset derived from meterological station observations.
-The first UKCP product for review is the UCKP convection-permitting dataset, on a 2.2km grid. Therefore, we are resmapling the 1km grid using bilenear interpolation to 2.2km grid extent.
-We have ran this seperately in both r and python. The aim of this doc is to:
Resampling script here The 2.2km grid was derived from a reprojected (to BNG) UKCP 2.2km .nc file
-In resampling it resampled the Sea as xx so replacing those vals as NA
-r1 <- paste0(dd,"TestData.R-python/Resampled_HADs_tasmax.2000.01.tif")
-r1 <- rast(r1)#Contains 31 layers for each day of Jan
-
-#In the resampling, the method used seemed to have relable all Sea values as '1.000000e+20' so relabelling them here (but to be checked as to why they've been valued like this in the resampling)
-r1[r1 > 200] = NA
-
-#check the crs
-crs(r1, proj=T)
-## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +a=6377563.396 +rf=299.324961266495 +units=m +no_defs"
-
-#Plot to check
-plot(r1$tasmax_1)
Resampling script here
-THIS UPDATED 17/02/23
-py.pros.tasmax <- list.files(paste0(dd,"Processed/HadsUKgrid/resampled_2.2km_newgrid/tasmax/day"))
-r2 <- py.pros.tasmax[grepl("200001", py.pros.tasmax)] #Same file as resampled above
-r2 <- paste0(paste0(dd, "Processed/HadsUKgrid/resampled_2.2km_newgrid/tasmax/day"),"/",r2)
-r2 <- rast(r2)
-crs(r2, proj=T) #check crs## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +a=6377563.396 +rf=299.324961266495 +units=m +no_defs"## Ok so interesting is missing a crs slot on read - I wonder why this is? This could cause future problem potentially?
-
-plot(r2$tasmax_1)
f <- paste0(dd, "Raw/HadsUKgrid/tasmax/day/")
-hads.tasmax <- list.files(f)
-
-hads.tasmax2 <- hads.tasmax[grepl("200001", hads.tasmax )] #Same file as resampled above
-og <- paste0(f, hads.tasmax2)
-
-og <- rast(og)
-crs(og, proj=T)
-
-## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +a=6377563.396 +rf=299.324961266495 +units=m +no_defs"
-
-plot(og$tasmax_1)
### 1d. UKCP example
For comparing the grids
-f <- paste0(dd,"Processed/UKCP2.2_Reproj/tasmax_bng2/01/latest/tasmax_rcp85_land-cpm_uk_2.2km_01_day_19991201-20001130.tif")
-ukcp <- rast(f)
-ukcp.r <- ukcp$`tasmax_rcp85_land-cpm_uk_2.2km_01_day_19991201-20001130_31`
-
-crs(ukcp.r, proj=T)
-
-## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
-
-#plot(ukcp.r)Just comparing by cropping to Scotland (bbox created here)
-scotland <- vect("~/Library/CloudStorage/OneDrive-TheAlanTuringInstitute/CLIM-RECAL/clim-recal/data/Scotland/Scotland.bbox.shp")Crop extents to be the same
-#Noticed the crop takes longer on r2_c - for investigation!
-
-b <- Sys.time()
-r1_c <- terra::crop(r1, scotland, snap="in")
-e <- Sys.time()
-e-b
-
-## Time difference of 0.02198005 secs
-plot(r1_c$tasmax_1)
b <- Sys.time()
-r2_c <- terra::crop(r2, scotland, snap="in")
-e <- Sys.time()
-e-b
-## Time difference of 33.57785 secs
-plot(r2_c$tasmax_1)
og_c <- terra::crop(og, scotland, snap="in")
-plot(og_c$tasmax_1)
Ok there are some differences that I can see from the plot between the two resampled files!
## Cropping to a small area to compare with the same orginal HADS file
-i <- rast()
-ext(i) <- c(200000, 210000, 700000, 710000)r1_ci <- crop(r1_c, i)
-plot(r1_ci$tasmax_1)
#Get number of cells in cropped extent
-cells <- cells(r1_ci)
-
-#get coords for all cells (for comparing above)
-r.reproj_c_xy <- sapply(cells, function(i){xyFromCell(r1_ci, i)})
-
-r.reproj_c_xy## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
-## [1,] 200935.7 203135.7 205335.7 207535.7 209735.7 200935.7 203135.7 205335.7
-## [2,] 709531.7 709531.7 709531.7 709531.7 709531.7 707331.7 707331.7 707331.7
-## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
-## [1,] 207535.7 200935.7 203135.7 205335.7 207535.7 200935.7 203135.7 205335.7
-## [2,] 707331.7 705131.7 705131.7 705131.7 705131.7 702931.7 702931.7 702931.7
-## [,17] [,18] [,19] [,20]
-## [1,] 209735.7 200935.7 203135.7 209735.7
-## [2,] 702931.7 700731.7 700731.7 700731.7ext(r1_ci)## SpatExtent : 199835.67457102, 210835.67457102, 699631.658882901, 710631.658882901 (xmin, xmax, ymin, ymax)r2_ci <- crop(r2_c, i)
-plot(r2_ci$tasmax_1)
ext(r2_ci)## SpatExtent : 199842.521629267, 210842.521629267, 699702.676089679, 710702.676089679 (xmin, xmax, ymin, ymax)og_ci <- crop(og_c, i)
-ext(og_c)## SpatExtent : 6000, 470000, 531000, 1220000 (xmin, xmax, ymin, ymax)plot(og_ci$tasmax_1)
ukcp_c <- terra::crop(ukcp.r, i)
-plot(ukcp_c$`tasmax_rcp85_land-cpm_uk_2.2km_01_day_19991201-20001130_31`)
ext(ukcp_c)## SpatExtent : 199835.67457102, 210835.67457102, 699631.658882901, 710631.658882901 (xmin, xmax, ymin, ymax)#Get number of cells in cropped extent
-cells <- cells(ukcp_c)
-
-#get coords for all cells (for comparing above)
-ukcp_c_xy <- sapply(cells, function(i){xyFromCell(ukcp_c, i)})
-
-ukcp_c_xy## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
-## [1,] 200935.7 203135.7 205335.7 207535.7 209735.7 200935.7 203135.7 205335.7
-## [2,] 709531.7 709531.7 709531.7 709531.7 709531.7 707331.7 707331.7 707331.7
-## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
-## [1,] 207535.7 209735.7 200935.7 203135.7 205335.7 207535.7 209735.7 200935.7
-## [2,] 707331.7 707331.7 705131.7 705131.7 705131.7 705131.7 705131.7 702931.7
-## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
-## [1,] 203135.7 205335.7 207535.7 209735.7 200935.7 203135.7 205335.7 207535.7
-## [2,] 702931.7 702931.7 702931.7 702931.7 700731.7 700731.7 700731.7 700731.7
-## [,25]
-## [1,] 209735.7
-## [2,] 700731.7all(ukcp_c_xy, r.reproj_c_xy)## Warning in all(ukcp_c_xy, r.reproj_c_xy): coercing argument of type 'double' to
-## logical
-
-## Warning in all(ukcp_c_xy, r.reproj_c_xy): coercing argument of type 'double' to
-## logical
-
-## [1] TRUERuth C E Bowyer 2023-06-13
-rm(list=ls())
-
-library(tidyverse)
-library(ggplot2)
-library(RColorBrewer)Script to identify the mean, 2nd highest and 2nd lowers daily tasmax per UKCP18 CPM run.
-These runs will be the focus of initial bias correction focus
-Data is tasmax runs converted to dataframe using sript ‘ConvertingAllCPMdataTOdf.R’, with files later renamed.Then daily means for historical periods and future periods were calculated using ‘calc.mean.sd.daily.R’ and summaries saved as .csv
-In retrospect the conversion to df might not have been necessary/the most resource efficient, see comment here:https://tmieno2.github.io/R-as-GIS-for-Economists/turning-a-raster-object-into-a-data-frame.html – this was tested and using terra::global to calculate the raster-wide mean was less efficient
Update 13.05.23 - Adding in infill data, mean to be calculated over the whole time period
-As of June 2023, the tasmax-as-dataframe and tasmax daily means and the df data is located in vmfileshare/Interim/tasmax_dfs/
There is an error in the naming convention - Y00_Y20 should be Y01 to reflect the infill data time period (although this does cover a breif period of 2000) - to be updated in future
-Runs <- c("01", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "15")
-
-files <- list.files("/Users/rbowyer/Library/CloudStorage/OneDrive-TheAlanTuringInstitute/tempdata/")
-
-files <- files[grepl(".csv", files)]
-fp <- paste0("/Users/rbowyer/Library/CloudStorage/OneDrive-TheAlanTuringInstitute/tempdata/", files)# Creating objects for names and filepath for each of the timer periods, for easy loading
-names <- gsub("df.avs_|.csv|df.", "", files)
-i <- c("hist", "Y00_Y20","Y21_Y40", "Y41_Y60", "Y61_Y80")
-
-namesL <- lapply(i, function(i){
- n <- names[grepl(i, names)]
- })
-
-names(namesL) <- paste0("names_",i)
-list2env(namesL, .GlobalEnv)## <environment: R_GlobalEnv>dfL <- lapply(i, function(i){
- fp <- fp[grepl(i, fp)]
- dfs <- lapply(fp, read.csv)
- n <- namesL[[paste0("names_",i)]]
- names(dfs) <- n
- return(dfs)
- })
-
-names(dfL) <- paste0("dfs_", i)
-list2env(dfL, .GlobalEnv)## <environment: R_GlobalEnv>Y <- rep(c(1981:2000), each=360)
-
-dfs_hist <- lapply(names_hist, function(i){
- df <- dfs_hist[[i]]
- names(df) <- c("day", "mean", "sd")
- df$model <- i
- df$dn <- 1:nrow(df)
- df$Y <- Y
- return(df)
-})
-
-#Create a single df in long form of Runs for the historical period
-historical_means <- dfs_hist %>% reduce(rbind)ggplot(historical_means) +
- geom_line(aes(x=dn, y=mean, group=model, colour=model)) +
-
- theme_bw() + xlab("Day (Historical 1980 - 2000)") +
- ylab("Daily mean max temp (tasmax) oC") +
-
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank(),
- legend.position = "none") +
- facet_wrap(.~ model, ncol=3)
#Create a pallete specific to the runs so when reordered maintain the same colours
-historical_means$model <- as.factor(historical_means$model)
-c <- brewer.pal(12, "Paired")
-my_colours <- setNames(c, levels(historical_means$model))historical_means %>%
- mutate(model = fct_reorder(model, mean, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean), y=mean, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(historical_means, aes(sample=mean, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
#Aggregating to year for annual average
-
-historical_means$Yf <- as.factor(historical_means$Y)
-
-historical_means_y <- historical_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(mean.annual=mean(mean, na.rm=T), sd.annual=sd(mean, na.rm = T))ggplot(historical_means_y) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (Historical 1980 - 2000)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
# Plotting with SDs in geom_ribbon to see if anything wildely different
-ggplot(historical_means_y) +
- geom_ribbon(aes(as.numeric(Yf), y=mean.annual,
- ymin = mean.annual - sd.annual,
- ymax= mean.annual + sd.annual,
- fill=model), alpha=0.4) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (Historical 1980 - 2000)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank()) + facet_wrap(.~model, ncol=3)
historical_means_y %>%
- mutate(model = fct_reorder(model, mean.annual, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean.annual), y=mean.annual, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max daily max temp oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(historical_means_y, aes(sample=mean.annual, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
historical_max_y <- historical_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(max=max(mean, na.rm=T))ggplot(historical_max_y) +
- geom_line(aes(x = as.numeric(Yf), y=max,
- color=model)) +
- theme_bw() + xlab("Year (Historical 1980 - 2000)") +
- ylab("Annual max of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
historical_max_y %>%
- mutate(model = fct_reorder(model, max, .fun='median')) %>%
- ggplot(aes(x=reorder(model, max), y=max, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max of mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
The daily max is quite different than the means - something to bear in mind but interesting to think about - eg Run 4 here has the 2nd lowest spread of max max temp, but is selected above based on means
-Daily means:
-#-1 removes the intercept to compare coefficients of all Runs
-av1 <- aov(mean ~ model - 1, historical_means)
-av1$coefficients[order(av1$coefficients)]## modelhist_Run10 modelhist_Run02 modelhist_Run05 modelhist_Run07 modelhist_Run09
-## 9.89052 11.06863 11.21424 11.22048 11.27647
-## modelhist_Run04 modelhist_Run03 modelhist_Run08 modelhist_Run01 modelhist_Run06
-## 11.29057 11.35848 11.45257 11.45414 11.86451
-## modelhist_Run11 modelhist_Run12
-## 11.99148 12.31870Annual means:
-av2 <- aov(mean.annual ~ model - 1, historical_means_y)
-av2$coefficients[order(av2$coefficients)]## modelhist_Run10 modelhist_Run02 modelhist_Run05 modelhist_Run07 modelhist_Run09
-## 9.89052 11.06863 11.21424 11.22048 11.27647
-## modelhist_Run04 modelhist_Run03 modelhist_Run08 modelhist_Run01 modelhist_Run06
-## 11.29057 11.35848 11.45257 11.45414 11.86451
-## modelhist_Run11 modelhist_Run12
-## 11.99148 12.31870Max of means:
-av3 <- aov(max ~ model - 1, historical_max_y)
-av3$coefficients[order(av3$coefficients)]## modelhist_Run10 modelhist_Run04 modelhist_Run02 modelhist_Run05 modelhist_Run03
-## 18.12329 18.81126 18.90054 19.01801 19.10454
-## modelhist_Run09 modelhist_Run01 modelhist_Run08 modelhist_Run07 modelhist_Run11
-## 19.23705 19.31541 19.44439 19.54981 19.57548
-## modelhist_Run06 modelhist_Run12
-## 19.88375 20.47650Max vals are different but based on means then selection would be Run 02 (2nd lowest), Run 04 & Run 03, and Run 11 (2nd lowest)
-Y <- rep(c(2021:2040), each=360)
-
-
-dfs_Y21_Y40 <- lapply(names_Y21_Y40, function(i){
- df <- dfs_Y21_Y40[[i]]
- names(df) <- c("day", "mean", "sd")
- df$model <- i
- df$dn <- 1:nrow(df)
- df$Y <- Y
- return(df)
-})
-
-#Create a single df in long form of Runs for the Y21_Y40 period
-Y21_Y40_means <- dfs_Y21_Y40 %>% reduce(rbind)ggplot(Y21_Y40_means) +
- geom_line(aes(x=dn, y=mean, group=model, colour=model)) +
- # Removing sd ribbon for ease of viewing
- #geom_ribbon(aes(x =dn, ymin = mean - sd, ymax= mean + sd), alpha=0.4) +
- theme_bw() + xlab("Daily (1980 - 2000)") +
- ylab("Daily mean max temp (tasmax) oC") +
- #scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank(),
- legend.position = "none") +
- facet_wrap(.~ model, ncol=3) + guides(fill = FALSE)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
-## of ggplot2 3.3.4.
-## This warning is displayed once every 8 hours.
-## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
-## generated.
#Create a pallete specific to the runs so when reordered maintain the same colours
-Y21_Y40_means$model <- as.factor(Y21_Y40_means$model)
-c <- brewer.pal(12, "Paired")
-my_colours <- setNames(c, levels(Y21_Y40_means$model))Y21_Y40_means %>%
- mutate(model = fct_reorder(model, mean, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean), y=mean, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(Y21_Y40_means, aes(sample=mean, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
#Aggregating to year for annual average
-
-Y21_Y40_means$Yf <- as.factor(Y21_Y40_means$Y)
-
-Y21_Y40_means_y <- Y21_Y40_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(mean.annual=mean(mean, na.rm=T), sd.annual=sd(mean, na.rm = T))ggplot(Y21_Y40_means_y) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (2021 - 2040)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
# Plotting with SDs in geom_ribbon to see if anything wildely different
-ggplot(Y21_Y40_means_y) +
- geom_ribbon(aes(as.numeric(Yf), y=mean.annual,
- ymin = mean.annual - sd.annual,
- ymax= mean.annual + sd.annual,
- fill=model), alpha=0.4) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (2021 - 2040)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank()) + facet_wrap(.~model, ncol=3)
Y21_Y40_means_y %>%
- mutate(model = fct_reorder(model, mean.annual, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean.annual), y=mean.annual, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max daily max temp oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(Y21_Y40_means_y, aes(sample=mean.annual, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
Y21_Y40_max_y <- Y21_Y40_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(max=max(mean, na.rm=T))ggplot(Y21_Y40_max_y) +
- geom_line(aes(x = as.numeric(Yf), y=max,
- color=model)) +
- theme_bw() + xlab("Year (2021 - 2040)") +
- ylab("Annual max of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
Y21_Y40_max_y %>%
- mutate(model = fct_reorder(model, max, .fun='median')) %>%
- ggplot(aes(x=reorder(model, max), y=max, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max of mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
Daily means:
-#-1 removes the intercept to compare coefficients of all Runs
-av1 <- aov(mean ~ model - 1, Y21_Y40_means)
-av1$coefficients[order(av1$coefficients)]## modelY21_Y40_Run10 modelY21_Y40_Run05 modelY21_Y40_Run07 modelY21_Y40_Run02
-## 10.93136 12.36223 12.64493 12.67791
-## modelY21_Y40_Run09 modelY21_Y40_Run03 modelY21_Y40_Run08 modelY21_Y40_Run01
-## 12.72584 12.85999 12.92934 13.03640
-## modelY21_Y40_Run04 modelY21_Y40_Run12 modelY21_Y40_Run06 modelY21_Y40_Run11
-## 13.07768 13.20011 13.38047 13.60076Annual means:
-av2 <- aov(mean.annual ~ model - 1, Y21_Y40_means_y)
-av2$coefficients[order(av2$coefficients)]## modelY21_Y40_Run10 modelY21_Y40_Run05 modelY21_Y40_Run07 modelY21_Y40_Run02
-## 10.93136 12.36223 12.64493 12.67791
-## modelY21_Y40_Run09 modelY21_Y40_Run03 modelY21_Y40_Run08 modelY21_Y40_Run01
-## 12.72584 12.85999 12.92934 13.03640
-## modelY21_Y40_Run04 modelY21_Y40_Run12 modelY21_Y40_Run06 modelY21_Y40_Run11
-## 13.07768 13.20011 13.38047 13.60076Max of means
-av3 <- aov(max ~ model - 1, Y21_Y40_max_y)
-av3$coefficients[order(av3$coefficients)]## modelY21_Y40_Run10 modelY21_Y40_Run02 modelY21_Y40_Run09 modelY21_Y40_Run03
-## 19.29044 20.69596 20.82538 21.05558
-## modelY21_Y40_Run05 modelY21_Y40_Run07 modelY21_Y40_Run08 modelY21_Y40_Run01
-## 21.09128 21.22942 21.33484 21.37443
-## modelY21_Y40_Run04 modelY21_Y40_Run06 modelY21_Y40_Run12 modelY21_Y40_Run11
-## 21.49363 21.98667 22.09476 22.65178Based on means then selection would be Run 02 (2nd lowest), Run 04 & Run 03, and Run 11 (2nd lowest)
-Based on this period, the seelction would be: Run 05, Run 03, Run 08, Run 06 (so definetly Run 3 but others to be discussed)
-Y <- rep(c(2061:2080), each=360)
-
-
-dfs_Y61_Y80 <- lapply(names_Y61_Y80, function(i){
- df <- dfs_Y61_Y80[[i]]
- names(df) <- c("day", "mean", "sd")
- df$model <- i
- df$dn <- 1:nrow(df)
- df$Y <- Y
- return(df)
-})
-
-#Create a single df in long form of Runs for the Y61_Y80 period
-Y61_Y80_means <- dfs_Y61_Y80 %>% reduce(rbind)ggplot(Y61_Y80_means) +
- geom_line(aes(x=dn, y=mean, group=model, colour=model)) +
- # Removing sd ribbon for ease of viewing
- #geom_ribbon(aes(x =dn, ymin = mean - sd, ymax= mean + sd), alpha=0.4) +
- theme_bw() + xlab("Day (2060 - 2080)") +
- ylab("Daily mean max temp (tasmax) oC") +
- #scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank(),
- legend.position = "none") +
- facet_wrap(.~ model, ncol=3) + guides(fill = FALSE)
#Create a pallete specific to the runs so when reordered maintain the same colours
-Y61_Y80_means$model <- as.factor(Y61_Y80_means$model)
-c <- brewer.pal(12, "Paired")
-my_colours <- setNames(c, levels(Y61_Y80_means$model))Y61_Y80_means %>%
- mutate(model = fct_reorder(model, mean, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean), y=mean, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(Y61_Y80_means, aes(sample=mean, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
#Aggregating to year for annual average
-
-Y61_Y80_means$Yf <- as.factor(Y61_Y80_means$Y)
-
-Y61_Y80_means_y <- Y61_Y80_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(mean.annual=mean(mean, na.rm=T), sd.annual=sd(mean, na.rm = T))ggplot(Y61_Y80_means_y) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (2061 - 2080)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
# Plotting with SDs in geom_ribbon to see if anything wildely different
-ggplot(Y61_Y80_means_y) +
- geom_ribbon(aes(as.numeric(Yf), y=mean.annual,
- ymin = mean.annual - sd.annual,
- ymax= mean.annual + sd.annual,
- fill=model), alpha=0.4) +
- geom_line(aes(x = as.numeric(Yf), y=mean.annual,
- color=model)) +
- theme_bw() + xlab("Year (2061 - 2080)") +
- ylab("Annual mean of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank()) + facet_wrap(.~model, ncol=3)
Y61_Y80_means_y %>%
- mutate(model = fct_reorder(model, mean.annual, .fun='median')) %>%
- ggplot(aes(x=reorder(model, mean.annual), y=mean.annual, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max daily max temp oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
ggplot(Y61_Y80_means_y, aes(sample=mean.annual, colour=factor(model))) +
- stat_qq() +
- stat_qq_line()+
- theme_bw()+
- scale_color_manual(values = my_colours) +
- facet_wrap(.~model, ncol=3)
Y61_Y80_max_y <- Y61_Y80_means %>%
- group_by(Yf, model) %>%
- dplyr::summarise(max=max(mean, na.rm=T))ggplot(Y61_Y80_max_y) +
- geom_line(aes(x = as.numeric(Yf), y=max,
- color=model)) +
- theme_bw() + xlab("Year (2061 - 2080)") +
- ylab("Annual max of mean daily max temp (tasmax) oC") +
- scale_fill_brewer(palette = "Paired", name = "") +
- scale_colour_brewer(palette = "Paired", name = "") +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
Y61_Y80_max_y %>%
- mutate(model = fct_reorder(model, max, .fun='median')) %>%
- ggplot(aes(x=reorder(model, max), y=max, fill=model)) +
- geom_boxplot() + theme_bw() +
- ylab("Annual max of mean daily max temp (tasmax) oC") + xlab("model") +
- scale_fill_manual(values = my_colours) +
- theme(axis.text.x=element_blank(),
- axis.ticks.x=element_blank())
Daily means:
-#-1 removes the intercept to compare coefficients of all Runs
-av1 <- aov(mean ~ model - 1, Y61_Y80_means)
-av1$coefficients[order(av1$coefficients)]## modelY61_Y80_Run10 modelY61_Y80_Run05 modelY61_Y80_Run01 modelY61_Y80_Run08
-## 12.70342 13.87016 14.55815 14.65973
-## modelY61_Y80_Run04 modelY61_Y80_Run09 modelY61_Y80_Run03 modelY61_Y80_Run12
-## 14.69527 14.76917 14.79545 14.87939
-## modelY61_Y80_Run07 modelY61_Y80_Run02 modelY61_Y80_Run11 modelY61_Y80_Run06
-## 14.94320 15.01577 15.11392 15.11814Annual means:
-av2 <- aov(mean.annual ~ model - 1, Y61_Y80_means_y)
-av2$coefficients[order(av2$coefficients)]## modelY61_Y80_Run10 modelY61_Y80_Run05 modelY61_Y80_Run01 modelY61_Y80_Run08
-## 12.70342 13.87016 14.55815 14.65973
-## modelY61_Y80_Run04 modelY61_Y80_Run09 modelY61_Y80_Run03 modelY61_Y80_Run12
-## 14.69527 14.76917 14.79545 14.87939
-## modelY61_Y80_Run07 modelY61_Y80_Run02 modelY61_Y80_Run11 modelY61_Y80_Run06
-## 14.94320 15.01577 15.11392 15.11814Max of means
-av3 <- aov(max ~ model - 1, Y61_Y80_max_y)
-av3$coefficients[order(av3$coefficients)]## modelY61_Y80_Run10 modelY61_Y80_Run05 modelY61_Y80_Run03 modelY61_Y80_Run04
-## 21.83290 23.32972 23.88512 23.98220
-## modelY61_Y80_Run02 modelY61_Y80_Run01 modelY61_Y80_Run08 modelY61_Y80_Run06
-## 23.98610 24.03094 24.13232 24.41824
-## modelY61_Y80_Run12 modelY61_Y80_Run09 modelY61_Y80_Run07 modelY61_Y80_Run11
-## 24.48810 24.53152 24.77651 25.09102Runs suggested by this slice are Run 05, Run 09, Run 03 and Run 11
-Run 3 and 5 suggested above
-The result per time slice suggest different runs, aside from run 5
-Update 13.05.23 - Adding in the infill data, and taking the anova result across the whole time period
-Y <- rep(c(2001:2020), each=360)
-
-dfs_Y00_Y20 <- lapply(names_Y00_Y20, function(i){
- df <- dfs_Y00_Y20[[i]]
- names(df) <- c("day", "mean", "sd")
- df$model <- i
- df$dn <- 1:nrow(df)
- df$Y <- Y
- df$Yf <- as.factor(df$Y)
- return(df)
-})
-
-
-Y <- rep(c(2041:2060), each=360)
-
-dfs_Y41_Y60 <- lapply(names_Y41_Y60, function(i){
- df <- dfs_Y41_Y60[[i]]
- names(df) <- c("day", "mean", "sd")
- df$model <- i
- df$dn <- 1:nrow(df)
- df$Y <- Y
- df$Yf <- as.factor(df$Y)
- return(df)
-})
-
-
-#Create a single df in long form as above
-Y00_Y20_means <- dfs_Y00_Y20 %>% reduce(rbind)
-Y41_Y60_means <- dfs_Y41_Y60 %>% reduce(rbind)Assessing what the combined times slices suggest via anova
-#-1 removes the intercept to compare coefficients of all Runs
-all.means <- rbind(historical_means, Y00_Y20_means, Y21_Y40_means, Y41_Y60_means, Y61_Y80_means)
-
-x <- as.character(all.means$model)
-all.means$model <- substr(x, nchar(x)-4, nchar(x))
-
-
-av1 <- aov(mean ~ model - 1, all.means)
-av1$coefficients[order(av1$coefficients)]## modelRun10 modelRun05 modelRun09 modelRun04 modelRun03 modelRun07 modelRun08
-## 11.12464 12.48165 12.79216 12.89910 12.91685 12.91894 12.95115
-## modelRun02 modelRun01 modelRun12 modelRun06 modelRun11
-## 12.95347 12.97947 13.38267 13.40644 13.61157# As above, creating annual means
-infill.L <- list(Y00_Y20_means, Y41_Y60_means)
-
-infill.L_y <- lapply(infill.L, function(x){
- means_y <- x %>%
- group_by(Yf, model) %>%
- dplyr::summarise(mean.annual=mean(mean, na.rm=T), sd.annual=sd(mean, na.rm = T))})## `summarise()` has grouped output by 'Yf'. You can override using the `.groups`
-## argument.
-## `summarise()` has grouped output by 'Yf'. You can override using the `.groups`
-## argument.all.means_y <- rbind(historical_means_y,
- infill.L_y[[1]],
- Y21_Y40_means_y,
- infill.L_y[[2]],
- Y61_Y80_means_y)
-
-x <- as.character(all.means_y$model)
-all.means_y$model <- substr(x, nchar(x)-4, nchar(x))
-
-av2 <- aov(mean.annual ~ model - 1, all.means_y)
-av2$coefficients[order(av2$coefficients)]## modelRun10 modelRun05 modelRun09 modelRun04 modelRun03 modelRun07 modelRun08
-## 11.12464 12.48165 12.79216 12.89910 12.91685 12.91894 12.95115
-## modelRun02 modelRun01 modelRun12 modelRun06 modelRun11
-## 12.95347 12.97947 13.38267 13.40644 13.61157Updated June 13th 2023 result
-Considering all together, suggests: Runs 05, Run07, Run08 and Run06
- - -clim_recal.debiasing.debias_wrapper
Wrapper for running preprocess_data.py and run_cmethods.py
| Name | +Description | +
|---|---|
| BaseRunConfig | +Manage creating command line scripts to run debiasing cli. |
+
| RunConfig | +Manage creating command line scripts to run debiasing cli. |
+
| RunConfigType | +Parameters needed for a model run. | +
clim_recal.debiasing.debias_wrapper.BaseRunConfig(self, command_dir=COMMAND_DIR_DEFAULT, run_prefix=RUN_PREFIX_DEFAULT, preprocess_data_file=PREPROCESS_FILE_NAME, run_cmethods_file=CMETHODS_FILE_NAME, data_path=DATA_PATH_DEFAULT, mod_folder=MOD_FOLDER_DEFAULT, obs_folder=OBS_FOLDER_DEFAULT, preprocess_out_folder=PREPROCESS_OUT_FOLDER_DEFAULT, cmethods_out_folder=CMETHODS_OUT_FOLDER_DEFAULT, calib_date_start=CALIB_DATE_START_DEFAULT, calib_date_end=CALIB_DATE_END_DEFAULT, valid_date_start=VALID_DATE_START_DEFAULT, valid_date_end=VALID_DATE_END_DEFAULT, processors=PROCESSORS_DEFAULT, date_format_str=CLI_DATE_FORMAT_STR, date_split_str=DATE_FORMAT_SPLIT_STR)
Manage creating command line scripts to run debiasing cli.
clim_recal.debiasing.debias_wrapper.RunConfig(self, command_dir=COMMAND_DIR_DEFAULT, run_prefix=RUN_PREFIX_DEFAULT, preprocess_data_file=PREPROCESS_FILE_NAME, run_cmethods_file=CMETHODS_FILE_NAME, data_path=DATA_PATH_DEFAULT, mod_folder=MOD_FOLDER_DEFAULT, obs_folder=OBS_FOLDER_DEFAULT, preprocess_out_folder=PREPROCESS_OUT_FOLDER_DEFAULT, cmethods_out_folder=CMETHODS_OUT_FOLDER_DEFAULT, calib_date_start=CALIB_DATE_START_DEFAULT, calib_date_end=CALIB_DATE_END_DEFAULT, valid_date_start=VALID_DATE_START_DEFAULT, valid_date_end=VALID_DATE_END_DEFAULT, processors=PROCESSORS_DEFAULT, date_format_str=CLI_DATE_FORMAT_STR, date_split_str=DATE_FORMAT_SPLIT_STR, variable=VariableOptions.default(), run=RunOptions.default(), region=RegionOptions.default(), method=MethodOptions.default())
Manage creating command line scripts to run debiasing cli.
| Name | +Description | +
|---|---|
| command_path | +Return command path relative to running tests. | +
| run_prefix_tuple | +Split self.run_prefix by ’ ’ to a tuple. |
+
| Name | +Description | +
|---|---|
| calib_dates_to_str | +Return date range as str from calib_date_start to calib_date_end. |
+
| cmethods_out_path | +Return path to save cmethods results. | +
| mod_path | +Return region estimates path. | +
| obs_path | +Return region observations path. | +
| preprocess_out_path | +Return path to save results. | +
| to_cli_preprocess_str | +Generate a command line interface str as a test example. | +
| to_cli_preprocess_tuple | +Generate a tuple of str for a command line command. |
+
| to_cli_preprocess_tuple_strs | +Generate a command line interface str tuple a test example. |
+
| to_cli_run_cmethods_str | +Generate a command line interface str as a test example. | +
| to_cli_run_cmethods_tuple | +Generate a tuple of str for a command line command. |
+
| to_cli_run_cmethods_tuple_strs | +Generate a command line interface str tuple a test example. |
+
| valid_dates_to_str | +Return date range as str from valid_date_start to valid_date_end. |
+
| yield_mod_folder | +Iterable of all Paths in self.mod_folder. |
+
| yield_obs_folder | +Iterable of all Paths in self.obs_folder. |
+
| yield_preprocess_out_folder | +Iterable of all Paths in self.preprocess_out_folder. |
+
clim_recal.debiasing.debias_wrapper.RunConfig.calib_dates_to_str(start_date, end_date, in_format_str=None, out_format_str=None, split_str=None)
Return date range as str from calib_date_start to calib_date_end.
>>> config: RunConfig = RunConfig()
+>>> config.calib_dates_to_str('20100101', '20100330')
+'20100101-20100330'
+>>> config.calib_dates_to_str(date(2010, 1, 1), '20100330')
+'20100101-20100330'
+>>> config.calib_dates_to_str(date(2010, 1, 1), '20100330', split_str="_")
+'20100101_20100330'clim_recal.debiasing.debias_wrapper.RunConfig.cmethods_out_path(region=None, run=None)
Return path to save cmethods results.
+>>> config: RunConfig = RunConfig()
+>>> config.cmethods_out_path()
+PosixPath('/.../ClimateData/Debiased/three.cities.cropped/Manchester/05')
+>>> config.cmethods_out_path(region='Glasgow', run='07')
+PosixPath('/.../ClimateData/Debiased/three.cities.cropped/Glasgow/07')clim_recal.debiasing.debias_wrapper.RunConfig.mod_path(region=None)
Return region estimates path.
+>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> config.mod_path()
+PosixPath('/.../ClimateData/Cropped/three.cities/CPM/Manchester')
+>>> config.mod_path('Glasgow')
+PosixPath('/.../ClimateData/Cropped/three.cities/CPM/Glasgow')clim_recal.debiasing.debias_wrapper.RunConfig.obs_path(region=None)
Return region observations path.
+>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> config.obs_path()
+PosixPath('/.../ClimateData/Cropped/three.cities/Hads.updated360/Manchester')
+>>> config.obs_path('Glasgow')
+PosixPath('/.../ClimateData/Cropped/three.cities/Hads.updated360/Glasgow')clim_recal.debiasing.debias_wrapper.RunConfig.preprocess_out_path(region=None, run=None, variable=None)
Return path to save results.
+>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> config.preprocess_out_path()
+PosixPath('.../ClimateData/Cropped/three.cities/Preprocessed/Manchester/05/tasmax')
+>>> config.preprocess_out_path(region='Glasgow', run='07')
+PosixPath('.../ClimateData/Cropped/three.cities/Preprocessed/Glasgow/07/tasmax')clim_recal.debiasing.debias_wrapper.RunConfig.to_cli_preprocess_str(variable=None, run=None, region=None, calib_start=None, calib_end=None, valid_start=None, valid_end=None)
Generate a command line interface str as a test example.
+>>> config: RunConfig = RunConfig()
+>>> config.to_cli_preprocess_str() == CLI_PREPROCESS_DEFAULT_COMMAND_STR_CORRECT
+True
+>>> CLI_PREPROCESS_DEFAULT_COMMAND_STR_CORRECT
+'python preprocess_data.py --mod /.../CPM/Manchester...'clim_recal.debiasing.debias_wrapper.RunConfig.to_cli_preprocess_tuple(variable=None, run=None, region=None, calib_start=None, calib_end=None, valid_start=None, valid_end=None)
Generate a tuple of str for a command line command.
This will leave Path objects uncoverted. See self.to_cli_preprocess_tuple_strs for passing to a terminal.
>>> config: RunConfig = RunConfig()
+>>> command_str_tuple: tuple[str, ...] = config.to_cli_preprocess_tuple()
+>>> assert command_str_tuple == CLI_PREPROCESS_DEFAULT_COMMAND_TUPLE_CORRECTclim_recal.debiasing.debias_wrapper.RunConfig.to_cli_preprocess_tuple_strs(variable=None, run=None, region=None, calib_start=None, calib_end=None, valid_start=None, valid_end=None)
Generate a command line interface str tuple a test example.
>>> config: RunConfig = RunConfig()
+>>> command_str_tuple: tuple[str, ...] = config.to_cli_preprocess_tuple_strs()
+>>> assert command_str_tuple == CLI_PREPROCESS_DEFAULT_COMMAND_TUPLE_STR_CORRECTclim_recal.debiasing.debias_wrapper.RunConfig.to_cli_run_cmethods_str(region=None, run=None, variable=None, method=None, input_data_path=None, cmethods_out_path=None, processors=None)
Generate a command line interface str as a test example.
+>>> config: RunConfig = RunConfig()
+>>> config.to_cli_run_cmethods_str() == CLI_CMETHODS_DEFAULT_COMMAND_STR_CORRECT
+True
+>>> CLI_CMETHODS_DEFAULT_COMMAND_STR_CORRECT
+'python run_cmethods.py...--method quantile_delta_mapping...'clim_recal.debiasing.debias_wrapper.RunConfig.to_cli_run_cmethods_tuple(region=None, run=None, variable=None, method=None, input_data_path=None, cmethods_out_path=None, processors=None)
Generate a tuple of str for a command line command.
This will leave Path objects uncoverted. See self.to_cli_run_cmethods_tuple_strs for passing to a terminal.
>>> config: RunConfig = RunConfig()
+>>> command_str_tuple: tuple[str, ...] = config.to_cli_run_cmethods_tuple()
+>>> assert command_str_tuple == CLI_CMETHODS_DEFAULT_COMMAND_TUPLE_CORRECTclim_recal.debiasing.debias_wrapper.RunConfig.to_cli_run_cmethods_tuple_strs(region=None, run=None, variable=None, method=None, input_data_path=None, cmethods_out_path=None, processors=None)
Generate a command line interface str tuple a test example.
>>> config: RunConfig = RunConfig()
+>>> command_str_tuple: tuple[str, ...] = config.to_cli_run_cmethods_tuple_strs()
+>>> assert command_str_tuple == CLI_CMEHTODS_DEFAULT_COMMAND_TUPLE_STR_CORRECTclim_recal.debiasing.debias_wrapper.RunConfig.valid_dates_to_str(start_date, end_date, in_format_str=None, out_format_str=None, split_str=None)
Return date range as str from valid_date_start to valid_date_end.
>>> config: RunConfig = RunConfig()
+>>> config.valid_dates_to_str('20100101', '20100330')
+'20100101-20100330'
+>>> config.valid_dates_to_str(date(2010, 1, 1), '20100330')
+'20100101-20100330'
+>>> config.valid_dates_to_str(date(2010, 1, 1), '20100330', split_str="_")
+'20100101_20100330'clim_recal.debiasing.debias_wrapper.RunConfig.yield_mod_folder(region=None)
Iterable of all Paths in self.mod_folder.
>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> len(tuple(config.yield_mod_folder())) == MOD_FOLDER_FILES_COUNT_CORRECT
+Trueclim_recal.debiasing.debias_wrapper.RunConfig.yield_obs_folder(region=None)
Iterable of all Paths in self.obs_folder.
>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> len(tuple(config.yield_obs_folder())) == OBS_FOLDER_FILES_COUNT_CORRECT
+Trueclim_recal.debiasing.debias_wrapper.RunConfig.yield_preprocess_out_folder(region=None, run=None, variable=None)
Iterable of all Paths in self.preprocess_out_folder.
>>> if not is_data_mounted:
+... pytest.skip('requires "vmfileshare/ClimateData" mounted')
+>>> config: RunConfig = RunConfig()
+>>> (len(tuple(config.yield_preprocess_out_folder())) ==
+... PREPROCESS_OUT_FOLDER_FILES_COUNT_CORRECT)
+Trueclim_recal.debiasing.debias_wrapper.RunConfigType()
Parameters needed for a model run.
+ + +