Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update.rpipeline #109

Closed
wants to merge 35 commits into from
Closed
Changes from 1 commit
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
7dd6114
Renaming file for better accuracy
RuthBowyer Aug 15, 2023
6d8875b
Renaming file for clarity
RuthBowyer Aug 15, 2023
e9f7a85
Updating func and appling BC qmapQuant to three cities
RuthBowyer Aug 15, 2023
687eb12
Moving some files to new folder
RuthBowyer Aug 24, 2023
05ad4c9
Script updates for grouping etc
RuthBowyer Aug 24, 2023
1695c31
Correcting precip error in Hads read in
RuthBowyer Aug 25, 2023
ad2845b
Correcting precip error in Hads read in
RuthBowyer Aug 25, 2023
17cca21
Merge branch 'LCAT.data' of https://github.com/alan-turing-institute/…
RuthBowyer Aug 29, 2023
b77383f
Updating fnc for scotland
RuthBowyer Aug 29, 2023
0c828f7
Updating fnc for scotland
RuthBowyer Aug 29, 2023
fd5642b
Start the BC for LCAT
RuthBowyer Aug 31, 2023
470fb08
Including fnc for Scotland
RuthBowyer Sep 5, 2023
dff994f
Deriving seasonal means
RuthBowyer Sep 5, 2023
c15beab
Loading obs data
RuthBowyer Sep 7, 2023
89df88e
Simplifying code - this might be a conflict later though
RuthBowyer Sep 7, 2023
b195855
Adding metrics, sorting out fig generation
RuthBowyer Sep 12, 2023
4e8be76
Adding in all runs rather than one and extracting rasters for general…
RuthBowyer Sep 12, 2023
049e749
Adding loops to deal with the four runs simultaneously
RuthBowyer Sep 12, 2023
5791e97
Updating with more figs
RuthBowyer Sep 15, 2023
bd9d793
conversion of val dfs to run over all runs
RuthBowyer Sep 19, 2023
894efea
update read_crop.fn.R - this could have been pushed already?
RuthBowyer Sep 19, 2023
68d7aca
Added density plot, rmse for each run
RuthBowyer Sep 21, 2023
282a75e
Cleaning up files for LCAT
RuthBowyer Sep 22, 2023
f90ac12
Updates for R pipeline
RuthBowyer Sep 22, 2023
be16f8f
Adding pre-processing to within method (qmap in this eg)
RuthBowyer Sep 26, 2023
879e608
Sep'ing script for conversion to df
RuthBowyer Sep 26, 2023
44175ac
updating for easier R pipeline
RuthBowyer Sep 28, 2023
4ca30d4
Update to script to convert city crops to data frame
RuthBowyer Sep 28, 2023
580d36f
Add notes and fixing function
RuthBowyer Sep 29, 2023
2e6fa1d
Small updated to wet.day and corrected related issue preventing run
RuthBowyer Oct 3, 2023
8678c07
Rewriting BC assessment to be more general for input
RuthBowyer Oct 4, 2023
70c4bb0
Resaving to avoid confusion
RuthBowyer Oct 4, 2023
2b90ca0
Sorting out raster plotting
RuthBowyer Oct 8, 2023
972502d
Updating metrics so as to work with new data format
RuthBowyer Oct 9, 2023
ad8851f
tidy markdown add kable
RuthBowyer Oct 9, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Deriving seasonal means
RuthBowyer committed Sep 5, 2023

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
commit dff994f61d780486fb5ada56bd525044c3141700
99 changes: 72 additions & 27 deletions R/LCAT/Assessing.BC.data.RMD
Original file line number Diff line number Diff line change
@@ -46,7 +46,7 @@ The objects within this R list are as follows:
- 'qm1.proj.a' - bias corrected values for the validation/projection period, values fitted with linear interpolation
- 'qm1.proj.b' - bias corrected values for the validation/projection period, values fitted with tricubic interpolation

## **1. Bias Correction Assessment**
## **1. Bias Correction Assessment: trends**

### **London - tasmax = Run 08**

@@ -60,6 +60,8 @@ London <- readRDS(paste0(dd,"/Debiased/R/QuantileMapping/resultsLRun08_UKI_tasma

### **1b. Check trends**

The next set of chunks visualise the data by converting back to raster, and by looking at the trends of data across all time periods

```{r convert to df and raster}
## Load a source raster to extract the crs
@@ -70,7 +72,7 @@ rast <- rast(rp)
crs <- crs(rast)
## Convert from matix to df, transpose, create x and y cols
## Convert from matix to df, transpose, create x and y cols - when run in chunk this works fine but for some reason can throw an error when run otherwise
London.df <- lapply(London, function(x){
df <- as.data.frame(t(x))
rn <- row.names(df) #The x_y coords were saves as rownames
@@ -100,7 +102,7 @@ tm_shape(London.rasts$t.cal[[1]]) + tm_raster(title="Calibration, 1980-12-01")
tm_shape(London.rasts$qm1.hist.a[[1]]) + tm_raster(title="Calibration, bias corrected, linear 1980-12-01")
tm_shape(London.rasts$qm1.hist.b[[1]]) + tm_raster(title="Calibration, bias corrected, tricubic 1980-12-01")
```
#### *Annual trends - Calibration period*
#### *Annual trends - Calibration period - daily mean*

```{r}
@@ -123,14 +125,18 @@ London.dfg <- lapply(names(London.df), function(i){
dfx_g <- dfx %>% purrr::reduce(rbind)
})
names(London.dfg) <- names(London.df)
names(London.dfg) <- c("obs.daymeans", "raw.cal.daymeans",
"raw.proj.daymeans", "bc.a.cal.daymeans",
"bc.b.cal.daymeans", "bc.a.proj.daymeans",
"bc.b.proj.daymeans")
```

```{r}
#Add a day index to align the cal and obs
London.dfg.calp <- London.dfg[c("t.obs", "t.cal", "qm1.hist.a", "qm1.hist.b")]
London.dfg.calp <- London.dfg[c("obs.daymeans", "raw.cal.daymeans",
"bc.b.cal.daymeans", "bc.a.cal.daymeans")]
London.dfg.calp <- lapply(London.dfg.calp, function(x){
x$dayi <- 1:nrow(x)
@@ -162,38 +168,67 @@ ggplot(London.dfg.calp_mm, aes(dayi, value, group=variable, colour=variable)) +
```

#### *Annual trends - Calibration period*

#### *Annual trends - Calibration period - seasonal mean*

#Annotate season based on month index
Annotate season based on month index - the dates have different formats depending on the input data (ie hads vs cpm) so am pulling out the necessary to adjust sep

```{r}
proj.raw.df.g$season <- ifelse(grepl("-12-|-01-|-02-", proj.raw.df.g$dmy), "Winter",
ifelse(grepl("-03-|-04-|-05-", proj.raw.df.g$dmy), "Spring",
ifelse(grepl("-06-|-07-|-08-", proj.raw.df.g$dmy), "Summer", "Autumn")))
#Hads/obs df
obs.daymeans.df <- London.dfg$obs.daymeans
x <- obs.daymeans.df$day
obs.daymeans.df$season <- ifelse(grepl("1231_|0131_|0228_|0229_", x), "Winter",
ifelse(grepl("0331_|0430_|0531_", x), "Spring",
ifelse(grepl("0630_|0731_|0831_", x), "Summer", "Autumn")))
#A note here - the seasons should each have 90 days but seemingly Winter and Autumn have 89 and Spring and Summer have 91 - this is due to how the manual aligning worked out and should be updated when the hads data is re-run
proj.raw.df.g$year <- as.numeric(sub("-.*", "", proj.raw.df.g$dmy))
# Mutate to a seasonal mean
obs.seasonal.mean.df <- obs.daymeans.df %>%
group_by(season_year) %>%
mutate(mean.seasonal = mean(t.obs.mean),
sd.high.seasonal = mean.seasonal + sd(t.obs.mean),
sd.low.seasonal = mean.seasonal - sd(t.obs.mean))
obs.seasonal.mean.df <- obs.seasonal.mean.df %>%
dplyr::select(season_year:sd.low.seasonal) %>% distinct()
#Grouping variable for later vars
obs.seasonal.mean.df$model <- "obs"
```

```{r}
#lapply needs to needed
#Create a season_year var than considers the same Winter season across 2 years
## i.e. - Jan 2021 is considered as Winter 2020
proj.raw.df.g$season_year <- ifelse(proj.raw.df.g$season != "Winter"| grepl("-12-", proj.raw.df.g$dmy),
paste0(proj.raw.df.g$season, "_", proj.raw.df.g$year), paste0(proj.raw.df.g$season,"_", proj.raw.df.g$year-1))
London.dfg[c("raw.cal.daymeans", "bc.b.cal.daymeans", "bc.a.cal.daymeans")]
#lapply for remaining dfs
x <- df$day
#The CPM days are consecutive 1 - 360 by year
winter <- paste0(30, "_", 1:90, collapse="|")
spring <- paste0(30, "_", 91:180, collapse="|")
summer <- paste0(30, "_", 181:270, collapse="|")
```
#HERE - sorting out below season
```{r}
df$season <- ifelse(grepl(winter, x), "Winter",
ifelse(grepl(spring, x), "Spring",
ifelse(grepl(summer, x), "Summer", "Autumn")))
#Calculate seasonal mean and SD
seasonal.mean <- proj.raw.df.g %>%
group_by(season_year) %>% mutate(mean.seasonal = mean(mean),
sd.high.seasonal = mean.seasonal + sd(mean),
sd.low.seasonal = mean.seasonal - sd(mean))
# Mutate to a seasonal mean
obs.seasonal.mean.df <- obs.daymeans.df %>%
group_by(season_year) %>%
mutate(mean.seasonal = mean(t.obs.mean),
sd.high.seasonal = mean.seasonal + sd(t.obs.mean),
sd.low.seasonal = mean.seasonal - sd(t.obs.mean))
#Remove daily vals to avoid confusion
seasonal.mean[c("mean", "sd.high", "sd.low")] <- NULL
obs.seasonal.mean.df <- obs.seasonal.mean.df %>%
dplyr::select(season_year:sd.low.seasonal) %>% distinct()
#Remove duplicate values
seasonal.mean <- distinct(seasonal.mean, season_year, .keep_all=T) #160 seasons
obs.seasonal.mean.df$model <- "obs"
```


@@ -270,7 +305,17 @@ ggplot(dfg_sm_s) +
```

### Metrics
#### *Annual trends - seasonal max*

I think visualising the daily data is not mega helpful, but now grouping to season and calculating the seasonal maxima vals (i.e. rather than means above)

#### *Validation period - annual trends - seasonal mean*

#### *Validation period - annual trends - seasonal max*


## **2. Bias Correction Assessment: Metrics**

Add in HADs data
Add in HADs data from the cal period

### mean by cell