diff --git a/.Rbuildignore b/.Rbuildignore index a61a71e..8ff2853 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -15,3 +15,11 @@ ^pkgdown$ ^doc$ ^Meta$ +^COMPLIANCE\.yaml$ +^Compare_CHELSA\.R$ +^DD_Processing_Scripts$ +^archive$ +^data_processing$ +^logo\.svg$ +^logo_icon\.png$ +^man-roxygen$ \ No newline at end of file diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 6af27d8..1490109 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -1,8 +1,4 @@ on: - push: - branches: - - main - - devl pull_request: branches: - main diff --git a/DESCRIPTION b/DESCRIPTION index ddc57b0..78694de 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,55 +1,63 @@ Package: climr Title: Downscaling climate data in R -Version: 0.1.0.9000 -Date: 23-05-2024 +Version: 0.1.1 +Date: 19-08-2024 Authors@R: c( - person("Kiri","Daust", email = "kiri.daust@gov.bc.ca", role = c("aut", "cre")), - person("Colin", "Mahony", email = "Colin.Mahony@gov.bc.ca", role = c("aut"), + person("Kiri", "Daust", , "kiri.daust@gov.bc.ca", role = c("aut", "cre")), + person("Colin", "Mahony", , "Colin.Mahony@gov.bc.ca", role = "aut", comment = c(ORCID = "0000-0002-6111-5675")), - person("Bruno", "Tremblay", email = "bruno@boostao.ca", role = c("aut"), + person("Bruno", "Tremblay", , "bruno@boostao.ca", role = "aut", comment = c(ORCID = "0000-0002-2945-356X")), - person("Ceres", "Barros", email = "ceres.barros@nrcan-rncan.gc.ca", role = c("aut"), - comment = c(ORCID = "0000-0003-4036-977X")), - person("Francois", "Bornais", email = "francois@boostao.ca", role = c("ctb")), - person(family = "Province of British Columbia", role = c("cph", "fnd"))) -Description: `climr` is an R package that builds on the downscaling concepts - operationalized in the ClimateNA tool (climatena.ca) (Wang et al. 2016). - It provides downscaling of observational and simulated climate data using change-factor - downscaling, a simple method that adds low-spatial-resolution - climate anomalies to a high-spatial-resolution reference climatological map, with additional - elevation adjustment for temperature. Elevation-adjusted monthly values of basic climate - elements (temperature and precipitation) are then used to estimate derived variables - (e.g., degree-days) based on published equations and parameters from Wang et al. 2016. - `climr` is designed to be fast and to minimize local data storage requirements. To do so, - it uses a remote PostGIS database, and optionally caches data locally. + person("Ceres", "Barros", , "ceres.barros@nrcan-rncan.gc.ca", role = "aut", + comment = c(ORCID = "0000-0003-4036-977X")), + person("Francois", "Bornais", , "francois@boostao.ca", role = "ctb"), + person(, "Province of British Columbia", role = c("cph", "fnd")) + ) +Description: Builds on the downscaling concepts operationalized in the + 'ClimateNA' tool () (Wang et al. 2016 + ). It provides downscaling of + observational and simulated climate data using change-factor + downscaling, a simple method that adds low-spatial-resolution climate + anomalies to a high-spatial-resolution reference climatological map, + with additional elevation adjustment for temperature. + Elevation-adjusted monthly values of basic climate elements + (temperature and precipitation) are then used to estimate derived + variables (e.g., degree-days) based on published equations and + parameters from Wang et al. 2016. This package is designed to be fast + and to minimize local data storage requirements. To do so, it uses a + remote 'PostGIS' database, and optionally caches data locally. License: Apache License (== 2) -Encoding: UTF-8 -Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +Depends: + R (>= 4.0) Imports: data.table, DBI, + magrittr, methods, pool, RPostgres, - terra, + scales, sf, stinepack, + terra, uuid, - scales, - magrittr + dplyr, + abind, + dbplyr Suggests: ggplot2, knitr, parallel, plotly, - rmarkdown, remotes, + rmarkdown, testthat (>= 3.0.0), utils, withr -Depends: - R (>= 4.0) +VignetteBuilder: + knitr Config/testthat/edition: 3 +Encoding: UTF-8 LazyData: true -VignetteBuilder: knitr +Roxygen: list(markdown = TRUE) +RoxygenNote: 7.3.1 diff --git a/NAMESPACE b/NAMESPACE index 6abd306..c896314 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -20,13 +20,15 @@ export(list_gcms) export(list_obs_periods) export(list_obs_years) export(list_refmaps) -export(list_run) +export(list_runs_historic) +export(list_runs_ssp) export(list_ssps) export(list_vars) export(pgGetTerra) export(plot_bivariate) export(plot_timeSeries) export(plot_timeSeries_input) +import(abind) import(data.table) import(uuid) importFrom(DBI,dbGetQuery) @@ -48,13 +50,19 @@ importFrom(data.table,setDTthreads) importFrom(data.table,setcolorder) importFrom(data.table,setkey) importFrom(data.table,setorder) +importFrom(dplyr,collect) +importFrom(dplyr,mutate) +importFrom(dplyr,sql) +importFrom(dplyr,tbl) importFrom(grDevices,hcl.colors) importFrom(grDevices,palette) importFrom(graphics,axis) +importFrom(graphics,box) importFrom(graphics,legend) importFrom(graphics,lines) importFrom(graphics,par) importFrom(graphics,points) +importFrom(graphics,polygon) importFrom(graphics,text) importFrom(graphics,title) importFrom(magrittr,"%>%") @@ -66,6 +74,7 @@ importFrom(sf,st_as_sf) importFrom(sf,st_join) importFrom(stats,as.formula) importFrom(stats,complete.cases) +importFrom(stats,smooth.spline) importFrom(stinepack,stinterp) importFrom(terra,"crs<-") importFrom(terra,as.list) @@ -92,6 +101,7 @@ importFrom(terra,xres) importFrom(terra,yres) importFrom(tools,R_user_dir) importFrom(utils,askYesNo) +importFrom(utils,data) importFrom(utils,head) importFrom(utils,tail) importFrom(uuid,UUIDgenerate) diff --git a/NEWS.md b/NEWS.md index 6f497ad..251eefb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,96 +1,136 @@ +# `climr` 0.1.1 + +## Bug Fixes + +- Modified the Hargreaves method for calculation of monthly solar radiation to allow for calculation of Eref and CMD above the arctic circle in a way that matches ClimateNA results. +- Fixed an edge case in the tiling that resulted in NA values in southern Mexico if full North American extent was queried. +- Fixed bug in caching where cache would fail due to incorrect folder name in certain cases. +- Fixed a bug in `plot_downscale()` that affected colors when `pal = "gcms"` and labels when `endlabel = "gcms"` +- Fixed a bug where climr didn't download the correct run names for historic GCMs +- Modified the equations for frost-free period so that they are confined to a range of 0-365. + +## Enhancements + +- Redesigned the database to substantially improve speed of downloading gcm ssp timeseries, especially for a small number of points. +- Added `run_nm` parameter to `downscale` and `input_gcm_*` functions so users can specify desired run(s). Also added `list_runs_ssp()` and `list_runs_historic()` + +## Known issues + +- There are discrepancies between climr and ClimateNA in values of frost-free period (eFFP, bFFP, and FFP) and frost-free days (NFFD). For FFP see issue [#297](https://github.com/bcgov/climr/issues/297); For NFFD see issue [#295](https://github.com/bcgov/climr/issues/295). We will attempt to resolve these differences in a future version. + # `climr` 0.1.0 + ## Implementation of naming conventions -* We overhauled the naming of functions, parameters, and options to make them more intuitive and internally consistent. You will need to revise the climr code in your workflows to accommodate these changes. A table of correspondence specifying the changes is located at ./data-raw/namingChanges.csv and is viewable by calling `data(name_changes)`. -* We changed the climate variable naming convention so that the climate element and the time of year are always separated by an underscore. e.g., Tmin01 becomes Tmin_01; DD_0_01 becomes DDsub0_01. The variables table called by `data(variables)` now has a field "Code_ClimateNA" with the variable codes used by ClimateBC/NA to allow users to crosswalk the two conventions. + +- We overhauled the naming of functions, parameters, and options to make them more intuitive and internally consistent. You will need to revise the climr code in your workflows to accommodate these changes. A table of correspondence specifying the changes is located at ./data-raw/namingChanges.csv and is viewable by calling `data(name_changes)`. +- We changed the climate variable naming convention so that the climate element and the time of year are always separated by an underscore. e.g., Tmin01 becomes Tmin_01; DD_0_01 becomes DDsub0_01. The variables table called by `data(variables)` now has a field "Code_ClimateNA" with the variable codes used by ClimateBC/NA to allow users to crosswalk the two conventions. ## User Actions Required -* To implement this new version, users must clear their cache of climate data by running the following line of code: cache_clear(). We will do our best to avoid the need for wholesale cache-clearing in the future. + +- To implement this new version, users must clear their cache of climate data by running the following line of code: cache_clear(). We will do our best to avoid the need for wholesale cache-clearing in the future. ## Bug Fixes -* Fixed an error in the calculation of Hogg's climatic moisture index (CMI). This error was inherited from an unreported error in Equation 3 of Hogg (1997). -* Fixed an error in the calculation of precipitation as snow. The PAS parameters in ClimateNA and `climr` differ from those originally published in Wang et al. (2016). -* fixed an out-of-bounds error that affected user queries close to the coastline. -* Added missing future 2015-2100 time series for the GFDL-ESM4 climate model. + +- Fixed an error in the calculation of Hogg's climatic moisture index (CMI). This error was inherited from an unreported error in Equation 3 of Hogg (1997).\ +- Fixed an error in the calculation of precipitation as snow. The PAS parameters in ClimateNA and `climr` differ from those originally published in Wang et al. (2016). +- fixed an out-of-bounds error that affected user queries close to the coastline. +- Added missing future 2015-2100 time series for the GFDL-ESM4 climate model. ## Enhancements -* Added `plot_timeSeries()` and `plot_timeSeries_input()` functions to generate plots of 20th and 21st century climate change for user-selected locations and climate variables. -* Added a 1901-2022 observational time series of for the combined Climatic Research Unit TS dataset (for Temperature) and Global Precipitation Climatology Centre dataset (for precipitation). -* Extended the ClimateNA observational time series to 1901-2023. -* Added a vignette (article) providing guidance for climate model ensemble selection and emissions scenario selection. -* Added a vignette (article) on the methods used to select the 13 global climate models provided by climr, and the 8-model ensemble recommended for most purposes. + +- Added `plot_timeSeries()` and `plot_timeSeries_input()` functions to generate plots of 20th and 21st century climate change for user-selected locations and climate variables. +- Added a 1901-2022 observational time series of for the combined Climatic Research Unit TS dataset (for Temperature) and Global Precipitation Climatology Centre dataset (for precipitation). +- Extended the ClimateNA observational time series to 1901-2023. +- Added a vignette (article) providing guidance for climate model ensemble selection and emissions scenario selection. +- Added a vignette (article) on the methods used to select the 13 global climate models provided by climr, and the 8-model ensemble recommended for most purposes. ## Known issues -* Downloads of time series take a long time. We are looking into ways to speed this up, but until then we recommend users dedicate some time prior to analysis to cache their time series of interest for their areas of interest in a batch. Once the time series are cached, they don't need to be downloaded again. -* Related to the issue of time series download speed, the `plot_timeSeries_input()` function can take >1hr to run for the first time it is called for a location. + +- Downloads of time series take a long time. We are looking into ways to speed this up, but until then we recommend users dedicate some time prior to analysis to cache their time series of interest for their areas of interest in a batch. Once the time series are cached, they don't need to be downloaded again. +- Related to the issue of time series download speed, the `plot_timeSeries_input()` function can take \>1hr to run for the first time it is called for a location. # `climr` 0.0.4 + ## Bug Fixes -* Updated future timeseries data to include full 2015-2100 period and added missing models (built some cool bash scrips using `parallel` to quickly to `raster2pgsql` conversion for large numbers of files). -* Updated historic modelled timeseries to extend to December 31, 2014. -* Restructured naming scheme for timeseries data, updated internal `dbnames` table, and updated `postgresql` functions to allow hyphens in table names. -* Reprocessed future GCM periods to include all of North America. -* Fixed bug in `plot_bivariate()` for focal periods after 2001-2020. -* Fixed caching issue where it would fail for very larger numbers of layers (>65000) by saving as .gri binary files. + +- Updated future timeseries data to include full 2015-2100 period and added missing models (built some cool bash scrips using `parallel` to quickly to `raster2pgsql` conversion for large numbers of files). +- Updated historic modelled timeseries to extend to December 31, 2014. +- Restructured naming scheme for timeseries data, updated internal `dbnames` table, and updated `postgresql` functions to allow hyphens in table names. +- Reprocessed future GCM periods to include all of North America. +- Fixed bug in `plot_bivariate()` for focal periods after 2001-2020. +- Fixed caching issue where it would fail for very larger numbers of layers (\>65000) by saving as .gri binary files. ## Enhancements -* Added checks for bounding box projection. + +- Added checks for bounding box projection. # `climr` 0.0.3 + ## Enhancements -* new tests comparing to reference outputs -* code further streamlined -* new messages warn user about meaningless `downscale`/`climr_downscale` argument combinations -* argument options in `climr_downscale(..., which_normal)` now match the options of `normal_input(..., normal)` -* add `plot_bivariate()` function to generate plots showing climate model ensemble variation in recent and future climate change. -* new functions `list_historic_ts` and `list_gcm_ts` to get available years for historic/future time series + +- new tests comparing to reference outputs +- code further streamlined +- new messages warn user about meaningless `downscale`/`climr_downscale` argument combinations +- argument options in `climr_downscale(..., which_normal)` now match the options of `normal_input(..., normal)` +- add `plot_bivariate()` function to generate plots showing climate model ensemble variation in recent and future climate change. +- new functions `list_historic_ts` and `list_gcm_ts` to get available years for historic/future time series ## Behaviour changes -* `xyz` (argument to `climr_downscale` and `downscale`) and `in_xyz` (argument to `get_bb`), must now be a 4 column `data.table` (or coercible class) with `lon`, `lat`, `elev` and `id` columns. All other columns are ignored and NOT returned. Column order no longer matters. + +- `xyz` (argument to `climr_downscale` and `downscale`) and `in_xyz` (argument to `get_bb`), must now be a 4 column `data.table` (or coercible class) with `lon`, `lat`, `elev` and `id` columns. All other columns are ignored and NOT returned. Column order no longer matters. ## Bugfixes -* cache fixes -* fixing geographical checks to get highest resolution beyond BC, Canada -* fixing `historic_input_ts` to get only queried years -* `get_bb` follows column names + +- cache fixes +- fixing geographical checks to get highest resolution beyond BC, Canada +- fixing `historic_input_ts` to get only queried years +- `get_bb` follows column names ## Other -* `climr_downscale` now accepts `...` to pass arguments to `downscale` + +- `climr_downscale` now accepts `...` to pass arguments to `downscale` # `climr` 0.0.2 ## Bugfixes -* fixed temperature values of composite anomalies -* name of composite anomalies changed to "normal_composite" in `normal_input(..., normal)`. + +- fixed temperature values of composite anomalies +- name of composite anomalies changed to "normal_composite" in `normal_input(..., normal)`. ## Documentation -* added vignettes -* `pkgdown` website for `climr` is live + +- added vignettes +- `pkgdown` website for `climr` is live # `climr` 0.0.1 ## Enhancements -* continuous testing implemented via GitHub Actions -* code was cleaned up following `tidyverse` syntax recommendations -* internal function definitions are now avoided -* improvements to function documentation -* removal of deprecated functions -* increased code coverage -* some code streamlining -* added new composite climatologies of Western Canada and Western US. + +- continuous testing implemented via GitHub Actions +- code was cleaned up following `tidyverse` syntax recommendations +- internal function definitions are now avoided +- improvements to function documentation +- removal of deprecated functions +- increased code coverage +- some code streamlining +- added new composite climatologies of Western Canada and Western US. ## Bugfixes -* fixed examples -* fixed incomplete changing of package name (`climRpnw` to `climr`) -* added missing pkg imports -* fixed caching problem where cached objects were not being retrieved when the PostGIS server was unavailable -* fixed model names in PostGIS server, which fixed bugs in `gcm_hist_input` and `gcm_ts_input`. + +- fixed examples +- fixed incomplete changing of package name (`climRpnw` to `climr`) +- added missing pkg imports +- fixed caching problem where cached objects were not being retrieved when the PostGIS server was unavailable +- fixed model names in PostGIS server, which fixed bugs in `gcm_hist_input` and `gcm_ts_input`. ## Dependency changes -* `methods` removed from Imports -* `sf` added to Imports + +- `methods` removed from Imports +- `sf` added to Imports ## Other -* old data base access functions removed + +- old data base access functions removed # `climr` 0.0.0.9990 diff --git a/R/append_clim_vars.R b/R/append_clim_vars.R index b55e880..dde278c 100644 --- a/R/append_clim_vars.R +++ b/R/append_clim_vars.R @@ -150,43 +150,43 @@ append_clim_vars <- function(dt, vars) { set(dt, j = "PET_12", value = calc_PET(v("Tave_12"), v("Tmin_12"), v("Tmax_12"), v("elev"))) }, "CMI" = function() { - set(dt, j = "CMI", value = (v("PPT") - (v("PET_01") + v("PET_02") + v("PET_03") + v("PET_04") + v("PET_05") + v("PET_06") + v("PET_07") + v("PET_08") + v("PET_09") + v("PET_10") + v("PET_11") + v("PET_12")))/ 10) + set(dt, j = "CMI", value = (v("PPT") - (v("PET_01") + v("PET_02") + v("PET_03") + v("PET_04") + v("PET_05") + v("PET_06") + v("PET_07") + v("PET_08") + v("PET_09") + v("PET_10") + v("PET_11") + v("PET_12"))) / 10) }, "CMI_01" = function() { - set(dt, j = "CMI_01", value = (v("PPT_01") - v("PET_01"))/ 10) + set(dt, j = "CMI_01", value = (v("PPT_01") - v("PET_01")) / 10) }, "CMI_02" = function() { - set(dt, j = "CMI_02", value = (v("PPT_02") - v("PET_02"))/ 10) + set(dt, j = "CMI_02", value = (v("PPT_02") - v("PET_02")) / 10) }, "CMI_03" = function() { - set(dt, j = "CMI_03", value = (v("PPT_03") - v("PET_03"))/ 10) + set(dt, j = "CMI_03", value = (v("PPT_03") - v("PET_03")) / 10) }, "CMI_04" = function() { - set(dt, j = "CMI_04", value = (v("PPT_04") - v("PET_04"))/ 10) + set(dt, j = "CMI_04", value = (v("PPT_04") - v("PET_04")) / 10) }, "CMI_05" = function() { - set(dt, j = "CMI_05", value = (v("PPT_05") - v("PET_05"))/ 10) + set(dt, j = "CMI_05", value = (v("PPT_05") - v("PET_05")) / 10) }, "CMI_06" = function() { - set(dt, j = "CMI_06", value = (v("PPT_06") - v("PET_06"))/ 10) + set(dt, j = "CMI_06", value = (v("PPT_06") - v("PET_06")) / 10) }, "CMI_07" = function() { - set(dt, j = "CMI_07", value = (v("PPT_07") - v("PET_07"))/ 10) + set(dt, j = "CMI_07", value = (v("PPT_07") - v("PET_07")) / 10) }, "CMI_08" = function() { - set(dt, j = "CMI_08", value = (v("PPT_08") - v("PET_08"))/ 10) + set(dt, j = "CMI_08", value = (v("PPT_08") - v("PET_08")) / 10) }, "CMI_09" = function() { - set(dt, j = "CMI_09", value = (v("PPT_09") - v("PET_09"))/ 10) + set(dt, j = "CMI_09", value = (v("PPT_09") - v("PET_09")) / 10) }, "CMI_10" = function() { - set(dt, j = "CMI_10", value = (v("PPT_10") - v("PET_10"))/ 10) + set(dt, j = "CMI_10", value = (v("PPT_10") - v("PET_10")) / 10) }, "CMI_11" = function() { - set(dt, j = "CMI_11", value = (v("PPT_11") - v("PET_11"))/ 10) + set(dt, j = "CMI_11", value = (v("PPT_11") - v("PET_11")) / 10) }, "CMI_12" = function() { - set(dt, j = "CMI_12", value = (v("PPT_12") - v("PET_12"))/ 10) + set(dt, j = "CMI_12", value = (v("PPT_12") - v("PET_12")) / 10) }, "NFFD_01" = function() { set(dt, j = "NFFD_01", value = calc_NFFD(1, v("Tmin_01"))) @@ -711,11 +711,11 @@ append_clim_vars <- function(dt, vars) { } # Remove unwanted variables - j_out <- names(dt)[!names(dt) %in% c("id","DATASET", "GCM", "SSP", "RUN", "PERIOD", vars)] + j_out <- names(dt)[!names(dt) %in% c("id", "DATASET", "GCM", "SSP", "RUN", "PERIOD", vars)] if (length(j_out)) { set(dt, j = j_out, value = NULL) } # Reorder to match vars - setcolorder(dt, c(names(dt)[names(dt) %in% c("id","DATASET", "GCM", "SSP", "RUN", "PERIOD")], vars)) + setcolorder(dt, c(names(dt)[names(dt) %in% c("id", "DATASET", "GCM", "SSP", "RUN", "PERIOD")], vars)) } diff --git a/R/calc_Eref_CMD.R b/R/calc_Eref_CMD.R index 408af26..2662261 100644 --- a/R/calc_Eref_CMD.R +++ b/R/calc_Eref_CMD.R @@ -62,6 +62,7 @@ calc_S0_I <- function(d, tm, latitude) { XLR <- latitude / 57.2958 Z <- -tan(XLR) * tan(DEC) OM <- -atan(Z / sqrt(-Z * Z + 1)) + pi / 2 + OM[!is.finite(OM)] <- if(d %in% 91:274) 3.1 else 0 # NB this is a modification of Hargreaves program to provide finite values above the arctic circle. # CALCULATE THE DAILY EXTRATERRESTRIAL RADIATION IN LANGLEYS/DAY DL <- OM / 0.1309 RAL <- 120 * (DL * sin(XLR) * sin(DEC) + 7.639 * cos(XLR) * cos(DEC) * sin(OM)) / ES diff --git a/R/calc_FFP.R b/R/calc_FFP.R index d633f6f..4000690 100644 --- a/R/calc_FFP.R +++ b/R/calc_FFP.R @@ -22,10 +22,12 @@ calc_bFFP <- function(td, NFFD, t_min_list) { tmin4 <- t_min_list[[4]] tmin6 <- t_min_list[[6]] - 352.1358994 + -0.021715653 * tmin4^2 + -3.542187618 * tmin6 + + res <- 352.1358994 + -0.021715653 * tmin4^2 + -3.542187618 * tmin6 + 0.020359471 * tmin6^2 - 4.897998097 * td + 0.033521327 * td^2 - 2.164862277 * NFFD + 0.006767633 * NFFD^2 - 0.00000929 * NFFD^3 + 0.043516586 * (td * NFFD) - 0.00000253 * (td * NFFD)^2 + res[res < 0] <- 0 + return(res) } #' Calculate eFFP @@ -53,9 +55,11 @@ calc_eFFP <- function(NFFD, t_min_list) { tmin10 <- t_min_list[[10]] tmin11 <- t_min_list[[11]] - 243.7752209 + 4.134210825 * tmin9 - 0.162876448 * tmin9^2 + + res <- 243.7752209 + 4.134210825 * tmin9 - 0.162876448 * tmin9^2 + 1.248649021 * tmin10 + 0.145073612 * tmin10^2 + 0.004319892 * tmin10 + -0.005753127 * tmin10^2 - 0.06296471 * NFFD + 0.000399177 * NFFD^2 + res[res > 365] <- 365 + return(res) } #' Calculate Frost-Free Period @@ -71,5 +75,8 @@ calc_eFFP <- function(NFFD, t_min_list) { #' } #' @noRd calc_FFP <- function(bFFP, eFFP) { - eFFP - bFFP + res <- eFFP - bFFP + res[res < 0] <- 0 + res[res > 365] <- 365 + return(res) } diff --git a/R/climr-package.R b/R/climr-package.R index 9c28599..0512e4f 100644 --- a/R/climr-package.R +++ b/R/climr-package.R @@ -5,6 +5,12 @@ ## usethis namespace: end NULL +.onAttach <- function(libname, pkgname) { + packageStartupMessage( + "Welcome to climr 0.1.1!" + ) +} + # On load, instantiate either as new or from cache #' @importFrom data.table fwrite #' @importFrom RPostgres dbGetQuery @@ -13,7 +19,6 @@ NULL .onLoad <- function(libname, pkgname) { rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") - packageStartupMessage("climr version 0.1.0 includes an overhaul of the naming conventions for variables, functions, parameters, and options. Call `data(name_changes)` for a table specifying the correspondence between old and new names. If you have used a previous version of climr you MUST call `cache_clear()` before using this one.") dbCon <- data_connect() on.exit(try(pool::poolClose(dbCon)), add = TRUE) @@ -25,8 +30,8 @@ NULL gcm_ts_runs <- dbGetQuery(dbCon, "select distinct mod, scenario, run from esm_layers_ts order by mod, scenario, run;") gcm_hist_runs <- dbGetQuery(dbCon, "select distinct mod, run from esm_layers_hist order by mod, run;") fwrite(gcm_period_runs, file.path(rInfoPath, "gcm_periods.csv")) - fwrite(gcm_period_runs, file.path(rInfoPath, "gcm_ts.csv")) - fwrite(gcm_period_runs, file.path(rInfoPath, "gcm_hist.csv")) + fwrite(gcm_ts_runs, file.path(rInfoPath, "gcm_ts.csv")) + fwrite(gcm_hist_runs, file.path(rInfoPath, "gcm_hist.csv")) } } diff --git a/R/data-lists.R b/R/data-lists.R index e0f9304..b687e3e 100644 --- a/R/data-lists.R +++ b/R/data-lists.R @@ -1,9 +1,9 @@ -#' List available runs, global circulation models, periods and climate scenarios +#' List available runs, global climate models (GCMs), time periods and scenarios (SSPs) #' #' @return a character vector. #' #' @description -#' `list_gcms` lists available global circulation models. +#' `list_gcms` lists available global climate models. #' #' @rdname data-option-lists #' @export @@ -17,7 +17,7 @@ list_gcms <- function() { } #' @description -#' `list_ssps` lists available shared socioeconomic pathways. +#' `list_ssps` lists available greenhouse gas concentration scenarios (SSP-RCPs). #' #' @rdname data-option-lists #' @export @@ -27,7 +27,7 @@ list_ssps <- function() { } #' @description -#' `list_gcm_periods` lists available periods. +#' `list_gcm_periods` lists available 20-year normal periods for GCM simulations. #' #' @rdname data-option-lists #' @export @@ -36,17 +36,34 @@ list_gcm_periods <- function() { c("2001_2020", "2021_2040", "2041_2060", "2061_2080", "2081_2100") } + #' @description -#' `list_run` lists available runs for a given GCM. -#' -#' @template dbCon -#' @param gcms Character vector to specify requested GCMs -#' @importFrom RPostgres dbGetQuery -#' +#' lists available runs for a given GCM/ssp. +#' @param gcm Name of GCM. Must be one of the elements in list_gcms(). +#' @param ssp Name of scenario Must be one of the elements in list_ssps(). +#' @importFrom data.table fread +#' @importFrom tools R_user_dir +#' +#' @rdname data-option-lists +#' @export +list_runs_ssp <- function(gcm, ssp){ + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") + runs <- fread(file.path(rInfoPath, "gcm_periods.csv")) + runs[mod == gcm & scenario == ssp, run] +} + +#' @description +#' lists available runs from the historical simulation (1851-2014) for a specified GCM. +#' @param gcm Name of GCM +#' @importFrom data.table fread +#' @importFrom tools R_user_dir +#' #' @rdname data-option-lists #' @export -list_run <- function(dbCon, gcms) { - sort(dbGetQuery(dbCon, paste0("SELECT DISTINCT run FROM esm_layers_period WHERE mod IN ('", paste(gcms, collapse = "','", "')")))[, 1]) +list_runs_historic <- function(gcm){ + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") + runs <- fread(file.path(rInfoPath, "gcm_hist.csv")) + runs[mod == gcm, run] } #' @description @@ -68,7 +85,7 @@ list_refmaps <- function() { #' @description -#' `list_obs_periods` lists available observational periods +#' `list_obs_periods` lists available normal periods for observational climate data #' #' @rdname data-option-lists #' @export @@ -77,7 +94,7 @@ list_obs_periods <- function() { } #' @description -#' `list_vars` lists climate variables +#' `list_vars` lists available climate variables #' #' @param set character. One of All, Monthly, Seasonal, Annual, or any combination thereof. Defaults to "All". #' @param only_extra logical. Should Tmin, Tmax and PPT be excluded? Defaults to FALSE. @@ -102,7 +119,7 @@ list_vars <- function(set = c("All", "Monthly", "Seasonal", "Annual"), only_extr #' @description -#' `list_obs_years` lists available years for obs time series +#' `list_obs_years` lists available years for time series of observational climate data #' #' @rdname data-option-lists #' @export @@ -111,7 +128,7 @@ list_obs_years <- function() { } #' @description -#' `list_gcm_ssp_years` lists available years for future projections' time series +#' `list_gcm_ssp_years` lists available years for time series of global climate model future simulations #' #' @rdname data-option-lists #' @export @@ -120,7 +137,7 @@ list_gcm_ssp_years <- function() { } #' @description -#' `list_gcm_hist_years` lists available years for obs projections' time series +#' `list_gcm_hist_years` lists available years for time series of global climate model historical simulations #' #' @rdname data-option-lists #' @export diff --git a/R/dbGetRaster.R b/R/dbGetRaster.R index 26b7863..37be005 100644 --- a/R/dbGetRaster.R +++ b/R/dbGetRaster.R @@ -28,7 +28,7 @@ pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73, # }else{ # name1 <- name # } - + name1 <- name nameque <- paste(name1, collapse = ".") namechar <- gsub("'", "''", paste(gsub('^"|"$', "", name1), collapse = ".")) @@ -96,11 +96,11 @@ pgGetTerra <- function(conn, name, tile, rast = "rast", bands = 37:73, boundary_ls <- list() if (length(x_seq) < 2 | length(y_seq) < 2) { - boundary_ls[["11"]] <- boundary + boundary_ls[["1_1"]] <- boundary } else { for (i in 1:(length(x_seq) - 1)) { for (j in 1:(length(y_seq) - 1)) { - boundary_ls[[paste0(i, j)]] <- c(x_seq[i + 1], x_seq[i], y_seq[j + 1], y_seq[j]) + boundary_ls[[paste0(i,"_", j)]] <- c(x_seq[i + 1], x_seq[i], y_seq[j + 1], y_seq[j]) } } } diff --git a/R/downscale.R b/R/downscale.R index 54c600f..db8146c 100644 --- a/R/downscale.R +++ b/R/downscale.R @@ -37,8 +37,8 @@ #' series of observational climate data. Default `NULL`. See [`list_obs_years()`] #' for available years. #' @param obs_ts_dataset character. The dataset to use for observational time series data. Options -#' are `"climatena"` for the ClimateNA gridded time series or `"cru.gpcc"` for the combined Climatic -#' Research Unit TS dataset (for temperature) and Global Precipitation Climatology Centre dataset +#' are `"climatena"` for the ClimateNA gridded time series or `"cru.gpcc"` for the combined Climatic +#' Research Unit TS dataset (for temperature) and Global Precipitation Climatology Centre dataset #' (for precipitation). Defaults to `NULL`. #' @param gcms character. Vector of global climate model names. Options #' are [`list_gcms()`]. Defaults to `NULL`. @@ -52,8 +52,9 @@ #' historical scenario. See [`list_gcm_hist_years()`] for available years. #' Defaults to `NULL`. #' @template max_run +#' @template run_nm #' @param cache logical. Cache data locally? Default `TRUE` -#' @param ... other arguments passed to [`downscale_core()`]. Namely: `return_refperiod`, +#' @param ... other arguments passed to [`downscale_core()`]. Namely: `return_refperiod`, #' `vars`, `out_spatial` and `plot` #' @return `data.table` of downscaled climate variables for each location. @@ -77,7 +78,7 @@ #' ## historic observational time series #' vars <- c("PPT", "CMD", "Tave_07") #' climate_norms_hist <- downscale( -#' xyz = in_xyz, +#' xyz = in_xyz, #' which_refmap = "auto", #' return_refperiod = TRUE, #' obs_periods = "2001_2020", @@ -96,7 +97,7 @@ #' ## future projections for annual variables from three models #' climate_norms_fut <- downscale( #' xyz = in_xyz, which_refmap = "auto", -#' gcms = list_gcms()[c(1,5,6)], +#' gcms = list_gcms()[c(1, 5, 6)], #' ssps = list_ssps()[2], #' gcm_periods = list_gcm_periods()[1:2], #' # gcm_ssp_years = 2020:2060, @@ -105,21 +106,22 @@ #' ) #' #' @export -downscale <- function(xyz, which_refmap = "auto", - obs_periods = NULL, - obs_years = NULL, - obs_ts_dataset = NULL, - gcms = NULL, ssps = NULL, - gcm_periods = NULL, gcm_ssp_years = NULL, - gcm_hist_years = NULL, max_run = 0L, - cache = TRUE, ...) { +downscale <- function(xyz, which_refmap = "auto", + obs_periods = NULL, + obs_years = NULL, + obs_ts_dataset = NULL, + gcms = NULL, ssps = NULL, + gcm_periods = NULL, gcm_ssp_years = NULL, + gcm_hist_years = NULL, max_run = 0L, + run_nm = NULL, + cache = TRUE, ...) { message("Welcome to climr!") ## checks - .checkClimrDwnsclArgs( + .checkDwnsclArgs( xyz, which_refmap, obs_periods, obs_years, obs_ts_dataset, gcms, ssps, gcm_periods, gcm_ssp_years, - gcm_hist_years, max_run + gcm_hist_years, max_run, run_nm ) expectedCols <- c("lon", "lat", "elev", "id") @@ -155,15 +157,16 @@ downscale <- function(xyz, which_refmap = "auto", # xyz <- as.data.frame(xyz) message("Getting normals...") - if (which_refmap == "refmap_climatena") { - reference <- input_refmap(dbCon = dbCon, reference = "normal_na", bbox = thebb, cache = cache) - } else if (which_refmap == "refmap_prism") { - reference <- input_refmap(dbCon = dbCon, reference = "normal_bc", bbox = thebb, cache = cache) - } else if (which_refmap == "refmap_climr") { - reference <- input_refmap(dbCon = dbCon, reference = "normal_composite", bbox = thebb, cache = cache) + if(which_refmap %in% c("refmap_climatena","refmap_prism","refmap_climr")){ + reference <- input_refmap(dbCon = dbCon, reference = which_refmap, bbox = thebb, cache = cache) } else { # message("Normals not specified, using highest resolution available for each point") - bc_outline <- rast(system.file("extdata", "wna_outline.tif", package = "climr")) + rastFile <- system.file("extdata", "wna_outline.tif", package = "climr") + ## if package is loaded with devtools::load_all, file won't be found and we need to pass .libPaths + if (rastFile == "") { + rastFile <- system.file("extdata", "wna_outline.tif", package = "climr", lib.loc = .libPaths()) + } + bc_outline <- rast(rastFile) pnts <- extract(bc_outline, xyz[, .(lon, lat)], method = "simple") bc_ids <- xyz[["id"]][!is.na(pnts$PPT_01)] if (length(bc_ids) >= 1) { @@ -171,9 +174,9 @@ downscale <- function(xyz, which_refmap = "auto", xyz <- xyz[!is.na(pnts$PPT_01), ] thebb_bc <- get_bb(xyz) message("for BC...") - reference <- input_refmap(dbCon = dbCon, reference = "normal_bc", bbox = thebb_bc, cache = cache) + reference <- input_refmap(dbCon = dbCon, reference = "refmap_prism", bbox = thebb_bc, cache = cache) } else { - reference <- input_refmap(dbCon = dbCon, reference = "normal_na", bbox = thebb, cache = cache) + reference <- input_refmap(dbCon = dbCon, reference = "refmap_climatena", bbox = thebb, cache = cache) } } @@ -182,8 +185,10 @@ downscale <- function(xyz, which_refmap = "auto", obs_periods <- input_obs(dbCon, bbox = thebb, period = obs_periods, cache = cache) } if (!is.null(obs_years)) { - obs_years <- input_obs_ts(dbCon, dataset = obs_ts_dataset, - bbox = thebb, years = obs_years, cache = cache) + obs_years <- input_obs_ts(dbCon, + dataset = obs_ts_dataset, + bbox = thebb, years = obs_years, cache = cache + ) } if (!is.null(gcms)) { @@ -194,6 +199,7 @@ downscale <- function(xyz, which_refmap = "auto", ssps = ssps, period = gcm_periods, max_run = max_run, + run_nm = run_nm, cache = cache ) } else { @@ -205,7 +211,9 @@ downscale <- function(xyz, which_refmap = "auto", ssps = ssps, years = gcm_ssp_years, max_run = max_run, - cache = cache + cache = cache, + run_nm = run_nm, + fast = TRUE ) } else { gcm_ssp_ts <- NULL @@ -215,6 +223,7 @@ downscale <- function(xyz, which_refmap = "auto", bbox = thebb, gcms = gcms, years = gcm_hist_years, max_run = max_run, + run_nm = run_nm, cache = cache ) } else { @@ -224,7 +233,7 @@ downscale <- function(xyz, which_refmap = "auto", gcm_ssp_periods <- gcm_ssp_ts <- gcm_hist_ts <- NULL } - message("Downscaling!!") + message("Downscaling...") results <- downscale_core( xyz = xyz, refmap = reference, @@ -253,7 +262,7 @@ downscale <- function(xyz, which_refmap = "auto", results_na <- downscale_core( xyz = na_xyz, - reference = reference, + refmap = reference, obs = obs_periods, obs_ts = obs_years, gcms = gcm_ssp_periods, @@ -277,72 +286,112 @@ downscale <- function(xyz, which_refmap = "auto", #' #' @return NULL #' @noRd -.checkClimrDwnsclArgs <- function(xyz, which_refmap = NULL, obs_periods = NULL, obs_years = NULL, - obs_ts_dataset = NULL, gcms = NULL, ssps = list_ssps(), gcm_periods = NULL, gcm_ssp_years = NULL, - gcm_hist_years = NULL, max_run = 0L) { - if(is.null(ssps) & (!is.null(gcm_periods) | !is.null(gcm_ssp_years))){ +.checkDwnsclArgs <- function(xyz, which_refmap = NULL, obs_periods = NULL, obs_years = NULL, + obs_ts_dataset = NULL, gcms = NULL, ssps = list_ssps(), gcm_periods = NULL, gcm_ssp_years = NULL, + gcm_hist_years = NULL, max_run = 0L, run_nm = NULL) { + if (is.null(ssps) & (!is.null(gcm_periods) | !is.null(gcm_ssp_years))) { stop("ssps must be specified") } - if(!is.null(ssps)) ssps <- match.arg(ssps, list_ssps(), several.ok = TRUE) - + if(!is.null(run_nm) & max_run > 1){ + warning("max_run is > 0, but run_nm is specified. Only named runs will be returned.") + } + + if (!is.null(ssps)) { + notSupportedSsps <- setdiff(ssps, list_ssps()) + if (length(notSupportedSsps)) { + stop( + "The following SSPs are not supported:", + "\n ", paste(notSupportedSsps, collapse = ", "), + "\n Please see 'list_ssps' for list of supported SSPs." + ) + } + } + if (!is.null(which_refmap)) { which_refmap <- match.arg(which_refmap, c("auto", list_refmaps())) } - + if (!is.null(obs_periods)) { - obs_periods <- match.arg(obs_periods, list_obs_periods(), several.ok = TRUE) + notSupportedPeriods <- setdiff(obs_periods, list_obs_periods()) + if (length(notSupportedPeriods)) { + stop( + "The following observed periods are not supported:", + "\n ", paste(notSupportedPeriods, collapse = ", "), + "\n Please see 'list_obs_periods' for list of supported periods." + ) + } } - + if (!is.null(obs_years)) { - if (!all(obs_years %in% 1901:2023)) { - stop("'obs_years' must be in 1901:2023") - } + if (!all(obs_years %in% list_obs_years())) { + stop( + "'obs_years' must be in ", range(list_obs_years())[1], ":", + range(list_obs_years())[2] + ) + } } - - if(!is.null(obs_ts_dataset)){ - if(any(!obs_ts_dataset %in% c("cru.gpcc","climatena"))){ - stop("obs_ts_dataset must be cru.gpcc, climatena, or both") + + if (!is.null(obs_ts_dataset)) { + if (any(!obs_ts_dataset %in% c("cru.gpcc", "climatena"))) { + stop("obs_ts_dataset must be 'cru.gpcc', 'climatena', or 'both'") } - if(is.null(obs_years)){ + if (is.null(obs_years)) { stop("'obs_years' must be specified") } } - + if (!is.null(gcms)) { - gcms <- match.arg(gcms, list_gcms(), several.ok = TRUE) + notSupportedGCMs <- setdiff(gcms, list_gcms()) + if (length(notSupportedGCMs)) { + stop( + "The following GCMs are not supported:", + "\n ", paste(notSupportedGCMs, collapse = ", "), + "\n Please see 'list_gcms' for list of supported GCMs." + ) + } } - + if (!is.null(gcm_periods)) { - gcm_periods <- match.arg(gcm_periods, list_gcm_periods(), several.ok = TRUE) + notSupportedGCMPs <- setdiff(gcm_periods, list_gcm_periods()) + if (length(notSupportedGCMPs)) { + stop( + "The following projected periods are not supported:", + "\n ", paste(notSupportedGCMPs, collapse = ", "), + "\n Please see 'list_gcm_periods' for list of supported periods." + ) + } } - + if (!is.null(gcm_ssp_years)) { - if (!all(gcm_ssp_years %in% 2015:2100)) { - stop("'gcm_ssp_years' must be in 2015:2100") + if (!all(gcm_ssp_years %in% list_gcm_ssp_years())) { + stop( + "'gcm_ssp_years' must be in ", range(list_gcm_ssp_years())[1], ":", + range(list_gcm_ssp_years())[2] + ) } } - + msg <- "'max_run' must be 0 or larger" if (!inherits(max_run, c("integer", "numeric"))) { stop(msg) } else if (max_run < 0) { stop(msg) } - + ## check for "silly" parameter combinations if (!is.null(gcms) & - all(is.null(gcm_hist_years), is.null(gcm_ssp_years), is.null(gcm_periods), is.null(ssps))) { + all(is.null(gcm_hist_years), is.null(gcm_ssp_years), is.null(gcm_periods), is.null(ssps))) { message("'gcms' will be ignored, since 'gcm_hist_years', 'gcm_ssp_years', 'gcm_periods' and 'ssps' are missing") } - + if (is.null(gcms) & - any(!is.null(gcm_hist_years), !is.null(gcm_ssp_years), !is.null(gcm_periods), !is.null(ssps))) { + any(!is.null(gcm_hist_years), !is.null(gcm_ssp_years), !is.null(gcm_periods), !is.null(ssps))) { message("'gcms' is missing. 'gcm_hist_years', 'gcm_ssp_years', 'gcm_periods' and 'ssps' will be ignored") } - + if ((!is.null(max_run) | max_run > 0) & - is.null(gcms)) { + is.null(gcms)) { message("'gcms' is missing. 'max_run' will be ignored") } return(invisible(NULL)) diff --git a/R/downscale_core.R b/R/downscale_core.R index d1f5384..054bd87 100644 --- a/R/downscale_core.R +++ b/R/downscale_core.R @@ -3,10 +3,10 @@ #' @description #' `downscale_core()` is the engine for [`downscale()`]. #' It takes user-supplied high- and low-resolution rasters as input and downscales to user-specified point locations. -#' While less user-friendly than [`downscale()`], `downscale_core()` is more flexible in that users can supply their -#' own raster inputs. For example, a user could supply their own high-resolution climate map, instead of what is -#' available in climr, as the input to `refmap`. Another example is in downscaling a uniform warming level, as shown -#' in the example for this function. +#' While less user-friendly than [`downscale()`], `downscale_core()` is more flexible in that users can supply their +#' own raster inputs. For example, a user could supply their own high-resolution climate map, instead of what is +#' available in climr, as the input to `refmap`. Another example is in downscaling a uniform warming level, as shown +#' in the example for this function. #' #' @details #' We recommend [`downscale()`] for most purposes. @@ -31,7 +31,7 @@ #' @param out_spatial logical. Should a SpatVector be returned instead of a #' `data.frame`. #' @param plot character. If `out_spatial` is TRUE, the name of a variable to plot. -#' If the variable exists in `reference`, then its reference values will also be plotted. +#' If the variable exists in `reference`, then its reference values will also be plotted. #' Otherwise, reference January total precipitation (PPT01) values will be plotted. #' Defaults to no plotting (NULL). #' @@ -49,41 +49,43 @@ #' #' @export #' @examples -#' ## +#' ## #' library(terra) #' xyz <- data.frame(lon = runif(10, -130, -106), lat = runif(10, 37, 50), elev = runif(10), id = 1:10) #' #' ## get bounding box based on input points #' thebb <- get_bb(xyz) -#' +#' #' ## get database connection #' dbCon <- data_connect() -#' -#' # obtain the climatena 1961-1990 normals for the study area. +#' +#' # obtain the climatena 1961-1990 normals for the study area. #' refmap <- input_refmap(dbCon, thebb, reference = "refmap_climatena") -#' -#' # obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. +#' +#' # obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. #' gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], period = list_gcm_periods()[2]) -#' +#' #' # downscale the GCM data #' gcm_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = gcm_raw, vars = c("MAT", "PAS")) -#' +#' #' # create an input of uniform warming of 2 degrees Celsius and no precipitation change, for use as a null comparison to the GCM warming #' null <- gcm_raw #' use the gcm input object as a template #' names(null) <- "null_2C" -#' names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep="_")) -#' for(var in names(null[[1]])){ values(null[[1]][[var]]) <- if(length(grep("PPT", var)==1)) 1 else 2 } #' repopulate with the null values -#' +#' names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep = "_")) +#' for (var in names(null[[1]])) { +#' values(null[[1]][[var]]) <- if (length(grep("PPT", var) == 1)) 1 else 2 +#' } #' repopulate with the null values +#' #' # downscale the null values for variables of interest #' null_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = null, vars = c("MAT", "PAS")) #' pool::poolClose(dbCon) -#' +#' downscale_core <- function(xyz, refmap, gcms = NULL, obs = NULL, gcm_ssp_ts = NULL, - gcm_hist_ts = NULL, obs_ts = NULL, return_refperiod = TRUE, - vars = sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), - ppt_lr = FALSE, nthread = 1L, out_spatial = FALSE, plot = NULL) { + gcm_hist_ts = NULL, obs_ts = NULL, return_refperiod = TRUE, + vars = sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + ppt_lr = FALSE, nthread = 1L, out_spatial = FALSE, plot = NULL) { ## checks - .checkDwnsclArgs( + .checkDwnsclCoreArgs( xyz, refmap, gcms, obs, gcm_ssp_ts, gcm_hist_ts, obs_ts, return_refperiod, out_spatial, plot, vars ) @@ -363,7 +365,7 @@ downscale_ <- function(xyz, refmap, gcms, gcm_ssp_ts, gcm_hist_ts, labels <- nm normal_ <- res # Reshape (melt / dcast) to obtain final form - #ref_dt <- tstrsplit(nm, "_") + # ref_dt <- tstrsplit(nm, "_") ref_dt <- data.table(VAR = nm) # setDT(ref_dt) # setnames(ref_dt, c("VAR")) @@ -469,38 +471,40 @@ process_one_climaterast <- function(climaterast, res, xyz, timeseries = FALSE, # Cropping will reduce the size of data to load in memory climaterast <- crop(climaterast, ex, snap = "out") - gc(reset = TRUE) ## free unused memory - + gc(reset = TRUE) ## free unused memory + climaterast <- try(extract(x = climaterast, y = xyz[, .(lon, lat)], method = "bilinear")) - + ## we may have run out of memory if there are MANY rasters - ## attempt to get only unique raster cell values + ## attempt to get only unique raster cell values ## (i.e. xyz may be at higher res than the climaterast leading to extracting the same values many times) if (is(climaterast, "try-error")) { if (grepl("bad_alloc", climaterast)) { message("System is out of memory to extract climate values for the supplied coordinates") - stop("Insufficient memory to downscale climate data for these many points/climate layers.\n", - " Try reducing number of points/layers.") - } - } - + stop( + "Insufficient memory to downscale climate data for these many points/climate layers.\n", + " Try reducing number of points/layers." + ) + } + } + # else { Ceres not sure what this is for but it's always causing fails # stop("Climate value extraction failed.", # "\n Please contact developers with a reproducible example and the error:\n", - # climaterast) + # climaterast) # } - + # Create match set to match with res names - - - labels <- vapply( - strsplit(nm, "_"), - \(x) { - paste0(x[2:3], collapse = "_") - }, - character(1) - ) - + + + labels <- vapply( + strsplit(nm, "_"), + function(x) { + paste0(x[2:3], collapse = "_") + }, + character(1) + ) + if (type %in% c("obs")) { ## Create match set to match with res names labels <- nm @@ -525,12 +529,12 @@ process_one_climaterast <- function(climaterast, res, xyz, timeseries = FALSE, } setDT(ref_dt) - if (type %in% c("obs","obs_ts")) { + if (type %in% c("obs", "obs_ts")) { if (timeseries) { setnames(ref_dt, c("DATASET", "VAR", "MONTH", "PERIOD")) set(ref_dt, j = "variable", value = nm) } else { - setnames(ref_dt, c("VAR","MONTH")) + setnames(ref_dt, c("VAR", "MONTH")) set(ref_dt, j = "variable", value = nm) set(ref_dt, j = "PERIOD", value = "2001_2020") } @@ -686,7 +690,7 @@ unpackRasters <- function(ras) { if (!inherits(xyz$id, colTypes)) { stop("'xyz$id' must be an column of type ", paste(colTypes, collapse = ", ")) } - + return(xyz) } @@ -697,22 +701,29 @@ unpackRasters <- function(ras) { #' #' @return NULL #' @noRd -.checkDwnsclArgs <- function(xyz, refmap, gcms = NULL, obs = NULL, gcm_ssp_ts = NULL, gcm_hist_ts = NULL, - obs_ts = NULL, return_refperiod = FALSE, - out_spatial = FALSE, plot = NULL, vars = list_vars()) { - vars <- match.arg(vars, list_vars(), several.ok = TRUE) - +.checkDwnsclCoreArgs <- function(xyz, refmap, gcms = NULL, obs = NULL, gcm_ssp_ts = NULL, gcm_hist_ts = NULL, + obs_ts = NULL, return_refperiod = FALSE, + out_spatial = FALSE, plot = NULL, vars = list_vars()) { + notSupportedVars <- setdiff(vars, list_vars()) + if (length(notSupportedVars)) { + stop( + "The following variables are not supported:", + "\n ", paste(notSupportedVars, collapse = ", "), + "\n Please see 'list_vars' for list of supported variables." + ) + } + if (!return_refperiod %in% c(TRUE, FALSE)) { stop("'return_refperiod' must be TRUE or FALSE") } if (!out_spatial %in% c(TRUE, FALSE)) { stop("'out_spatial' must be TRUE or FALSE") } - + plot <- if (!is.null(plot)) { - match.arg(plot,list_vars()) + match.arg(plot, list_vars()) } - + if (!isTRUE(attr(refmap, "builder") == "climr")) { stop( "Please use `input_refmap` function to create `refmap`.", diff --git a/R/gcm.R b/R/gcm.R index 6a32d6a..21dd9c2 100644 --- a/R/gcm.R +++ b/R/gcm.R @@ -14,6 +14,7 @@ #' @template period #' @template max_run #' @template cache +#' @template run_nm #' #' @details #' This function returns a list with one slot for each requested GCM. Rasters inside the list contain anomalies for all requested SSPs, runs, and periods. @@ -47,32 +48,33 @@ #' pool::poolClose(dbCon) #' @rdname gcms-input-data #' @export -input_gcms <- function(dbCon, bbox = NULL, gcms = list_gcms(), ssps = list_ssps(), period = list_gcm_periods(), max_run = 0L, cache = TRUE) { +input_gcms <- function(dbCon, bbox = NULL, gcms = list_gcms(), ssps = list_ssps(), period = list_gcm_periods(), max_run = 0L, cache = TRUE, run_nm = NULL) { ## checks if (!is.null(bbox)) { .check_bb(bbox) } - + gcms <- match.arg(gcms, list_gcms(), several.ok = TRUE) ssps <- match.arg(ssps, list_ssps(), several.ok = TRUE) period <- match.arg(period, list_gcm_periods(), several.ok = TRUE) - + if (!is(max_run, "numeric")) { stop("please pass a numeric value to 'max_runs'") } - + if (!is(cache, "logical")) { stop("please pass a logical value to 'cache'") } - + # Load each file individually + select layers res <- sapply(gcms, process_one_gcm2, - ssps = ssps, period = period, - bbox = bbox, dbnames = dbnames, dbCon = dbCon, - max_run = max_run, cache = cache, USE.NAMES = TRUE, simplify = FALSE + ssps = ssps, period = period, + bbox = bbox, dbnames = dbnames, dbCon = dbCon, + max_run = max_run, cache = cache, run_nm = run_nm, + USE.NAMES = TRUE, simplify = FALSE ) attr(res, "builder") <- "climr" - + # Return a list of SpatRasters, one element for each model return(res) } @@ -89,6 +91,7 @@ input_gcms <- function(dbCon, bbox = NULL, gcms = list_gcms(), ssps = list_ssps( #' See [`list_gcm_hist_years()`] for available years. #' @template max_run #' @template cache +#' @template run_nm #' #' @seealso [list_gcm_periods()], [`list_gcm_periods()`] #' @@ -107,21 +110,21 @@ input_gcms <- function(dbCon, bbox = NULL, gcms = list_gcms(), ssps = list_ssps( #' @rdname gcms-input-data #' @export input_gcm_hist <- function(dbCon, bbox = NULL, gcms = list_gcms(), - years = 1901:2014, max_run = 0L, cache = TRUE) { + years = 1901:2014, max_run = 0L, cache = TRUE, run_nm = NULL) { ## checks if (!is.null(bbox)) { .check_bb(bbox) } - + # Load each file individually + select layers res <- sapply(gcms, process_one_gcm3, - years = years, - dbCon = dbCon, bbox = bbox, dbnames = dbnames_hist, - max_run = max_run, cache = cache, USE.NAMES = TRUE, simplify = FALSE + years = years, + dbCon = dbCon, bbox = bbox, dbnames = dbnames_hist, + max_run = max_run, cache = cache, run_nm = run_nm, USE.NAMES = TRUE, simplify = FALSE ) - res <- res[!sapply(res, is.null)] ##remove NULL + res <- res[!sapply(res, is.null)] ## remove NULL attr(res, "builder") <- "climr" - + # Return a list of SpatRasters, one element for each model return(res) } @@ -144,6 +147,8 @@ input_gcm_hist <- function(dbCon, bbox = NULL, gcms = list_gcms(), #' See [`list_gcm_ssp_years()`] for available years. #' @template max_run #' @template cache +#' @template run_nm +#' @param fast Logical. Should we use the faster method of downloading data from the database using arrays instead of Postgis rasters? #' #' @return A `list` of `SpatRasters`, each with possibly multiple layers, that can #' be used with [`downscale_core()`]. @@ -155,28 +160,41 @@ input_gcm_hist <- function(dbCon, bbox = NULL, gcms = list_gcms(), #' @importFrom terra rast writeRaster ext nlyr #' @importFrom utils head #' @importFrom RPostgres dbGetQuery +#' @importFrom dplyr tbl sql collect mutate #' @import uuid #' @import data.table +#' @import abind #' #' @rdname gcms-input-data #' @export input_gcm_ssp <- function(dbCon, bbox = NULL, gcms = list_gcms(), ssps = list_ssps(), - years = 2020:2030, max_run = 0L, cache = TRUE) { + years = 2020:2030, max_run = 0L, cache = TRUE, run_nm = NULL, fast = TRUE) { + ## checks if (!is.null(bbox)) { .check_bb(bbox) } - + if (nrow(dbnames_ts) < 1) stop("That isn't a valid GCM") # Load each file individually + select layers - res <- sapply(gcms, process_one_gcm4, - ssps = ssps, period = years, - dbnames = dbnames_ts, bbox = bbox, dbCon = dbCon, - max_run = max_run, cache = cache, USE.NAMES = TRUE, simplify = FALSE - ) + if(fast){ + res <- sapply(gcms, process_one_gcmts_fast, + ssps = ssps, period = years, + dbnames = dbnames_ts_fast, bbox = bbox, dbCon = dbCon, + max_run = max_run, cache = cache, run_nm = run_nm, USE.NAMES = TRUE, simplify = FALSE + ) + }else{ + res <- sapply(gcms, process_one_gcm4, + ssps = ssps, period = years, + dbnames = dbnames_ts, bbox = bbox, dbCon = dbCon, + max_run = max_run, cache = cache, run_nm = run_nm, USE.NAMES = TRUE, simplify = FALSE + ) + } + res <- res[!sapply(res, is.null)] ##remove NULL + attr(res, "builder") <- "climr" - + # Return a list of SpatRasters, one element for each model return(res) } @@ -223,30 +241,39 @@ list_unique <- function(files, col_num) { #' corresponding names in the PostGIS data base. See climr:::dbnames #' @template dbCon #' @template cache +#' @template run_nm #' #' @importFrom tools R_user_dir #' @importFrom data.table fread #' #' @return `SpatRaster` #' @noRd -process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbnames, dbCon, cache) { ## need to update to all GCMs +process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbnames, dbCon, cache, run_nm) { ## need to update to all GCMs gcmcode <- dbnames$dbname[dbnames$GCM == gcm_nm] # gcm_nm <- gsub("-", ".", gcm_nm) - + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") runs <- fread(file.path(rInfoPath, "gcm_periods.csv")) runs <- sort(unique(runs[mod == gcm_nm & scenario %in% ssps, run])) - sel_runs <- runs[1:(max_run + 1L)] + if(is.null(run_nm)){ + sel_runs <- runs[1:(max_run + 1L)] + }else{ + if(!run_nm %in% runs){ + stop("Run ", run_nm, "doesn't exist for this GCM.") + } + sel_runs <- run_nm + } + ## check cached needDownload <- TRUE - + cPath <- file.path(cache_path(), "gcms", gcmcode) - + if (dir.exists(cPath)) { bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) - + if (is(bnds, "try-error")) { ## try to get the data again message( @@ -257,28 +284,28 @@ process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbna needDownload <- FALSE } } - + if (!needDownload) { setorder(bnds, -numlay) - - spat_match <- lapply(1:nrow(bnds), FUN = \(x){ + + spat_match <- lapply(1:nrow(bnds), FUN = function(x){ if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] }) spat_match <- spat_match[!sapply(spat_match, is.null)] - + if (length(spat_match) > 0) { periods <- fread(file.path(cPath, "meta_period.csv")) ssps_cached <- fread(file.path(cPath, "meta_ssp.csv")) isin <- FALSE for (oldid in spat_match) { if (all(period %in% periods[uid == oldid, period]) & - all(ssps %in% ssps_cached[uid == oldid, ssps]) & - max_run <= bnds[uid == oldid, max_run]) { + all(ssps %in% ssps_cached[uid == oldid, ssps]) & + max_run <= bnds[uid == oldid, max_run]) { isin <- TRUE break } } - + if (isin) { message("Retrieving from cache...") gcm_rast <- rast(file.path(cPath, paste0(oldid, ".tif"))) @@ -297,7 +324,7 @@ process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbna needDownload <- TRUE } } - + if (needDownload) { q <- paste0( "select * from esm_layers_period where mod = '", gcm_nm, "' and scenario in ('", paste(ssps, collapse = "','"), @@ -307,14 +334,14 @@ process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbna layerinfo <- as.data.table(dbGetQuery(dbCon, q)) message("Downloading GCM anomalies") gcm_rast <- pgGetTerra(dbCon, gcmcode, tile = FALSE, bands = layerinfo$laynum, boundary = bbox) - layerinfo[,fullnm := paste(mod,var,month,scenario,run,period,sep = "_")] + layerinfo[, fullnm := paste(mod, var, month, scenario, run, period, sep = "_")] names(gcm_rast) <- layerinfo$fullnm - + if (cache) { message("Caching data...") uid <- UUIDgenerate() dir.create(cPath, recursive = TRUE, showWarnings = FALSE) - + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".tif"))) rastext <- ext(gcm_rast) t1 <- data.table( @@ -328,7 +355,7 @@ process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbna fwrite(t3, file = file.path(cPath, "meta_ssp.csv"), append = TRUE) } } - + return(gcm_rast) } @@ -343,30 +370,39 @@ process_one_gcm2 <- function(gcm_nm, ssps, bbox, period, max_run, dbnames = dbna #' @param dbnames `data.frame` with the list of available GCMs (historical projections) #' and their corresponding names in the PostGIS data base. See climr:::dbnames_hist #' @template cache +#' @template run_nm #' #' @importFrom tools R_user_dir #' @importFrom data.table fread #' #' @return `SpatRaster` #' @noRd -process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbnames_hist, cache) { ## need to update to all GCMs - if(gcm_nm %in% dbnames$GCM){ +process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbnames_hist, cache, run_nm) { ## need to update to all GCMs + if (gcm_nm %in% dbnames$GCM) { gcmcode <- dbnames$dbname[dbnames$GCM == gcm_nm] - + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") - + runs <- fread(file.path(rInfoPath, "gcm_hist.csv")) runs <- sort(unique(runs[mod == gcm_nm, run])) - sel_runs <- runs[1:(max_run + 1L)] + if(is.null(run_nm)){ + sel_runs <- runs[1:(max_run + 1L)] + }else{ + if(!run_nm %in% runs){ + stop("Run ", run_nm, "doesn't exist for this GCM.") + } + sel_runs <- run_nm + } + ## check cached needDownload <- TRUE - + cPath <- file.path(cache_path(), "gcmhist", gcmcode) - + if (dir.exists(cPath)) { bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) - + if (is(bnds, "try-error")) { ## try to get the data again message( @@ -377,26 +413,26 @@ process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbna needDownload <- FALSE } } - + if (!needDownload) { setorder(bnds, -numlay) - - spat_match <- lapply(1:nrow(bnds), FUN = \(x){ + + spat_match <- lapply(1:nrow(bnds), FUN = function(x){ if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] }) spat_match <- spat_match[!sapply(spat_match, is.null)] - + if (length(spat_match) > 0) { yeardat <- fread(file.path(cPath, "meta_years.csv")) isin <- FALSE for (oldid in spat_match) { ## see if any have all required variables if (all(years %in% yeardat[uid == oldid, years]) & - max_run <= bnds[uid == oldid, max_run]) { + max_run <= bnds[uid == oldid, max_run]) { isin <- TRUE break } } - + if (isin) { message("Retrieving from cache...") gcm_rast <- rast(file.path(cPath, paste0(oldid, ".tif"))) @@ -414,7 +450,7 @@ process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbna needDownload <- TRUE } } - + if (needDownload) { q <- paste0("select mod, var, month, run, year, laynum from esm_layers_hist where mod = '", gcm_nm, "' and year in ('", paste(years, collapse = "','"), "') and run in ('", paste(sel_runs, collapse = "','"), "')") # print(q) @@ -423,12 +459,12 @@ process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbna message("Downloading GCM anomalies") gcm_rast <- pgGetTerra(dbCon, gcmcode, tile = FALSE, bands = layerinfo$laynum, boundary = bbox) names(gcm_rast) <- layerinfo$fullnm - + if (cache) { message("Caching data...") uid <- UUIDgenerate() dir.create(cPath, recursive = TRUE, showWarnings = FALSE) - + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".tif"))) rastext <- ext(gcm_rast) t1 <- data.table( @@ -440,7 +476,7 @@ process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbna fwrite(t2, file = file.path(cPath, "meta_years.csv"), append = TRUE) } } - + return(gcm_rast) } } @@ -456,131 +492,305 @@ process_one_gcm3 <- function(gcm_nm, years, dbCon, bbox, max_run, dbnames = dbna #' @template bbox #' @template dbCon #' @template cache +#' @template run_nm #' #' @importFrom tools R_user_dir #' @importFrom data.table fread #' #' @return a `SpatRaster` #' @noRd -process_one_gcm4 <- function(gcm_nm, ssps, period, max_run, dbnames = dbnames_ts, bbox, dbCon, cache) { ## need to update to all GCMs - if(gcm_nm %in% dbnames$GCM){ - gcmcode <- dbnames$dbname[dbnames$GCM == gcm_nm] - - rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") - - runs <- fread(file.path(rInfoPath, "gcm_ts.csv")) - runs <- sort(unique(runs[mod == gcm_nm & scenario %in% ssps, run])) - if (length(runs) < 1) { - warning("That GCM isn't in our database yet.") - }else{ - sel_runs <- runs[1:(max_run + 1L)] - - ## check cached - needDownload <- TRUE - - cPath <- file.path(cache_path(), "gcmts", gcmcode) - - if (dir.exists(cPath)) { - bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) - - if (is(bnds, "try-error")) { - ## try to get the data again - message( - "Metadata file no longer exists or is unreadable.", - " Downloading the data again" +process_one_gcm4 <- function(gcm_nm, ssps, period, max_run, dbnames = dbnames_ts, bbox, dbCon, cache, run_nm) { ## need to update to all GCMs + if (gcm_nm %in% dbnames$GCM) { + gcmcode <- dbnames$dbname[dbnames$GCM == gcm_nm] + + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") + + runs <- fread(file.path(rInfoPath, "gcm_ts.csv")) + runs <- sort(unique(runs[mod == gcm_nm & scenario %in% ssps, run])) + if (length(runs) < 1) { + warning("That GCM isn't in our database yet.") + } else { + if(is.null(run_nm)){ + sel_runs <- runs[1:(max_run + 1L)] + }else{ + if(!run_nm %in% runs){ + stop("Run ", run_nm, "doesn't exist for this GCM.") + } + sel_runs <- run_nm + } + + ## check cached + needDownload <- TRUE + + cPath <- file.path(cache_path(), "gcmts", gcmcode) + + if (dir.exists(cPath)) { + bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) + + if (is(bnds, "try-error")) { + ## try to get the data again + message( + "Metadata file no longer exists or is unreadable.", + " Downloading the data again" + ) + } else { + needDownload <- FALSE + } + } + + + if (!needDownload) { + setorder(bnds, -numlay) + + spat_match <- lapply(1:nrow(bnds), FUN = function(x){ + if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] + }) + spat_match <- spat_match[!sapply(spat_match, is.null)] + + if (length(spat_match) > 0) { + periods <- fread(file.path(cPath, "meta_period.csv")) + ssps_cache <- fread(file.path(cPath, "meta_ssp.csv")) + isin <- FALSE + for (oldid in spat_match) { ## see if any have all required variables + if (all(period %in% periods[uid == oldid, period]) & + all(ssps %in% ssps_cache[uid == oldid, ssps]) & + max_run <= bnds[uid == oldid, max_run]) { + isin <- TRUE + break + } + } + + if (isin) { + message("Retrieving from cache...") + gcm_rast <- tryCatch( + suppressWarnings(rast(file.path(cPath, paste0(oldid, ".tif")))), + error = function(e) { + rast(file.path(cPath, paste0(oldid, ".grd"))) + } + ) + layinfo <- data.table(fullnm = names(gcm_rast)) + layinfo[, c("Mod", "Var", "Month", "Scenario", "Run", "Year") := tstrsplit(fullnm, "_")] + layinfo[, laynum := seq_along(fullnm)] + sel <- layinfo[Scenario %in% ssps & Year %in% period & Run %in% sel_runs, laynum] + gcm_rast <- gcm_rast[[sel]] + } else { + message("Not fully cached :( Will download more") + needDownload <- TRUE + } + } else { + message("Not fully cached :( Will download more") + needDownload <- TRUE + } + } + + if (needDownload) { + q <- paste0( + "select fullnm, laynum from esm_layers_ts where mod = '", gcm_nm, "' and scenario in ('", paste(ssps, collapse = "','"), + "') and period in ('", paste(period, collapse = "','"), "') and run in ('", paste(sel_runs, collapse = "','"), "')" ) - } else { - needDownload <- FALSE + # print(q) + layerinfo <- dbGetQuery(dbCon, q) + message("Downloading GCM anomalies") + message("Precip...") + gcm_rast_ppt <- pgGetTerra(dbCon, gsub("VAR", "ppt", gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) + names(gcm_rast_ppt) <- gsub("PPT", "PPT", layerinfo$fullnm) + message("Tmax...") + gcm_rast_tmax <- pgGetTerra(dbCon, gsub("VAR", "tmax", gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) + names(gcm_rast_tmax) <- gsub("PPT", "Tmax", layerinfo$fullnm) + message("Tmin...") + gcm_rast_tmin <- pgGetTerra(dbCon, gsub("VAR", "tmin", gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) + names(gcm_rast_tmin) <- gsub("PPT", "Tmin", layerinfo$fullnm) + gcm_rast <- c(gcm_rast_ppt, gcm_rast_tmax, gcm_rast_tmin) + + if (cache) { + message("Caching data...") + uid <- UUIDgenerate() + dir.create(cPath, recursive = TRUE, showWarnings = FALSE) + if (nlyr(gcm_rast) > 65500) { ## geotifs are limited to 65535 layers + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".grd")), filetype = "ENVI") + } else { + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".tif"))) + } + rastext <- ext(gcm_rast) + t1 <- data.table( + uid = uid, ymax = rastext[4], ymin = rastext[3], xmax = rastext[2], xmin = rastext[1], + numlay = nlyr(gcm_rast), max_run = max_run + ) + t2 <- data.table(uid = rep(uid, length(period)), period = period) + t3 <- data.table(uid = rep(uid, length(ssps)), ssps = ssps) + fwrite(t1, file = file.path(cPath, "meta_area.csv"), append = TRUE) + fwrite(t2, file = file.path(cPath, "meta_period.csv"), append = TRUE) + fwrite(t3, file = file.path(cPath, "meta_ssp.csv"), append = TRUE) + } } + + return(gcm_rast) } + } +} + + +#' Process one GCM time series at a time (faster version) +#' +#' @template gcm_nm +#' @template ssps +#' @template period +#' @template max_run +#' @param dbnames `data.frame` with the list of available GCMs (time series projections) +#' and their corresponding names in the PostGIS data base. See climr:::dbnames_ts +#' @template bbox +#' @template dbCon +#' @template cache +#' @template run_nm +#' +#' @importFrom tools R_user_dir +#' +#' @return a `SpatRaster` +#' @noRd +process_one_gcmts_fast <- function(gcm_nm, ssps, period, max_run, dbnames = dbnames_ts, bbox, dbCon, cache, run_nm) { + if(gcm_nm %in% dbnames$GCM){ + gcmcode <- dbnames$dbname[dbnames$GCM == gcm_nm] + gcmarray <- dbnames$dbarray[dbnames$GCM == gcm_nm] + rInfoPath <- file.path(R_user_dir("climr", "data"), "run_info") - - if (!needDownload) { - setorder(bnds, -numlay) + runs <- fread(file.path(rInfoPath, "gcm_ts.csv")) + runs <- sort(unique(runs[mod == gcm_nm & scenario %in% ssps, run])) + if (length(runs) < 1) { + warning("That GCM isn't in our database yet.") + }else{ + if(is.null(run_nm)){ + sel_runs <- runs[1:(max_run + 1L)] + }else{ + if(!run_nm %in% runs){ + stop("Run ", run_nm, "doesn't exist for this GCM.") + } + sel_runs <- run_nm + } - spat_match <- lapply(1:nrow(bnds), FUN = \(x){ - if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] - }) - spat_match <- spat_match[!sapply(spat_match, is.null)] - if (length(spat_match) > 0) { - periods <- fread(file.path(cPath, "meta_period.csv")) - ssps_cache <- fread(file.path(cPath, "meta_ssp.csv")) - isin <- FALSE - for (oldid in spat_match) { ## see if any have all required variables - if (all(period %in% periods[uid == oldid, period]) & - all(ssps %in% ssps_cache[uid == oldid, ssps]) & - max_run <= bnds[uid == oldid, max_run]) { - isin <- TRUE - break - } + ## check cached + needDownload <- TRUE + + cPath <- file.path(cache_path(), "gcmts", gcmcode) + + if (dir.exists(cPath)) { + bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) + + if (is(bnds, "try-error")) { + ## try to get the data again + message( + "Metadata file no longer exists or is unreadable.", + " Downloading the data again" + ) + } else { + needDownload <- FALSE } + } + + + if (!needDownload) { + setorder(bnds, -numlay) - if (isin) { - message("Retrieving from cache...") - gcm_rast <- tryCatch( - suppressWarnings(rast(file.path(cPath, paste0(oldid, ".tif")))), - error = function(e) { - rast(file.path(cPath, paste0(oldid, ".grd"))) + spat_match <- lapply(1:nrow(bnds), FUN = \(x){ + if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] + }) + spat_match <- spat_match[!sapply(spat_match, is.null)] + + if (length(spat_match) > 0) { + periods <- fread(file.path(cPath, "meta_period.csv")) + ssps_cache <- fread(file.path(cPath, "meta_ssp.csv")) + isin <- FALSE + for (oldid in spat_match) { ## see if any have all required variables + if (all(period %in% periods[uid == oldid, period]) & + all(ssps %in% ssps_cache[uid == oldid, ssps]) & + max_run <= bnds[uid == oldid, max_run]) { + isin <- TRUE + break } + } + + if (isin) { + message("Retrieving from cache...") + gcm_rast <- tryCatch( + suppressWarnings(rast(file.path(cPath, paste0(oldid, ".tif")))), + error = function(e) { + rast(file.path(cPath, paste0(oldid, ".grd"))) + } ) - layinfo <- data.table(fullnm = names(gcm_rast)) - layinfo[, c("Mod", "Var", "Month", "Scenario", "Run", "Year") := tstrsplit(fullnm, "_")] - layinfo[, laynum := seq_along(fullnm)] - sel <- layinfo[Scenario %in% ssps & Year %in% period & Run %in% sel_runs, laynum] - gcm_rast <- gcm_rast[[sel]] + layinfo <- data.table(fullnm = names(gcm_rast)) + layinfo[, c("Mod", "Var", "Month", "Scenario", "Run", "Year") := tstrsplit(fullnm, "_")] + layinfo[, laynum := seq_along(fullnm)] + sel <- layinfo[Scenario %in% ssps & Year %in% period & Run %in% sel_runs, laynum] + gcm_rast <- gcm_rast[[sel]] + } else { + message("Not fully cached :( Will download more") + needDownload <- TRUE + } } else { message("Not fully cached :( Will download more") needDownload <- TRUE } - } else { - message("Not fully cached :( Will download more") - needDownload <- TRUE } - } - - if (needDownload) { - q <- paste0( - "select fullnm, laynum from esm_layers_ts where mod = '", gcm_nm, "' and scenario in ('", paste(ssps, collapse = "','"), - "') and period in ('", paste(period, collapse = "','"), "') and run in ('", paste(sel_runs, collapse = "','"), "')" - ) - # print(q) - layerinfo <- dbGetQuery(dbCon, q) - message("Downloading GCM anomalies") - message("Precip...") - gcm_rast_ppt <- pgGetTerra(dbCon, gsub("VAR","ppt",gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) - names(gcm_rast_ppt) <- gsub("PPT", "PPT", layerinfo$fullnm) - message("Tmax...") - gcm_rast_tmax <- pgGetTerra(dbCon, gsub("VAR","tmax",gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) - names(gcm_rast_tmax) <- gsub("PPT", "Tmax", layerinfo$fullnm) - message("Tmin...") - gcm_rast_tmin <- pgGetTerra(dbCon, gsub("VAR","tmin",gcmcode), tile = FALSE, bands = layerinfo$laynum, boundary = bbox) - names(gcm_rast_tmin) <- gsub("PPT", "Tmin", layerinfo$fullnm) - gcm_rast <- c(gcm_rast_ppt, gcm_rast_tmax, gcm_rast_tmin) - if (cache) { - message("Caching data...") - uid <- UUIDgenerate() - dir.create(cPath, recursive = TRUE, showWarnings = FALSE) - if(nlyr(gcm_rast) > 65500){ ##geotifs are limited to 65535 layers - writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".grd")), filetype = "ENVI") + if (needDownload) { + template <- pgGetTerra(dbCon, name = gcmcode, tile = F, bands = 1, boundary = bbox) + + if(length(period) >= 79){ ##faster if almost all years are selected + results <- tbl(dbCon, sql(paste0("select cellid, ssp, year, run, vals from ",gcmarray," where cellid in (",paste0(values(template)[,1], collapse = ','),") + and ssp in ('",paste(ssps, collapse = "','"),"') and run in ('", paste(sel_runs, collapse = "','"), "')"))) }else{ - writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".tif"))) + results <- tbl(dbCon, sql(paste0("select cellid, ssp, year, run, vals from ",gcmarray," where cellid in (",paste0(values(template)[,1], collapse = ','),") + and year in ('",paste(period, collapse = "','"),"') and ssp in ('",paste(ssps, collapse = "','"),"') and run in ('", paste(sel_runs, collapse = "','"), "')"))) + } + + # dat <- + # results %>% + # mutate(vals = unnest(vals)) %>% + # collect() + dat <- collect(mutate(results,vals = unnest(vals))) + setDT(dat) + setorder(dat, cellid, year, ssp, run) + dat[, month := rep(sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + nrow(dat)/(12*3))] + + dat[, fullnm := paste(month, ssp, run, year, sep = "_")] + cell_nums = cells(template) + coords <- rowColFromCell(template, cell_nums) + cellcoords <- data.table(coords, values(template)) + setnames(cellcoords, c("row","col","cellid")) + dat[cellcoords, `:=`(row = i.row, col = i.col), on = "cellid"] + temp <- dat[,.(row,col,fullnm,vals)] + t2 <- dcast(temp, fullnm + row ~ col, value.var = "vals") + + t_array = split(as.data.frame(t2[,!c("fullnm","row")]), t2$fullnm) + t3 <- abind(t_array, along = 3) + gcm_rast <- rast(t3) + ext(gcm_rast) <- ext(template) + names(gcm_rast) <- paste0(gcm_nm,"_", names(t_array)) + + if (cache) { + message("Caching data...") + uid <- UUIDgenerate() + dir.create(cPath, recursive = TRUE, showWarnings = FALSE) + if(nlyr(gcm_rast) > 65500){ ##geotifs are limited to 65535 layers + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".grd")), filetype = "ENVI") + }else{ + writeRaster(gcm_rast, file.path(cPath, paste0(uid, ".tif"))) + } + rastext <- ext(gcm_rast) + t1 <- data.table( + uid = uid, ymax = rastext[4], ymin = rastext[3], xmax = rastext[2], xmin = rastext[1], + numlay = nlyr(gcm_rast), max_run = max_run + ) + t2 <- data.table(uid = rep(uid, length(period)), period = period) + t3 <- data.table(uid = rep(uid, length(ssps)), ssps = ssps) + fwrite(t1, file = file.path(cPath, "meta_area.csv"), append = TRUE) + fwrite(t2, file = file.path(cPath, "meta_period.csv"), append = TRUE) + fwrite(t3, file = file.path(cPath, "meta_ssp.csv"), append = TRUE) } - rastext <- ext(gcm_rast) - t1 <- data.table( - uid = uid, ymax = rastext[4], ymin = rastext[3], xmax = rastext[2], xmin = rastext[1], - numlay = nlyr(gcm_rast), max_run = max_run - ) - t2 <- data.table(uid = rep(uid, length(period)), period = period) - t3 <- data.table(uid = rep(uid, length(ssps)), ssps = ssps) - fwrite(t1, file = file.path(cPath, "meta_area.csv"), append = TRUE) - fwrite(t2, file = file.path(cPath, "meta_period.csv"), append = TRUE) - fwrite(t3, file = file.path(cPath, "meta_ssp.csv"), append = TRUE) } + + return(gcm_rast) } - - return(gcm_rast) - } } } diff --git a/R/globalVars.R b/R/globalVars.R index ef79b93..a8025a2 100644 --- a/R/globalVars.R +++ b/R/globalVars.R @@ -1,8 +1,11 @@ utils::globalVariables(c( - "..cols", "..expectedCols", "fullnm", "GCM", "id", "id_orig", "lat", "lon", + "..cols", "..expectedCols", "C", "Code", "DATASET", + "degree", "element1", "element2", + "fullnm", "GCM", "id", "id_orig", "lat", "lon", "laynum", "mod", "numlay", "period", "Period", "PERIOD", "Period1", "Period2", "run", "Run", "RUN", "scenario", "Scenario", "skip_if_not_installed", "SSP", - "xanom", "var", - "variables", "yanom", "Year", ".", "%>%" + "xanom", "var", "variables", "Year", "yanom", + "yeartime1", "yeartime2", + ".", "%>%" )) diff --git a/R/historic.R b/R/historic.R index 828c3f5..682f1e0 100644 --- a/R/historic.R +++ b/R/historic.R @@ -1,10 +1,10 @@ #' Retrieve observational anomalies. #' #' @description -#' `input_obs` produces anomalies of the average observed climate for a given **period**, +#' `input_obs` produces anomalies of the average observed climate for a given **period**, #' relative to the 1961-1990 reference period. The anomalies are calculated from the `"cru.gpcc"` dataset -#' which is the Climatic Research Unit TS dataset (for temperature) and Global Precipitation Climatology Centre dataset -#' (for precipitation). +#' which is the Climatic Research Unit TS dataset (for temperature) and Global Precipitation Climatology Centre dataset +#' (for precipitation). #' #' @return A `list` of `SpatRasters`, each with possibly multiple layers, that can #' be used with [`downscale_core()`]. Each element of the list corresponds to a particular period, and the @@ -35,7 +35,7 @@ input_obs <- function(dbCon, bbox = NULL, period = list_obs_periods(), cache = T if (!is.null(bbox)) { .check_bb(bbox) } - + dbnames2 <- structure(list( PERIOD = c("2001_2020"), dbname = c("historic_periods") @@ -126,8 +126,8 @@ input_obs <- function(dbCon, bbox = NULL, period = list_obs_periods(), cache = T #' `input_obs_ts` produces anomalies of observed climate for a given **time series** of individual years. #' #' @template dbCon -#' @param dataset Character. Which observational dataset to use? Options are `"climatena"` for the -#' ClimateNA gridded time series or `"cru.gpcc"` for the combined Climatic Research Unit TS dataset +#' @param dataset Character. Which observational dataset to use? Options are `"climatena"` for the +#' ClimateNA gridded time series or `"cru.gpcc"` for the combined Climatic Research Unit TS dataset #' (for temperature) and Global Precipitation Climatology Centre dataset (for precipitation). #' @template bbox #' @template cache @@ -151,15 +151,15 @@ input_obs_ts <- function(dbCon, dataset = c("cru.gpcc", "climatena"), bbox = NUL if (!is.null(bbox)) { .check_bb(bbox) } - + res <- sapply(dataset, process_one_historicts, - years = years, - dbnames = dbnames_hist_obs, bbox = bbox, dbCon = dbCon, - cache = cache, USE.NAMES = TRUE, simplify = FALSE + years = years, + dbnames = dbnames_hist_obs, bbox = bbox, dbCon = dbCon, + cache = cache, USE.NAMES = TRUE, simplify = FALSE ) - res <- res[!sapply(res, is.null)] ##remove NULL + res <- res[!sapply(res, is.null)] ## remove NULL attr(res, "builder") <- "climr" - + # Return a list of SpatRasters, one element for each model return(res) @@ -170,18 +170,18 @@ input_obs_ts <- function(dbCon, dataset = c("cru.gpcc", "climatena"), bbox = NUL } process_one_historicts <- function(dataset, years, dbCon, bbox, dbnames = dbnames_hist_obs, cache) { - if(dataset %in% dbnames$dataset){ + if (dataset %in% dbnames$dataset) { ts_name <- dataset dbcode <- dbnames$dbname[dbnames$dataset == dataset] - + ## check cached needDownload <- TRUE - + cPath <- file.path(cache_path(), "obs_ts", ts_name) - + if (dir.exists(cPath)) { bnds <- try(fread(file.path(cPath, "meta_area.csv")), silent = TRUE) - + if (is(bnds, "try-error")) { ## try to get the data again message( @@ -192,11 +192,11 @@ process_one_historicts <- function(dataset, years, dbCon, bbox, dbnames = dbname needDownload <- FALSE } } - + if (!needDownload) { setorder(bnds, -numlay) - - spat_match <- lapply(1:nrow(bnds), FUN = \(x){ + + spat_match <- lapply(1:nrow(bnds), FUN = function(x){ if (is_in_bbox(bbox, matrix(bnds[x, 2:5]))) bnds$uid[x] }) spat_match <- spat_match[!sapply(spat_match, is.null)] @@ -209,7 +209,7 @@ process_one_historicts <- function(dataset, years, dbCon, bbox, dbnames = dbname break } } - + if (isin) { message("Retrieving from cache...") hist_rast <- rast(file.path(cPath, paste0(oldid, ".tif"))) @@ -227,20 +227,20 @@ process_one_historicts <- function(dataset, years, dbCon, bbox, dbnames = dbname needDownload <- TRUE } } - + if (needDownload) { - q <- paste0("select var_nm, period, laynum from ",dbcode,"_layers where period in ('", paste(years, collapse = "','"), "')") + q <- paste0("select var_nm, period, laynum from ", dbcode, "_layers where period in ('", paste(years, collapse = "','"), "')") # print(q) layerinfo <- dbGetQuery(dbCon, q) message("Downloading obs anomalies") hist_rast <- pgGetTerra(dbCon, dbcode, tile = FALSE, bands = layerinfo$laynum, boundary = bbox) - names(hist_rast) <- paste(ts_name,layerinfo$var_nm,layerinfo$period, sep = "_") - + names(hist_rast) <- paste(ts_name, layerinfo$var_nm, layerinfo$period, sep = "_") + if (cache) { message("Caching data...") uid <- UUIDgenerate() dir.create(cPath, recursive = TRUE, showWarnings = FALSE) - + writeRaster(hist_rast, file.path(cPath, paste0(uid, ".tif"))) rastext <- ext(hist_rast) t1 <- data.table( diff --git a/R/name_changes.R b/R/name_changes.R index fc42ea2..14d6f8b 100644 --- a/R/name_changes.R +++ b/R/name_changes.R @@ -1,6 +1,6 @@ #' Changes to naming conventions applied in Version 0.1.0 #' -#' A table of correspondence between old and new names for functions, parameters, and options. +#' A table of correspondence between old and new names for functions, parameters, and options. #' For changes in climate variable names, see `data(variables)`. #' #' @format A `data.table` with columns: diff --git a/R/plot_bivariate.R b/R/plot_bivariate.R index 7a31445..6968f13 100644 --- a/R/plot_bivariate.R +++ b/R/plot_bivariate.R @@ -27,8 +27,8 @@ #' @param period_focal character. The 20-year period for which to plot the ensemble #' detail. options are `list_gcm_periods()`. #' @param ssp character. A single SSP-RCP scenario (representative concentration pathways paired with shared socioeconomic pathways). -#' Options are [`list_ssps()`]. Defaults to SSP2-4.5. -#' @param obs_period character. A single 20-year period for observed climate data. +#' Options are [`list_ssps()`]. Defaults to SSP2-4.5. +#' @param obs_period character. A single 20-year period for observed climate data. #' Options are `list_obs_periods()`. #' @inheritParams downscale #' @param legend_pos character. Position of the legend. Options are `c("bottomright", diff --git a/R/plot_timeSeries.R b/R/plot_timeSeries.R index 386fe1d..4c0f4ed 100644 --- a/R/plot_timeSeries.R +++ b/R/plot_timeSeries.R @@ -12,94 +12,94 @@ #' All global climate model anomalies are bias-corrected to the 1961-1990 reference period normals. #' #' @details -#' The input table `X` provides climate data for a single location or the average of multiple locations. +#' The input table `X` provides climate data for a single location or the average of multiple locations. #' The purpose of conducting the generation of the input table in a separate function is to allow users -#' to make multiple calls to [`plot_timeSeries()`] without needing to generate the inputs each time. -#' -#' Some combinations of `var1` and `var2` are not compatible or meaningful. -#' Examples of meaningful combinations are winter vs summer values of the same climate var or minimum vs. -#' maximum temperatures. -#' -#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr -#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they -#' don't need to be downloaded again. +#' to make multiple calls to [`plot_timeSeries()`] without needing to generate the inputs each time. #' -#' @param X A `data.table` object produced using the function [`plot_timeSeries_input()`]. This -#' table can include more models, scenarios, and variables than are used in individual calls to +#' Some combinations of `var1` and `var2` are not compatible or meaningful. +#' Examples of meaningful combinations are winter vs summer values of the same climate var or minimum vs. +#' maximum temperatures. +#' +#' Downloads of GCM time series take some time. The `plot_timeSeries_input()` function can take ~5 minutes +#' to run for the first time it is called for a location. Once the time series are cached, they +#' don't need to be downloaded again. +#' +#' @param X A `data.table` object produced using the function [`plot_timeSeries_input()`]. This +#' table can include more models, scenarios, and variables than are used in individual calls to #' [`plot_timeSeries()`]. #' @inheritParams downscale #' @param var1 character. A climate var. options are [`list_vars()`]. -#' @param var2 character. A second climate var to plot in combination with `var1`. +#' @param var2 character. A second climate var to plot in combination with `var1`. #' options are [`list_vars()`]. -#' @param showObserved logical. Plot a time series of observed climate. -#' @param showrange logical. Plot a shaded region indicating the minimum and maximum of the -#' selected ensemble of GCM simulations for each selected scenario. -#' @param yfit logical. Set the range of the y axis to the range of the visible data. If `FALSE` -#' the y axis is the range of all values of `var1` (and `var2` if applicable) in the -#' input table defined by `X`. +#' @param showObserved logical. Plot a time series of observed climate. +#' @param showrange logical. Plot a shaded region indicating the minimum and maximum of the +#' selected ensemble of GCM simulations for each selected scenario. +#' @param yfit logical. Set the range of the y axis to the range of the visible data. If `FALSE` +#' the y axis is the range of all values of `var1` (and `var2` if applicable) in the +#' input table defined by `X`. #' @param cex Numeric. The magnification factor for text size. Default is 1. -#' @param mar A numerical vector of length 4, giving the margin sizes in number of lines of text: c(bottom, left, +#' @param mar A numerical vector of length 4, giving the margin sizes in number of lines of text: c(bottom, left, #' top, right). The default is c(3,3,0.1,4). -#' @param showmean logical. Plot the ensemble mean time series. Multi-model ensemble means are -#' calculated from the mean of simulations for each model. -#' @param compile logical. Compile multiple global climate models into a multi-model ensemble. -#' If `FALSE` the single-model ensembles are plotted individually. -#' @param simplify logical. Simplify the ensemble range and mean using a smoothing spline. +#' @param showmean logical. Plot the ensemble mean time series. Multi-model ensemble means are +#' calculated from the mean of simulations for each model. +#' @param compile logical. Compile multiple global climate models into a multi-model ensemble. +#' If `FALSE` the single-model ensembles are plotted individually. +#' @param simplify logical. Simplify the ensemble range and mean using a smoothing spline. #' @param refline logical. Plot the 1961-1990 reference period mean for the selected var -#' and extend this line to the year 2100 as a visual reference. -#' @param refline.obs logical. Plot the 1961-1990 reference period mean for the observational data. -#' This should be the same as the reference line for the GCM time series. -#' @param pal character. color palette. Options are "scenario", for use when comparing scenarios, -#' and "gcms", for use when comparing GCMs. -#' @param label.endyear logical. Add a label of the final year of the observational time series. -#' @param endlabel character. Add a label to the end of each simulated time series. Options +#' and extend this line to the year 2100 as a visual reference. +#' @param refline.obs logical. Plot the 1961-1990 reference period mean for the observational data. +#' This should be the same as the reference line for the GCM time series. +#' @param pal character. color palette. Options are "scenario", for use when comparing scenarios, +#' and "gcms", for use when comparing GCMs. +#' @param label.endyear logical. Add a label of the final year of the observational time series. +#' @param endlabel character. Add a label to the end of each simulated time series. Options #' are "change", to indicate the change in year 2100 relative to the 1961-1990 baseline, or "gcms" -#' to indicate the global climate model. -#' @param yearmarkers logical. Add white points to the observational time series as a visual aid. +#' to indicate the global climate model. +#' @param yearmarkers logical. Add white points to the observational time series as a visual aid. #' @param yearlines logical. Add vertical lines on every fifth year as a visual reference -#' @param legend_pos character. Position of the legend. Viable options are c("bottomright", -#' "bottomleft", "topleft", and "topright"). +#' @param legend_pos character. Position of the legend. Viable options are `c("bottomright", +#' "bottomleft", "topleft", "topright")`. #' #' @return NULL. Draws a plot in the active graphics device. #' #' @examples -#' if(FALSE){ -#' # data frame of arbitrary points -#' my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +#' if (FALSE) { +#' # data frame of arbitrary points +#' my_points <- data.frame(lon = c(-127.7300, -127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) #' -#' # generate the input data -#' my_data <- plot_timeSeries_input(my_points) +#' # generate the input data +#' my_data <- plot_timeSeries_input(my_points) #' -#' # use the input to create a plot -#' plot_timeSeries(my_data, var1 = "Tmin_sm") +#' # use the input to create a plot +#' plot_timeSeries(my_data, var1 = "Tmin_sm") #' -#' # compare observational time series -#' plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) -#' -#' # compare mean daily minimum and maximum temperatures -#' plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") +#' # compare observational time series +#' plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) #' -#' # compare summer and winter temperatures (without simplifying the ensemble range) -#' plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) +#' # compare mean daily minimum and maximum temperatures +#' plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") #' -#' # compare global climate models -#' plot_timeSeries(my_data, gcms = list_gcms()[c(7,13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar=c(3,3,0.1,6), showObserved = FALSE) +#' # compare summer and winter temperatures (without simplifying the ensemble range) +#' plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) #' -#' # export plot to a temporary directory, including a title -#' figDir <- tempdir() -#' png( -#' filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", -#' width = 6, height = 5, pointsize = 10, res = 300 -#' ) -#' plot_timeSeries(my_data, var1 = "Tmin_sm", mar=c(3,3,2,4)) -#' title("Historical and projected summer night-time warming in the Bulkley Valley, BC") -#' dev.off() +#' # compare global climate models +#' plot_timeSeries(my_data, gcms = list_gcms()[c(7, 13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar = c(3, 3, 0.1, 6), showObserved = FALSE) +#' +#' # export plot to a temporary directory, including a title +#' figDir <- tempdir() +#' png( +#' filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", +#' width = 6, height = 5, pointsize = 10, res = 300 +#' ) +#' plot_timeSeries(my_data, var1 = "Tmin_sm", mar = c(3, 3, 2, 4)) +#' title("Historical and projected summer night-time warming in the Bulkley Valley, BC") +#' dev.off() #' } -#' #' #' @importFrom scales alpha #' @importFrom stinepack stinterp +#' @importFrom utils data +#' @importFrom graphics box #' #' @export @@ -109,299 +109,425 @@ plot_timeSeries <- function( var2 = NULL, showObserved = TRUE, obs_ts_dataset = "climatena", - gcms = list_gcms()[c(1,4,5,6,7,10,11,12)], + gcms = list_gcms()[c(1, 4, 5, 6, 7, 10, 11, 12)], ssps = list_ssps()[1:3], - showrange = TRUE, + showrange = TRUE, yfit = TRUE, cex = 1, - mar=c(3,3,0.1,4), + mar = c(3, 3, 0.1, 4), showmean = TRUE, compile = TRUE, simplify = TRUE, refline = FALSE, refline.obs = TRUE, pal = "scenario", - label.endyear = FALSE, - endlabel = "change", + label.endyear = FALSE, + endlabel = "change", yearmarkers = TRUE, yearlines = FALSE, legend_pos = "topleft") { + + ## checks if (!requireNamespace("scales", quietly = TRUE)) { stop("package scales must be installed to use this function") - } else if (!requireNamespace("stinepack", quietly = TRUE)) { + } + + if (!requireNamespace("stinepack", quietly = TRUE)) { stop("package stinepack must be installed to use this function") - } else { - data("variables", envir = environment()) - - # Scenario definitions - scenarios.selected <- c("historical", ssps) - scenarios <- c("historical", list_ssps()) - scenario.names <- c("Historical simulations", "SSP1-2.6", "SSP2-4.5", "SSP3-7.0", "SSP5-8.5") - pal.scenario <- c("gray60", "dodgerblue4", "seagreen", "darkorange3", "darkred") # these roughly match the IPCC standard scenario colours. - - # GCM color palette (from https://mk.bcgsc.ca/colorblind/) - pal.gcms <- c("#004949","#009292","#ff6db6","#ffb6db", "#490092","#006ddb","#b66dff","#6db6ff","#b6dbff", "#920000","#924900","#db6d00","#24ff24") - - # function for specifying the color - colSelect <- function(scenario, gcms){ - if(missing(gcms)){ - col = pal.scenario[which(scenarios==scenario)] - } else { - col = if(pal=="gcms") pal.gcms[which(list_gcms()==gcms)] else pal.scenario[which(scenarios==scenario)] - } - return(col) + } + + legopts <- c("bottomright", "bottomleft", "topleft", "topright") + if (!legend_pos %in% legopts) { + stop("'legend_pos' must be one of: ", + paste(legopts, collapse = ", ")) + } + + data("variables", envir = environment()) + + # Scenario definitions + scenarios.selected <- c("historical", ssps) + scenarios <- c("historical", list_ssps()) + scenario.names <- c("Historical simulations", "SSP1-2.6", "SSP2-4.5", "SSP3-7.0", "SSP5-8.5") + pal.scenario <- c("gray60", "dodgerblue4", "seagreen", "darkorange3", "darkred") # these roughly match the IPCC standard scenario colours. + + # GCM color palette (from https://mk.bcgsc.ca/colorblind/) + pal.gcms <- c("#004949", "#009292", "#ff6db6", "#ffb6db", "#490092", + "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", "#920000", + "#924900", "#db6d00", "#24ff24") + + # yeartime definitions + monthcodes <- as.character(c(paste0("0", 1:9), 10:12)) + seasonmonth.mat <- matrix(monthcodes[c(12, 1:11)], 4, byrow = TRUE) + seasons <- c("wt", "sp", "sm", "at") + season.names <- c("Winter", "Spring", "Summer", "Autumn") + yeartimes <- c(seasons, monthcodes) + yeartime.names <- c(season.names, month.name) + + # ensemble statistics definitions + ensstats <- c("ensmin", "ensmax", "ensmean") + + ## Assemble the data that will be used in the plot (for setting the ylim) + alldata <- X[, if (is.null(var2)) get(var1) else c(get(var1), get(var2))] # a vector of all potential data on the plot for setting the ylim (y axis range) + visibledata <- X[GCM %in% c(NA, gcms) & SSP %in% c(NA, ssps), + (if (is.null(var2)) get(var1) else c(get(var1), get(var2)))] # a vector of all visible data on the plot for setting the ylim (y axis range) + + # components of the var + nums <- if (is.null(var2)) 1 else 1:2 # nums is the set of var numbers (var1 or var2) (this is used later on as well) + for (num in nums) { + var <- get(paste("var", num, sep = "")) + assign(paste0("yeartime", num), as.character(variables[Code == var, "Code_Time"])) + assign(paste0("element", num), as.character(variables[Code == var, "Code_Element"])) + } + + # PLOT + par(mfrow = c(1, 1), mar = mar, mgp = c(1.75, 0.25, 0), cex = cex) + # y axis title. + if (is.null(var2)) { # if there is no second var + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"]) + } else if (element1 == element2) { # if both variables have the same element + ylab <- variables[Code == var1, "Element"] + } else if (yeartime1 == yeartime2) { # if both variables have the same yeartime + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"], "or", variables[Code == var2, "Element"]) + } else { # if variables1 and 2 have different elements and yeartimes + ylab <- paste(yeartime.names[which(yeartimes == yeartime1)], variables[Code == var1, "Element"], "or", variables[Code == var2, "Element"]) + } + plot(0, col = "white", xlim = c(1900, 2100), + ylim = range(if (yfit) visibledata else alldata, na.rm = TRUE), + xaxs = "i", xaxt = "n", tck = 0, xlab = "", ylab = ylab) + axis(1, at = seq(1850, 2100, 25), labels = seq(1850, 2100, 25), tck = 0) + + num <- 1 + for (num in nums) { + yeartime <- get(paste("yeartime", num, sep = "")) + element <- get(paste("element", num, sep = "")) + var <- get(paste("var", num, sep = "")) + + if (compile) { # this plots a single envelope for the ensemble as a whole + temp.data <- X[GCM %in% gcms, c("PERIOD", "SSP", "RUN", var), with = FALSE] + plot.ensemble(temp.data, + var = var, var2 = var2, + refline = refline, showmean = showmean, + endlabel = endlabel, element = element, + element1 = element1, element2 = element2, + compile = compile, yeartime.names = yeartime.names, + yeartimes = yeartimes, yeartime = yeartime, + gcm = NULL, pal = pal, pal.scenario = pal.scenario, + pal.gcms = pal.gcms, + scenarios.selected = scenarios.selected, scenarios = scenarios, + showrange = showrange, simplify = simplify) + } else { + for (gcm in gcms) { # this plots of individual GCM ensembles. + temp.data <- X[GCM == gcm, c("PERIOD", "SSP", "RUN", var), with = FALSE] + plot.ensemble(temp.data, + var = var, var2 = var2, + refline = refline, showmean = showmean, + endlabel = endlabel, element = element, + element1 = element1, element2 = element2, + compile = compile, yeartime.names = yeartime.names, + yeartimes = yeartimes, yeartime = yeartime, + gcm = gcm, pal = pal, pal.scenario = pal.scenario, + pal.gcms = pal.gcms, + scenarios.selected = scenarios.selected, scenarios = scenarios, + showrange = showrange, simplify = simplify) } - - # yeartime definitions - monthcodes <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12") - seasonmonth.mat <- matrix(monthcodes[c(12, 1:11)],4, byrow=TRUE) - seasons <- c("wt", "sp", "sm", "at") - season.names <- c("Winter", "Spring", "Summer", "Autumn") - yeartimes <- c(seasons, monthcodes) - yeartime.names <- c(season.names, month.name) - - # ensemble statistics definitions - ensstats <- c("ensmin", "ensmax", "ensmean") - - ## Assemble the data that will be used in the plot (for setting the ylim) - alldata <- X[, if(is.null(var2)) get(var1) else c(get(var1), get(var2))] # a vector of all potential data on the plot for setting the ylim (y axis range) - visibledata <- X[GCM%in%c(NA, gcms) & SSP%in%c(NA, ssps), (if(is.null(var2)) get(var1) else c(get(var1), get(var2)))] # a vector of all visible data on the plot for setting the ylim (y axis range) - - # components of the var - nums <- if(is.null(var2)) 1 else 1:2 #nums is the set of var numbers (var1 or var2) (this is used later on as well) - for(num in nums){ - var <- get(paste("var",num,sep="")) - assign(paste0("yeartime", num), as.character(variables[Code == var, "Code_Time"]) ) - assign(paste0("element", num), as.character(variables[Code == var, "Code_Element"]) ) + } + + # overlay the 5-year lines on top of all polygons + if (yearlines) { + for (n in seq(1905, 2095, 5)) { + lines(c(n, n), c(-9999, 9999), col = "grey", lty = 2) } - - # PLOT - par(mfrow=c(1,1), mar=mar, mgp=c(1.75, 0.25, 0), cex=cex) - # y axis title. - if(is.null(var2)){ #if there is no second var - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"]) - } else - if(element1==element2){ #if both variables have the same element - ylab <- variables[Code==var1, "Element"] - } else - if(yeartime1==yeartime2){ #if both variables have the same yeartime - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"], "or", variables[Code==var2, "Element"]) - } else { #if variables1 and 2 have different elements and yeartimes - ylab <- paste(yeartime.names[which(yeartimes==yeartime1)], variables[Code==var1, "Element"], "or", variables[Code==var2, "Element"]) - } - plot(0, col="white", xlim=c(1900, 2100), ylim=range(if(yfit) visibledata else alldata, na.rm = TRUE), xaxs="i", xaxt="n", tck=0, xlab="", ylab=ylab) - axis(1, at=seq(1850,2100,25), labels = seq(1850,2100,25), tck=0) - - num <- 1 - for(num in nums){ - yeartime <- get(paste("yeartime",num,sep="")) - element <- get(paste("element",num,sep="")) - var <- get(paste("var",num,sep="")) + } + + if (showObserved) { + # add in observations + obs.colors <- c("black", "blue", "red") + obs.options <- c("climatena", "cru.gpcc") ## , "era5" + for (obs.dataset in obs_ts_dataset) { # TODO update this code block once i know how the datasets are identified in the climr output + obs.color <- obs.colors[which(obs.options == obs.dataset)] + x.obs <- as.numeric(X[DATASET == obs.dataset & PERIOD %in% 1900:2100, "PERIOD"][[1]]) + y.obs <- X[DATASET == obs.dataset & PERIOD %in% 1900:2100, get(var)] + recent.obs <- mean(y.obs[which(x.obs %in% 2014:2023)], na.rm = TRUE) + baseline.obs <- mean(y.obs[which(x.obs %in% 1961:1990)], na.rm = TRUE) + end <- max(which(!is.na(y.obs))) + lines(x.obs[which(x.obs < 1951)], y.obs[which(x.obs < 1951)], lwd = 3, lty = 3, col = obs.color) + lines(x.obs[which(x.obs > 1949)], y.obs[which(x.obs > 1949)], lwd = 4, col = obs.color) - # function for plotting time series for gcms or compiled ensemble - plot.ensemble <- function(x) { - # scenario <- scenarios.selected[1] - temp.historical <- x[is.na(SSP), c("PERIOD", "RUN", var), with=FALSE] - x.historical <- as.numeric(temp.historical[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) - ensmin.historical <- temp.historical[RUN!="ensembleMean", .(min = min(get(var), na.rm=TRUE)), by = PERIOD][["min"]] - ensmax.historical <- temp.historical[RUN!="ensembleMean", .(max = max(get(var), na.rm=TRUE)), by = PERIOD][["max"]] - ensmean.historical <- temp.historical[RUN=="ensembleMean", .(mean = mean(get(var), na.rm=TRUE)), by = PERIOD][["mean"]] #calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs - - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])][-1]){ - temp <- x[SSP==scenario, c("PERIOD", "RUN", var), with=FALSE] - x.temp <- as.numeric(temp[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) - ensmin.temp <- temp[RUN!="ensembleMean", .(min = min(get(var), na.rm=TRUE)), by = PERIOD][["min"]] - ensmax.temp <- temp[RUN!="ensembleMean", .(max = max(get(var), na.rm=TRUE)), by = PERIOD][["max"]] - ensmean.temp <- temp[RUN=="ensembleMean", .(mean = mean(get(var), na.rm=TRUE)), by = PERIOD][["mean"]] #calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs - assign(paste("x", scenario, sep="."), c(x.historical[length(x.historical)], x.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmin", scenario, sep="."), c(ensmin.historical[length(ensmin.historical)], ensmin.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmax", scenario, sep="."), c(ensmax.historical[length(ensmax.historical)], ensmax.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - assign(paste("ensmean", scenario, sep="."), c(ensmean.historical[length(ensmean.historical)], ensmean.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. - } - - if(showrange) { - if(simplify==FALSE){ - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])]){ - x3 <- get(paste("x", scenario, sep=".")) - polygon(c(x3, rev(x3)), c(get(paste("ensmin", scenario, sep=".")), rev(get(paste("ensmax", scenario, sep=".")))), col=alpha(colSelect(scenario, gcms), 0.35), border=colSelect(scenario, gcms)) - } - } else { - scenarios.select <- scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])][-1] - for(scenario in scenarios.select){ - - if(scenario==scenarios.select[1]){ # we need to run spline through the historical/projected transition - x4 <- c(x.historical, get(paste("x", scenario, sep="."))[-1]) - y.ensmin <- c(ensmin.historical, get(paste("ensmin", scenario, sep="."))[-1]) - y.ensmax <- c(ensmax.historical, get(paste("ensmax", scenario, sep="."))[-1]) - s.ensmin <- smooth.spline(x4,y.ensmin, df=8) - s.ensmax <- smooth.spline(x4,y.ensmax, df=8) - subset.hist <- which(x4%in%x.historical) - subset.proj <- which(x4%in%get(paste("x", scenario, sep="."))) - polygon(c(s.ensmin$x[subset.hist], rev(s.ensmax$x[subset.hist])), c(s.ensmin$y[subset.hist], rev(s.ensmax$y[subset.hist])), col=alpha(if(pal=="gcms") colSelect(scenario, gcms) else pal.scenario[which(scenarios=="historical")], 0.35), border=NA) - polygon(c(s.ensmin$x[subset.proj], rev(s.ensmax$x[subset.proj])), c(s.ensmin$y[subset.proj], rev(s.ensmax$y[subset.proj])), col=alpha(colSelect(scenario, gcms), 0.35), border=NA) - } else { # this second routine uses interpolation splines so that the starting point for all scenarios is the same - x5 <- c(x.historical, get(paste("x", scenario, sep="."))[-1]) - y.ensmin2 <- c(ensmin.historical, get(paste("ensmin", scenario, sep="."))[-1]) - y.ensmax2 <- c(ensmax.historical, get(paste("ensmax", scenario, sep="."))[-1]) - s.ensmin2 <- smooth.spline(x5,y.ensmin2, df=8) - s.ensmax2 <- smooth.spline(x5,y.ensmax2, df=8) - knots.hist <- which(x5%in%c(seq(1860, 2000, 20), 2014)) - knots.proj <- which(x5%in%c(seq(2030, 2090, 20), 2100)) - s.ensmin3 <- stinterp(x5[c(knots.hist, knots.proj)],c(s.ensmin$y[knots.hist], s.ensmin2$y[knots.proj]), x5) - s.ensmax3 <- stinterp(x5[c(knots.hist, knots.proj)],c(s.ensmax$y[knots.hist], s.ensmax2$y[knots.proj]), x5) - polygon(c(s.ensmin3$x[subset.proj], rev(s.ensmax3$x[subset.proj])), c(s.ensmin3$y[subset.proj], rev(s.ensmax3$y[subset.proj])), col=alpha(colSelect(scenario, gcms), 0.35), border=NA) - } - } - } - } - - if(refline){ - ref.temp <- mean(ensmean.historical[which(x.historical%in%1961:1990)]) - lines(1961:1990, rep(ref.temp, 30), lwd=2) - lines(c(1990,2100), rep(ref.temp, 2), lty=2) - } - - # overlay the ensemble mean lines on top of all polygons - for(scenario in scenarios.selected[order(c(1,4,5,3,2)[which(scenarios%in%scenarios.selected)])]){ - - # calculate a spline through the time series (used for plotting and the text warming value) - if(scenario=="historical"){ # need to run spline through the historical/projected transition - x4 <- c(x.historical, get(paste("x", scenarios.selected[2], sep="."))) - y4 <- c(ensmean.historical, get(paste("ensmean", scenarios.selected[2], sep="."))) - } else { - x4 <- c(x.historical, get(paste("x", scenario, sep="."))) - y4 <- c(ensmean.historical, get(paste("ensmean", scenario, sep="."))) - } - s4 <- smooth.spline(x4,y4, df=10) - subset <- which(x4%in%get(paste("x", scenario, sep="."))) - - # plot the ensemble mean - if(showmean){ - if(simplify){ - lines(x=s4$x[subset], y=s4$y[subset], col=colSelect(scenario, gcms), lwd=2) - } else lines(x=get(paste("x", scenario, sep=".")), y=get(paste("ensmean", scenario, sep=".")), col=colSelect(scenario, gcms), lwd=2) - } - - # end labels - if(is.null(endlabel)==FALSE){ - if(scenario != "historical"){ - par(xpd=TRUE) - baseline <- mean(ensmean.historical[which(x.historical%in%1961:1990)]) - projected <- s4$y[length(s4$y)] - if(endlabel=="change"){ - if(element%in%c("PPT", "PAS", "CMD", "MAP", "MSP", "DDsub18", "DD18", "DDsub0", "DD5", "Eref")){ - change <- round(projected/baseline-1,2) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change*100,"%", sep="") else paste(change*100,"%", sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else if(element%in%c("Tave", "Tmin", "Tmax", "MCMT", "MWMT", "EXT", "EMT", "MAT")){ - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) bquote("+" * .(change) * degree * C) else bquote(.(change) * degree * C), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else if(element%in%c("RH")){ - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change,"%", sep="") else paste(change,"%", sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } else { - change <- round(projected-baseline,1) - if(is.na(change)==FALSE) text(2099,projected, if(change>0) paste("+",change, sep="") else paste(change, sep=""), col=colSelect(scenario, gcms), pos=4, font=2, cex=1) - } - } - if(endlabel=="gcms"){ - text(2099,projected, if(compile) "ensemble" else gcms, col=colSelect(scenario, gcms), pos=4, font=1, cex=1) - } - par(xpd=FALSE) - } - } - } - - # Text to identify the time of year - if(!is.null(var2)){ #if there is no second var - if(element1==element2){ - label <- yeartime.names[which(yeartimes==yeartime)] - } else { - label <- paste(yeartime.names[which(yeartimes==yeartime)], element) - } - temp <- get("ensmax.historical") - text(1925,mean(temp[10:40]), label, col="black", pos=3, font=2, cex=1) - } + if (yearmarkers) { + points(x.obs[which(x.obs > 1949)], + y.obs[which(x.obs > 1949)], + pch = 21, bg = "white", col = obs.color, cex = 0.4) + points(x.obs[which(x.obs > 1949)[seq(1, 999, 5)]], + y.obs[which(x.obs > 1949)[seq(1, 999, 5)]], + pch = 21, bg = "white", col = obs.color, cex = 0.7) + } + + if (label.endyear) { + points(x.obs[end], y.obs[end], pch = 16, cex = 1, col = obs.color) + text(x.obs[end], y.obs[end], x.obs[end], pos = 4, offset = 0.25, + col = obs.color, cex = 1, ) } - if(compile){ #this plots a single envelope for the ensemble as a whole - temp.data <- X[GCM%in%gcms, c("PERIOD", "SSP", "RUN", var), with=FALSE] - plot.ensemble(temp.data) + if (refline.obs) { + lines(1961:1990, rep(baseline.obs, 30), lwd = 1, col = obs.color) + lines(c(1990, 2100), rep(baseline.obs, 2), lty = 2, col = obs.color) + } + } + } + print(num) + } + + if (showObserved) { + # Sources legend + a <- if ("climatena" %in% obs_ts_dataset) 1 else NA + b <- if ("cru.gpcc" %in% obs_ts_dataset) 2 else NA + c <- if ("era5" %in% obs_ts_dataset) 3 else NA + d <- if (length(gcms > 0)) 4 else NA + s <- !is.na(c(a, b, c, d)) + legend.GCM <- if (length(gcms) > 1) { + paste("Simulated (", length(gcms), " GCMs)", sep = "") + } else { + paste("Simulated (", gcms, ")", sep = "") + } + legend(legend_pos, + title = "", + legend = c("Observed (ClimateNA)", "Observed (CRU/GPCC)", "Observed (ERA5)", legend.GCM)[s], + bty = "n", + lty = rep(1, 4)[s], + col = c(obs.colors, "gray")[s], + lwd = c(4, 4, 4, 2)[s], + pch = rep(NA, 4)[s], + pt.bg = rep(NA, 4)[s], + pt.cex = rep(NA, 4)[s] + ) + } + + # Scenario legend + if (pal == "gcms") { + s <- which(list_gcms() %in% gcms) + legend(ifelse(grepl("top", legend_pos), "top", "bottom"), + title = "GCMs", legend = gcms, bty = "n", + col = pal.gcms[s], pch = 22, pt.bg = alpha(pal.gcms[s], 0.35), pt.cex = 2 + ) + } else { + s <- rev(which(scenarios[-1] %in% scenarios.selected)) + legend(ifelse(grepl("top", legend_pos), "top", "bottom"), + title = "Scenarios", legend = c("Historical", scenario.names[-1][s]), bty = "n", + lty = rep(NA, 5)[c(1, s + 1)], col = pal.scenario[c(1, s + 1)], + lwd = rep(NA, 5)[c(1, s + 1)], pch = rep(22, 5)[c(1, s + 1)], + pt.bg = alpha(pal.scenario[c(1, s + 1)], 0.35), + pt.cex = rep(2, 5)[c(1, s + 1)] + ) + } + + # mtext(paste(" Created using climr (https://bcgov.github.io/climr/)"), side=1, line=1.5, adj=0.0, font=1, cex=1.1, col="gray") + + box() +} + + +#' Function for plotting time series for gcms or compiled ensemble +#' +#' @param x climate data +#' @param var TODO +#' @param scenarios.selected TODO +#' @param scenarios TODO +#' @param gcm TODO +#' @param pal.scenario TODO +#' @param pal.gcms TODO +#' @param element TODO +#' @param element1 TODO +#' @param element2 TODO +#' @param yeartime.names TODO +#' @param yeartimes TODO +#' @param yeartime TODO +#' @inheritParams plot_timeSeries +#' +#' @importFrom graphics polygon +#' @importFrom stats smooth.spline +plot.ensemble <- function(x, var, scenarios.selected, scenarios, + showrange = TRUE, simplify = TRUE, gcm = NULL, + pal, pal.scenario, pal.gcms, refline = FALSE, showmean = TRUE, + endlabel = "change", element, + compile = TRUE, var2 = NULL, element1, element2, + yeartime.names, yeartimes, yeartime) { + # scenario <- scenarios.selected[1] + temp.historical <- x[is.na(SSP), c("PERIOD", "RUN", var), with = FALSE] + x.historical <- as.numeric(temp.historical[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) + ensmin.historical <- temp.historical[RUN != "ensembleMean", .(min = min(get(var), na.rm = TRUE)), by = PERIOD][["min"]] + ensmax.historical <- temp.historical[RUN != "ensembleMean", .(max = max(get(var), na.rm = TRUE)), by = PERIOD][["max"]] + ensmean.historical <- temp.historical[RUN == "ensembleMean", .(mean = mean(get(var), na.rm = TRUE)), by = PERIOD][["mean"]] # calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs + + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])][-1]) { + temp <- x[SSP == scenario, c("PERIOD", "RUN", var), with = FALSE] + x.temp <- as.numeric(temp[, .(min = min(get(var))), by = PERIOD][["PERIOD"]]) + ensmin.temp <- temp[RUN != "ensembleMean", .(min = min(get(var), na.rm = TRUE)), by = PERIOD][["min"]] + ensmax.temp <- temp[RUN != "ensembleMean", .(max = max(get(var), na.rm = TRUE)), by = PERIOD][["max"]] + ensmean.temp <- temp[RUN == "ensembleMean", .(mean = mean(get(var), na.rm = TRUE)), by = PERIOD][["mean"]] # calculate mean using the single-model ensembleMean, to ensure that the mean isn't biased towards models with more runs + assign(paste("x", scenario, sep = "."), c(x.historical[length(x.historical)], x.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmin", scenario, sep = "."), c(ensmin.historical[length(ensmin.historical)], ensmin.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmax", scenario, sep = "."), c(ensmax.historical[length(ensmax.historical)], ensmax.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + assign(paste("ensmean", scenario, sep = "."), c(ensmean.historical[length(ensmean.historical)], ensmean.temp)) # add last year of historical runs. this is necessary to get a seamless transition from the historical polygon to the future polygon. + } + + if (showrange) { + if (isFALSE(simplify)) { + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])]) { + x3 <- get(paste("x", scenario, sep = ".")) + colSel <- colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) + polygon(c(x3, rev(x3)), + c(get(paste("ensmin", scenario, sep = ".")), rev(get(paste("ensmax", scenario, sep = ".")))), + col = alpha(colSel, 0.35), + border = colSel + ) + } + } else { + scenarios.select <- scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])][-1] + for (scenario in scenarios.select) { + if (scenario == scenarios.select[1]) { # we need to run spline through the historical/projected transition + x4 <- c(x.historical, get(paste("x", scenario, sep = "."))[-1]) + y.ensmin <- c(ensmin.historical, get(paste("ensmin", scenario, sep = "."))[-1]) + y.ensmax <- c(ensmax.historical, get(paste("ensmax", scenario, sep = "."))[-1]) + s.ensmin <- smooth.spline(x4, y.ensmin, df = 8) + s.ensmax <- smooth.spline(x4, y.ensmax, df = 8) + subset.hist <- which(x4 %in% x.historical) + subset.proj <- which(x4 %in% get(paste("x", scenario, sep = "."))) - } else for(gcm in gcms){ #this plots of individual GCM ensembles. - temp.data <- X[GCM==gcm, c("PERIOD", "SSP", "RUN", var), with=FALSE] - plot.ensemble(temp.data) + colSel <- colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) - print(gcms) + polygon(c(s.ensmin$x[subset.hist], rev(s.ensmax$x[subset.hist])), + c(s.ensmin$y[subset.hist], rev(s.ensmax$y[subset.hist])), + col = alpha(ifelse(pal == "gcms", colSel, pal.scenario[which(scenarios == "historical")]), 0.35), + border = NA + ) + polygon(c(s.ensmin$x[subset.proj], rev(s.ensmax$x[subset.proj])), + c(s.ensmin$y[subset.proj], rev(s.ensmax$y[subset.proj])), + col = alpha(colSel, 0.35), + border = NA + ) + } else { # this second routine uses interpolation splines so that the starting point for all scenarios is the same + x5 <- c(x.historical, get(paste("x", scenario, sep = "."))[-1]) + y.ensmin2 <- c(ensmin.historical, get(paste("ensmin", scenario, sep = "."))[-1]) + y.ensmax2 <- c(ensmax.historical, get(paste("ensmax", scenario, sep = "."))[-1]) + s.ensmin2 <- smooth.spline(x5, y.ensmin2, df = 8) + s.ensmax2 <- smooth.spline(x5, y.ensmax2, df = 8) + knots.hist <- which(x5 %in% c(seq(1860, 2000, 20), 2014)) + knots.proj <- which(x5 %in% c(seq(2030, 2090, 20), 2100)) + s.ensmin3 <- stinterp(x5[c(knots.hist, knots.proj)], c(s.ensmin$y[knots.hist], s.ensmin2$y[knots.proj]), x5) + s.ensmax3 <- stinterp(x5[c(knots.hist, knots.proj)], c(s.ensmax$y[knots.hist], s.ensmax2$y[knots.proj]), x5) + + colSel <- colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) + + polygon(c(s.ensmin3$x[subset.proj], rev(s.ensmax3$x[subset.proj])), + c(s.ensmin3$y[subset.proj], rev(s.ensmax3$y[subset.proj])), + col = alpha(colSel, 0.35), + border = NA + ) } + } + } + } + + if (refline) { + ref.temp <- mean(ensmean.historical[which(x.historical %in% 1961:1990)]) + lines(1961:1990, rep(ref.temp, 30), lwd = 2) + lines(c(1990, 2100), rep(ref.temp, 2), lty = 2) + } + + # overlay the ensemble mean lines on top of all polygons + for (scenario in scenarios.selected[order(c(1, 4, 5, 3, 2)[which(scenarios %in% scenarios.selected)])]) { + # calculate a spline through the time series (used for plotting and the text warming value) + if (scenario == "historical") { # need to run spline through the historical/projected transition + x4 <- c(x.historical, get(paste("x", scenarios.selected[2], sep = "."))) + y4 <- c(ensmean.historical, get(paste("ensmean", scenarios.selected[2], sep = "."))) + } else { + x4 <- c(x.historical, get(paste("x", scenario, sep = "."))) + y4 <- c(ensmean.historical, get(paste("ensmean", scenario, sep = "."))) + } + s4 <- smooth.spline(x4, y4, df = 10) + subset <- which(x4 %in% get(paste("x", scenario, sep = "."))) + + # plot the ensemble mean + colSel <- colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) + if (showmean) { + if (simplify) { + lines(x = s4$x[subset], y = s4$y[subset], col = colSel, lwd = 2) + } else { + lines(x = get(paste("x", scenario, sep = ".")), + y = get(paste("ensmean", scenario, sep = ".")), + col = colSel, lwd = 2) + } + } + + # end labels + if (!is.null(endlabel)) { + if (scenario != "historical") { + par(xpd = TRUE) + baseline <- mean(ensmean.historical[which(x.historical %in% 1961:1990)]) + projected <- s4$y[length(s4$y)] - # overlay the 5-year lines on top of all polygons - if(yearlines){ - for(n in seq(1905, 2095, 5)){ - lines(c(n, n), c(-9999, 9999), col="grey", lty=2) - } - } + colSel <- colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) - if(showObserved){ - # add in observations - obs.colors <- c("black", "blue", "red") - obs.options <- c("climatena", "cru.gpcc") ##, "era5" - for(obs.dataset in obs_ts_dataset){ #TODO update this code block once i know how the datasets are identified in the climr output - obs.color <- obs.colors[which(obs.options==obs.dataset)] - x.obs <- as.numeric(X[DATASET==obs.dataset & PERIOD%in%1900:2100, "PERIOD"][[1]]) - y.obs <- X[DATASET==obs.dataset & PERIOD%in%1900:2100, get(var)] - recent.obs <- mean(y.obs[which(x.obs%in%2014:2023)], na.rm=TRUE) - baseline.obs <- mean(y.obs[which(x.obs%in%1961:1990)], na.rm=TRUE) - end <- max(which(!is.na(y.obs))) - lines(x.obs[which(x.obs<1951)], y.obs[which(x.obs<1951)], lwd=3, lty=3, col=obs.color) - lines(x.obs[which(x.obs>1949)], y.obs[which(x.obs>1949)], lwd=4, col=obs.color) - if(yearmarkers){ - points(x.obs[which(x.obs>1949)], y.obs[which(x.obs>1949)], pch=21, bg="white", col=obs.color, cex=0.4) - points(x.obs[which(x.obs>1949)[seq(1,999,5)]], y.obs[which(x.obs>1949)[seq(1,999,5)]], pch=21, bg="white", col=obs.color, cex=0.7) + if (endlabel == "change") { + if (element %in% c("PPT", "PAS", "CMD", "MAP", "MSP", "DDsub18", "DD18", "DDsub0", "DD5", "Eref")) { + change <- round(projected / baseline - 1, 2) + if (is.na(change) == FALSE) { + txt <- ifelse(change > 0, paste0("+", change * 100, "%"), paste0(change * 100, "%")) + text(2099, projected, txt, col = colSel, pos = 4, font = 2, cex = 1) } - if(label.endyear){ - points(x.obs[end], y.obs[end], pch=16, cex=1, col=obs.color) - text(x.obs[end], y.obs[end], x.obs[end], pos= 4, offset = 0.25, col=obs.color, cex=1,) + } else if (element %in% c("Tave", "Tmin", "Tmax", "MCMT", "MWMT", "EXT", "EMT", "MAT")) { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) { + ## ifelse won't work with bquote. + txt <- if (change > 0) bquote("+" * .(change) * degree * C) else bquote(.(change) * degree * C) + text(2099, projected, txt, col = colSel, pos = 4, font = 2, cex = 1) } - if(refline.obs){ - lines(1961:1990, rep(baseline.obs, 30), lwd=1, col=obs.color) - lines(c(1990,2100), rep(baseline.obs, 2), lty=2, col=obs.color) + } else if (element %in% c("RH")) { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) { + txt <- ifelse(change > 0, paste0("+", change, "%"), paste(change, "%")) + text(2099, projected, txt, col = colSel, pos = 4, font = 2, cex = 1) + } + } else { + change <- round(projected - baseline, 1) + if (is.na(change) == FALSE) { + txt <- ifelse(change > 0, paste0("+", change), paste0(change)) + text(2099, projected, txt, col = colSel, pos = 4, font = 2, cex = 1) } } } - print(num) - } - - if(showObserved){ - # Sources legend - a <- if("climatena"%in%obs_ts_dataset) 1 else NA - b <- if("cru.gpcc"%in%obs_ts_dataset) 2 else NA - c <- if("era5"%in%obs_ts_dataset) 3 else NA - d <- if(length(gcms>0)) 4 else NA - s <- !is.na(c(a,b,c,d)) - legend.GCM <- if(length(gcms)>1) paste("Simulated (", length(gcms), " GCMs)", sep="") else paste("Simulated (", gcms, ")", sep="") - legend(legend_pos, title = "", legend=c("Observed (ClimateNA)", "Observed (CRU/GPCC)", "Observed (ERA5)", legend.GCM)[s], bty="n", - lty=c(1,1,1,1)[s], - col=c(obs.colors, "gray")[s], - lwd=c(4,4,4,2)[s], - pch=c(NA,NA,NA,NA)[s], - pt.bg = c(NA, NA,NA,NA)[s], - pt.cex=c(NA,NA,NA,NA)[s]) - } - - #Scenario legend - if(pal=="gcms"){ - s <- which(list_gcms()%in%gcms) - legend(c("top", "bottom")[if(length(grep("top", legend_pos))==1) 1 else 2], title = "GCMs", legend=gcms, bty="n", - col=pal.gcms[s], pch=22, pt.bg = alpha(pal.gcms[s], 0.35), pt.cex=2) - } else { - s <- rev(which(scenarios[-1]%in%scenarios.selected)) - legend(c("top", "bottom")[if(length(grep("top", legend_pos))==1) 1 else 2], title = "Scenarios", legend=c("Historical", scenario.names[-1][s]), bty="n", - lty=c(NA,NA,NA,NA,NA)[c(1,s+1)], col=pal.scenario[c(1,s+1)], lwd=c(NA,NA,NA,NA,NA)[c(1,s+1)], pch=c(22,22,22,22,22)[c(1,s+1)], pt.bg = alpha(pal.scenario[c(1,s+1)], 0.35), pt.cex=c(2,2,2,2,2)[c(1,s+1)]) + if (endlabel == "gcms") { + txt <- ifelse(compile, "ensemble", gcm) + text(2099, projected, txt, col = colSel, pos = 4, font = 1, cex = 1) + } + par(xpd = FALSE) } - - # mtext(paste(" Created using climr (https://bcgov.github.io/climr/)"), side=1, line=1.5, adj=0.0, font=1, cex=1.1, col="gray") - - box() - } + } + + # Text to identify the time of year + if (!is.null(var2)) { # if there is no second var + if (element1 == element2) { + label <- yeartime.names[which(yeartimes == yeartime)] + } else { + label <- paste(yeartime.names[which(yeartimes == yeartime)], element) + } + temp <- get("ensmax.historical") + text(1925, mean(temp[10:40]), label, col = "black", pos = 3, font = 2, cex = 1) + } +} + + +#' function for specifying the color +#' +#' @param scenario TODO +#' @inheritParams plot.ensemble + +colSelect <- function(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) { + if (is.null(gcm)) { + col <- pal.scenario[which(scenarios == scenario)] + } else { + col <- if (pal == "gcms") pal.gcms[which(list_gcms() == gcm)] else pal.scenario[which(scenarios == scenario)] + } + return(col) } diff --git a/R/plot_timeSeries_input.R b/R/plot_timeSeries_input.R index 582924e..e298317 100644 --- a/R/plot_timeSeries_input.R +++ b/R/plot_timeSeries_input.R @@ -1,62 +1,60 @@ #' Input data for the time series climate change plot #' #' @description -#' Input data for the [`plot_timeSeries()`] function. Since these inputs are time-consuming to generate, +#' Input data for the [`plot_timeSeries()`] function. Since these inputs are time-consuming to generate, #' the purpose of conducting the generation of the input table in a separate function is to allow users -#' to make multiple calls to [`plot_timeSeries()`] (e.g., for comparing different climate variables) +#' to make multiple calls to [`plot_timeSeries()`] (e.g., for comparing different climate variables) #' without needing to generate the inputs each time. -#' +#' #' @details #' This function generates standardized inputs for one or multiple locations at any spatial scale. -#' If multiple locations are specified, the output is the average of the climate variables for all locations. -#' -#' Downloads of GCM time series take a long time. The `plot_timeSeries_input()` function can take >1hr -#' to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -#' we recommend users dedicate some time to run this function in background. Once the time series are cached, they -#' don't need to be downloaded again. +#' If multiple locations are specified, the output is the average of the climate variables for all locations. +#' +#' Downloads of GCM time series take some time. The `plot_timeSeries_input()` function can take ~5 minutes +#' to run for the first time it is called for a location. Once the time series are cached, they +#' don't need to be downloaded again. #' #' @template xyz #' @inheritParams downscale #' @template vars #' #' @return `data.table` of average downscaled climate variables for all locations. -#' +#' #' @examples -#' if(FALSE){ -#' # data frame of arbitrary points -#' my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) -#' -#' # generate the input data -#' my_data <- plot_timeSeries_input(my_points) -#' -#' # use the input to create a plot -#' plot_timeSeries(my_data, variable1 = "Tmin_sm") +#' if (FALSE) { +#' # data frame of arbitrary points +#' my_points <- data.frame(lon = c(-127.7300, -127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +#' +#' # generate the input data +#' my_data <- plot_timeSeries_input(my_points) +#' +#' # use the input to create a plot +#' plot_timeSeries(my_data, var1 = "Tmin_sm") #' } -#' #' +#' #' #' @export plot_timeSeries_input <- function( - xyz, + xyz, gcms = list_gcms(), ssps = list_ssps(), max_run = 10, - obs_ts_dataset = c("cru.gpcc", "climatena"), + obs_ts_dataset = c("cru.gpcc", "climatena"), obs_years = list_obs_years(), - gcm_hist_years = list_gcm_hist_years(), - gcm_ssp_years = list_gcm_ssp_years(), - vars = list_vars() -) { - data <- downscale(xyz = xyz, - gcms = gcms, - ssps = ssps, - max_run = max_run, - obs_ts_dataset = obs_ts_dataset, - obs_years = obs_years, - gcm_hist_years = gcm_hist_years, - gcm_ssp_years = gcm_ssp_years, - vars = vars + gcm_hist_years = list_gcm_hist_years(), + gcm_ssp_years = list_gcm_ssp_years(), + vars = list_vars()) { + data <- downscale( + xyz = xyz, + gcms = gcms, + ssps = ssps, + max_run = max_run, + obs_ts_dataset = obs_ts_dataset, + obs_years = obs_years, + gcm_hist_years = gcm_hist_years, + gcm_ssp_years = gcm_ssp_years, + vars = vars ) data.agg <- data[, lapply(.SD, mean), by = .(GCM, SSP, RUN, PERIOD, DATASET), .SDcols = -c("id", "GCM", "SSP", "RUN", "PERIOD", "DATASET")] return(data.agg) } - diff --git a/R/refmap.R b/R/refmap.R index 6859c6c..5dfbac8 100644 --- a/R/refmap.R +++ b/R/refmap.R @@ -1,8 +1,8 @@ #' Retrieve reference period climate maps #' @description #' This function downloads (or retrieves from cache) monthly Tmin, Tmax, and PPT climatologies (maps of long-term average climate) -#' from a specified data source for the specified bounding box. -#' It is intended for use with [`downscale_core()`], but can also be used as stand-alone raster data. +#' from a specified data source for the specified bounding box. +#' It is intended for use with [`downscale_core()`], but can also be used as stand-alone raster data. #' #' @template dbCon #' @template bbox @@ -27,7 +27,7 @@ input_refmap <- function(dbCon, bbox, reference = "refmap_climatena", cache = TRUE) { ## checks if (is(reference, "character")) { - #match.arg(reference, list_refmaps()) ## temporarily disable + # match.arg(reference, list_refmaps()) ## temporarily disable } else { if (!is(reference, "SpatRaster")) { stop( @@ -40,16 +40,25 @@ input_refmap <- function(dbCon, bbox, reference = "refmap_climatena", cache = TR if (!is(cache, "logical")) { stop("please pass a logical value to 'cache'") } - + if (!is.null(bbox)) { .check_bb(bbox) } - ## check cached - ## check cached needDownload <- TRUE + if (!grepl("normal", reference)) { + rmap_nm <- switch(reference, + refmap_prism = "normal_bc", + refmap_climr = "normal_composite", + refmap_climatena = "normal_na", + auto = "normal_composite" + ) + } else { + rmap_nm <- reference + } + - cPath <- file.path(cache_path(), "reference", reference) + cPath <- file.path(cache_path(), "reference", rmap_nm) if (dir.exists(cPath)) { bnds <- try(fread(file.path(cPath, "meta_data.csv")), silent = TRUE) @@ -83,16 +92,6 @@ input_refmap <- function(dbCon, bbox, reference = "refmap_climatena", cache = TR } if (needDownload) { - if(!grepl("normal",reference)){ - rmap_nm <- switch(reference, - refmap_prism = "normal_bc", - refmap_climr = "normal_composite", - refmap_climatena = "normal_na", - auto = "normal_composite" - ) - }else{ - rmap_nm <- reference - } message("Downloading new data...") res <- pgGetTerra(dbCon, rmap_nm, tile = TRUE, boundary = bbox, bands = 1:73) diff --git a/R/sysdata.rda b/R/sysdata.rda index ec550cc..1d8f8d8 100644 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/R/utils.R b/R/utils.R index 45c8841..24f6261 100644 --- a/R/utils.R +++ b/R/utils.R @@ -37,13 +37,13 @@ is_in_bbox <- function(newbb, oldbb) { get_bb <- function(in_xyz) { .checkXYZ(copy(in_xyz)) thebb <- c(max(in_xyz[, "lat"]), min(in_xyz[, "lat"]), max(in_xyz[, "lon"]), min(in_xyz[, "lon"])) - + if (any(is.na(thebb))) { stop("Couldn't guess bounding box. Are there NA's in 'xyz'?") } - + .check_bb(thebb) - + return(thebb) } @@ -54,15 +54,17 @@ get_bb <- function(in_xyz) { #' #' @return NULL #' @noRd -.check_bb <- function (bbox) { +.check_bb <- function(bbox) { ## check projection - minLon <- -179.0625 + minLon <- -179.0625 maxLon <- -51.5625 minLat <- 14.375 maxLat <- 83.125 - - if (any((bbox[4] < minLon), (bbox[3] > maxLon), - (bbox[2] < minLat), (bbox[1] > maxLat))) { - stop("lon and lat are not in lat/long projection EPSG:4326") + + if (any( + (bbox[4] < minLon), (bbox[3] > maxLon), + (bbox[2] < minLat), (bbox[1] > maxLat) + )) { + stop("input fields lon and lat are not in lat/long coordinates, or extent is outside ext(-179.0625, -51.5625, 14.375, 83.125)") } } diff --git a/data-raw/derivedVariables/Variables_ClimateBC.csv b/data-raw/derivedVariables/variables.csv similarity index 99% rename from data-raw/derivedVariables/Variables_ClimateBC.csv rename to data-raw/derivedVariables/variables.csv index a5b5d6d..314bef3 100644 --- a/data-raw/derivedVariables/Variables_ClimateBC.csv +++ b/data-raw/derivedVariables/variables.csv @@ -159,7 +159,7 @@ PAS_09,September precipitation as snow,PAS,precipitation as snow,09,September,Mo PAS_10,October precipitation as snow,PAS,precipitation as snow,10,October,Monthly,ratio,mm,PAS10 PAS_11,November precipitation as snow,PAS,precipitation as snow,11,November,Monthly,ratio,mm,PAS11 PAS_12,December precipitation as snow,PAS,precipitation as snow,12,December,Monthly,ratio,mm,PAS12 -PPT,Precipitation,PPT,Precipitation,,Any,Annual,ratio,mm,PPT +PPT,Precipitation,PPT,Precipitation,,Any,Annual,ratio,mm,MAP PPT_at,Autumn precipitation,PPT,precipitation,at,Autumn,Seasonal,ratio,mm,PPT_at PPT_sm,Summer precipitation,PPT,precipitation,sm,Summer,Seasonal,ratio,mm,PPT_sm PPT_sp,Spring precipitation,PPT,precipitation,sp,Spring,Seasonal,ratio,mm,PPT_sp @@ -194,7 +194,7 @@ RH_10,October relative humidity,RH,relative humidity,10,October,Monthly,interval RH_11,November relative humidity,RH,relative humidity,11,November,Monthly,interval,%,RH11 RH_12,December relative humidity,RH,relative humidity,12,December,Monthly,interval,%,RH12 SHM,Summer heat:moisture index,SHM,Summer heat:moisture index,,Annual,Annual,ratio,,SHM -Tave,Mean temperature,Tave,Mean temperature,,Any,Annual,interval,deg. C,Tave +Tave,Mean temperature,Tave,Mean temperature,,Any,Annual,interval,deg. C,MAT Tave_at,Autumn mean temperature,Tave,mean temperature,at,Autumn,Seasonal,interval,deg. C,Tave_at Tave_sm,Summer mean temperature,Tave,mean temperature,sm,Summer,Seasonal,interval,deg. C,Tave_sm Tave_sp,Spring mean temperature,Tave,mean temperature,sp,Spring,Seasonal,interval,deg. C,Tave_sp @@ -212,7 +212,7 @@ Tave_10,October mean temperature,Tave,mean temperature,10,October,Monthly,interv Tave_11,November mean temperature,Tave,mean temperature,11,November,Monthly,interval,deg. C,Tave11 Tave_12,December mean temperature,Tave,mean temperature,12,December,Monthly,interval,deg. C,Tave12 TD,Continentality (MWMT-MCMT),TD,Continentality (MWMT-MCMT),,Annual,Annual,interval,deg. C,TD -Tmax,Mean daily maximum temperature,Tmax,Mean daily maximum temperature,,Any,Annual,interval,deg. C,Tmax +Tmax,Mean daily maximum temperature,Tmax,Mean daily maximum temperature,,Any,Annual,interval,deg. C, Tmax_at,Autumn mean daily maximum temperature,Tmax,mean daily maximum temperature,at,Autumn,Seasonal,interval,deg. C,Tmax_at Tmax_sm,Summer mean daily maximum temperature,Tmax,mean daily maximum temperature,sm,Summer,Seasonal,interval,deg. C,Tmax_sm Tmax_sp,Spring mean daily maximum temperature,Tmax,mean daily maximum temperature,sp,Spring,Seasonal,interval,deg. C,Tmax_sp @@ -229,7 +229,7 @@ Tmax_09,September mean daily maximum temperature,Tmax,mean daily maximum tempera Tmax_10,October mean daily maximum temperature,Tmax,mean daily maximum temperature,10,October,Monthly,interval,deg. C,Tmax10 Tmax_11,November mean daily maximum temperature,Tmax,mean daily maximum temperature,11,November,Monthly,interval,deg. C,Tmax11 Tmax_12,December mean daily maximum temperature,Tmax,mean daily maximum temperature,12,December,Monthly,interval,deg. C,Tmax12 -Tmin,Mean daily minimum temperature,Tmin,Mean daily minimum temperature,,Any,Annual,interval,deg. C,Tmin +Tmin,Mean daily minimum temperature,Tmin,Mean daily minimum temperature,,Any,Annual,interval,deg. C, Tmin_at,Autumn mean daily minimum temperature,Tmin,mean daily minimum temperature,at,Autumn,Seasonal,interval,deg. C,Tmin_at Tmin_sm,Summer mean daily minimum temperature,Tmin,mean daily minimum temperature,sm,Summer,Seasonal,interval,deg. C,Tmin_sm Tmin_sp,Spring mean daily minimum temperature,Tmin,mean daily minimum temperature,sp,Spring,Seasonal,interval,deg. C,Tmin_sp diff --git a/data-raw/internal_data.R b/data-raw/internal_data.R index 9caca72..e4ea254 100644 --- a/data-raw/internal_data.R +++ b/data-raw/internal_data.R @@ -69,6 +69,18 @@ dbnames_ts <- data.table( ) ) -usethis::use_data(param, dbnames, dbnames_hist, dbnames_ts, dbnames_hist_obs, +dbnames_ts_fast <- data.table( + GCM = c( + "ACCESS-ESM1-5","BCC-CSM2-MR", "CanESM5", + "CNRM-ESM2-1", "EC-Earth3", "GISS-E2-1-G", "INM-CM5-0", "GFDL-ESM4", + "IPSL-CM6A-LR", "MIROC6", "MPI-ESM1-2-HR", "MRI-ESM2-0", "UKESM1-0-LL" + ), + dbname = c("access_template", "bcc_template","canesm_template","cnrm_template", "ecearth3_template","giss_template","inm_template", + "gfdl_template", "ipsl_template","miroc6_template","mpi_template","mri_template","ukesm_template"), + dbarray = c("access_array", "bcc_array","canesm_array","cnrm_array", "ecearth3_array","giss_array","inm_array", + "gfdl_array", "ipsl_array","miroc6_array","mpi_array","mri_array","ukesm_array") +) + +usethis::use_data(param, dbnames, dbnames_hist, dbnames_ts, dbnames_hist_obs, dbnames_ts_fast, overwrite = TRUE, internal = TRUE ) diff --git a/data-raw/namingChanges.csv b/data-raw/namingChanges.csv index c0223a8..6168e84 100644 --- a/data-raw/namingChanges.csv +++ b/data-raw/namingChanges.csv @@ -51,3 +51,5 @@ function name,,gcm_hist_input(),input_gcm_hist() function name,,gcm_ts_input(),input_gcm_ssp() function name,,historic_input(),input_obs() function name,,historic_input_ts(),input_obs_ts() +variable,,DD_0,DDsub0 +variable,,DD_18,DDsub18 diff --git a/data-raw/variables.R b/data-raw/variables.R index 32a67fa..25858c7 100644 --- a/data-raw/variables.R +++ b/data-raw/variables.R @@ -2,6 +2,6 @@ library(data.table) -variables <- fread("data-raw/derivedVariables/Variables_ClimateBC.csv") +variables <- fread("data-raw/derivedVariables/variables.csv") usethis::use_data(variables, overwrite = TRUE, internal = FALSE) diff --git a/data/name_changes.rda b/data/name_changes.rda index fa50984..2dd69c0 100644 Binary files a/data/name_changes.rda and b/data/name_changes.rda differ diff --git a/data/variables.rda b/data/variables.rda index 49d391c..64bcb7e 100644 Binary files a/data/variables.rda and b/data/variables.rda differ diff --git a/data_processing/CreateDB.R b/data_processing/CreateDB.R index da51749..a2692bd 100644 --- a/data_processing/CreateDB.R +++ b/data_processing/CreateDB.R @@ -7,31 +7,31 @@ library(analogsea) library(ssh) -bc_rast <- rast("../Common_Files/composite_wna_wlrdem.tif")[[32]] -plot(bc_rast) -bc_rast[!is.na(bc_rast)] <- 1L -plot(bc_rast) -bc_rast2 <- terra::aggregate(bc_rast, fact = 4, fun = "max") -plot(bc_rast2) - -writeRaster(bc_rast2, filename = "inst/extdata/wna_outline.tif", datatype = "INT1U") +# bc_rast <- rast("../Common_Files/composite_wna_wlrdem.tif")[[32]] +# plot(bc_rast) +# bc_rast[!is.na(bc_rast)] <- 1L +# plot(bc_rast) +# bc_rast2 <- terra::aggregate(bc_rast, fact = 4, fun = "max") +# plot(bc_rast2) +# +# writeRaster(bc_rast2, filename = "inst/extdata/wna_outline.tif", datatype = "INT1U") ##process NA normals # Load normal files -dir_normal <- "../Common_Files/colin_climatology/" +dir_normal <- "../Common_Files/climatena_normals/Normal_1961_1990MP/" files <- list.files(dir_normal, full.names = T) -fnames <- list.files(dir_normal) +files <- files[grep("Tmin|Tmax|PPT",files)] d <- rast(files) -names(d) <- c("PPT01", "PPT02", "PPT03", "PPT04", "PPT05", "PPT06", "PPT07", "PPT08", "PPT09", "PPT10", - "PPT11", "PPT12", "Tmax01", "Tmax02", "Tmax03", "Tmax04", "Tmax05", "Tmax06", "Tmax07", - "Tmax08", "Tmax09", "Tmax10", "Tmax11", "Tmax12", "Tmin01", "Tmin02", "Tmin03", "Tmin04", - "Tmin05", "Tmin06", "Tmin07", "Tmin08", "Tmin09", "Tmin10", "Tmin11", "Tmin12") +# names(d) <- c("PPT01", "PPT02", "PPT03", "PPT04", "PPT05", "PPT06", "PPT07", "PPT08", "PPT09", "PPT10", +# "PPT11", "PPT12", "Tmax01", "Tmax02", "Tmax03", "Tmax04", "Tmax05", "Tmax06", "Tmax07", +# "Tmax08", "Tmax09", "Tmax10", "Tmax11", "Tmax12", "Tmin01", "Tmin02", "Tmin03", "Tmin04", +# "Tmin05", "Tmin06", "Tmin07", "Tmin08", "Tmin09", "Tmin10", "Tmin11", "Tmin12") -dem <- rast("../Common_Files/composite_WNA_dem.tif") +dem <- rast("../climatena/InputFiles/na4000.asc") #r <- r[[grep("PPT|Tmin|Tmax", names(r))]] lr <- lapse_rate( - normal = d, + reference = d, dem = dem, NA_replace = TRUE, nthread = 4, @@ -42,7 +42,7 @@ names(lr) <- paste0("lr_",names(lr)) # Actual writing terra::writeRaster( c(d, lr, dem), - file.path("../Common_Files/composite_wna_wlrdem.tif"), + file.path("../Common_Files/climatena_wlrdem.tif"), overwrite = TRUE, gdal="COMPRESS=NONE" ) diff --git a/data_processing/Test_Script.R b/data_processing/Test_Script.R index d25d2cf..cd1ab97 100644 --- a/data_processing/Test_Script.R +++ b/data_processing/Test_Script.R @@ -1,20 +1,60 @@ library(data.table) + library(terra) library(climr) + +wlr <- rast("../Common_Files/composite_WNA_dem.tif") +plot(wlr) + +dem <- rast("../Common_Files/dem_noram_lowres.tif") +test <- rast("../Common_Files/climatena_normals/Normal_1961_1990MP/Tmin07.asc") +plot(test) + +my_grid <- as.data.frame(dem, cells = TRUE, xy = TRUE) +colnames(my_grid) <- c("id", "lon", "lat", "elev") # rename column names to what climr expects +db <- data_connect() +bbox <- get_bb(my_grid) +bbox2 <- c(20,14.83,-80,-120) +refmap <- input_refmap(db, bbox) + +plot(refmap$Tmax_07) + pts <- data.frame(lon = c(-124.11, -125.11), lat = rep(48.82, 2), elev = rep(25,2), id = 1:2) bbox <- get_bb(pts[2,]) dbcon <- data_connect() -test <- normal_input(dbcon, bbox) +test <- input_refmap(dbcon, bbox) plot(test[[8]]) +dem <- rast("../Common_Files/dem_noram_lowres.tif") +test <- rast("../Common_Files/climatena_normals/Normal_1961_1990MP/Tmin07.asc") +plot(test) + +my_grid <- as.data.frame(dem, cells = TRUE, xy = TRUE) +colnames(my_grid) <- c("id", "lon", "lat", "elev") # rename column names to what climr expects +climr <- downscale( + xyz = my_grid, which_refmap = "refmap_climatena", gcms = list_gcms()[1], ssps = list_ssps()[1:3], + gcm_periods = list_gcm_periods(), run_nm = "r1i1p1f1", vars = "MAT" +) + +X <- rast(dem) +X[climr[, id]] <- climr$MAT +plot(X) + + +db <- data_connect() +bbox <- get_bb(my_grid) +bbox2 <- c(20,14.83,-80,-120) +refmap <- input_refmap(db, bbox) + + projected <- climr_downscale(pts[2,], - gcm_models = list_gcm()[c(4)], - ssp = list_ssp()[c(1,2)], - max_run = 3, - gcm_hist_years = 1851:2014, - gcm_ts_years = 2015:2100 + gcm_models = list_gcm()[c(4)], + ssp = list_ssp()[c(1,2)], + max_run = 3, + gcm_hist_years = 1851:2014, + gcm_ts_years = 2015:2100 ) dat <- fread("../climatena/Perioddat/Year_1905.ann") @@ -58,12 +98,18 @@ data <- downscale(xyz = pt, vars = list_vars() ) +"#E404D1" plot_timeSeries(data, var1 = "Tmax_08") -data <- plot_timeSeries_input(pt, gcms = list_gcms()[4], vars = "MAP") -plot_timeSeries(data, var1 = "MAP") +library(abind) +library(dplyr) +library(terra) +data <- plot_timeSeries_input(pt, gcms = list_gcms()[1], max_run = 5) +data <- downscale(pt, gcms = list_gcms()[1], gcm_ssp_years = list_gcm_ssp_years(), + ssps = list_ssps()[1:3], max_run = 5L) +plot_timeSeries(data, var1 = "MAT") data <- downscale(xyz = pt, - gcm_models = list_gcms()[2], + gcm_models = list_gcms()[5], ssp = list_ssps(), max_run = 10, historic_ts_dataset = c("cru.gpcc", "climatena"), diff --git a/data_processing/timeseries_speedup.R b/data_processing/timeseries_speedup.R new file mode 100644 index 0000000..549c4d2 --- /dev/null +++ b/data_processing/timeseries_speedup.R @@ -0,0 +1,224 @@ +library(terra) +library(data.table) +library(climr) +library(stringi) +library(DBI) +library(dplyr) +library(pool) +library(RPostgres) +conn <- dbPool(RPostgres::Postgres(), + dbname = 'climr', + host = '146.190.244.244', + port = 5432, + user = 'postgres', + password = 'climr2022') + +pts <- data.frame(lon = c(-124.11, -125.11), lat = rep(48.82, 2), elev = rep(25,2), id = 1:2) + +bbox <- get_bb(pts) +dbcon <- data_connect() +library(abind) +library(dplyr) +res <- input_gcm_ssp(dbcon, bbox, gcms = list_gcms()[2], years = 2020:2050, + max_run = 5L, cache = T, fast = T) + +input_gcm_ssp_fast <- function(dbCon, bbox = NULL, gcm_nm = list_gcms(), ssps = list_ssps(), + years = 2020:2030, max_run = 0L, cache = TRUE){ + gcmcode <- "gfdl_template" + gcmarray <- "gfdl_array" + ssps <- c("ssp245", "ssp370") + years = 2020:2030 + + template <- pgGetTerra(dbCon, name = gcmcode, tile = F, bands = 1, boundary = bbox) + + if(length(years) >= 79){ ##faster if almost all years are selected + results <- tbl(dbCon, sql(paste0("select cellid, ssp, year, run, vals from ",gcmarray," where cellid in (",paste0(values(template)[,1], collapse = ','),") + and ssp in ('",paste(ssps, collapse = "','"),"')"))) + }else{ + results <- tbl(dbCon, sql(paste0("select cellid, ssp, year, run, vals from ",gcmarray," where cellid in (",paste0(values(template)[,1], collapse = ','),") + and year in ('",paste(years, collapse = "','"),"') and ssp in ('",paste(ssps, collapse = "','"),"')"))) + } + + dat <- + results %>% + mutate(vals = unnest(vals)) %>% + collect() + setDT(dat) + setorder(dat, cellid, year, ssp, run) + dat[, month := rep(sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + nrow(dat)/(12*3))] + + dat[, fullnm := paste(month, ssp, run, year, sep = "_")] + cell_nums = cells(template) + coords <- rowColFromCell(template, cell_nums) + cellcoords <- data.table(coords, values(template)) + setnames(cellcoords, c("row","col","cellid")) + dat[cellcoords, `:=`(row = i.row, col = i.col), on = "cellid"] + temp <- dat[,.(row,col,fullnm,vals)] + t2 <- dcast(temp, fullnm + row ~ col, value.var = "vals") + + t_array = split(as.data.frame(t2[,!c("fullnm","row")]), t2$fullnm) + t3 <- abind(t_array, along = 3) + t_rast <- rast(t3) + ext(t_rast) <- ext(template) + names(t_rast) <- names(t_array) + return(t_rast) +} + +bbox <- get_bb(pts) +dbcon <- data_connect() + +res <- input_gcm_ssp_fast(dbCon = dbcon, bbox, ssps = ssps, years = 2020:2050) + +template <- pgGetTerra(conn, name = "gfdl_template", tile = F, bands = 1, boundary = bbox) +#plot(template) +results <- tbl(conn, sql(paste0("select cellid, ssp, year, run, vals from gfdl_array where cellid in (",paste0(values(template)[,1], collapse = ','),") order by cellid, year, ssp, run"))) +results <- tbl(conn, sql(paste0("select cellid, ssp, year, run, vals from gfdl_array where cellid in (",paste0(values(template)[,1], collapse = ','),")"))) + +approach_2 <- + results %>% + mutate(vals = unnest(vals)) %>% + collect() + +setDT(approach_2) +setorder(approach_2, cellid, year, ssp, run) +approach_2[, month := rep(sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + nrow(approach_1)/(12*3))] + +approach_1$month = rep(sort(sprintf(c("PPT_%02d", "Tmax_%02d", "Tmin_%02d"), sort(rep(1:12, 3)))), + nrow(approach_1)/(12*3)) +setDT(approach_1) + +approach_1[, fullnm := paste(month, ssp, run, year, sep = "_")] +cell_nums = cells(template) +coords <- rowColFromCell(template, cell_nums) +cellcoords <- data.table(coords, values(template)) +setnames(cellcoords, c("row","col","cellid")) +approach_1[cellcoords, `:=`(row = i.row, col = i.col), on = "cellid"] +temp <- approach_1[,.(row,col,fullnm,vals)] +t2 <- dcast(temp, fullnm + row ~ col, value.var = "vals") +library(abind) +t2[,fullnm := as.numeric(as.factor(fullnm))] +t_array = split(as.data.frame(t2[,!c("fullnm","row")]), t2$fullnm) +t3 <- abind(t_array, along = 3) +t_rds <- sds(x = t3) +t_rast <- rast(t3) +ext(t_rast) <- ext(template) +names(t_rast) <- names(t_array) +plot(t_rast) + +library(str2str) +library(torch) +Sys.setenv(CUDA=FALSE) +torch_tensor(1) +t_array <- d2a(temp, dim.nm = c("row","col","fullnm")) +t_array <- array(data = temp$vals, + dim=c(length(unique(df$x)), + length(unique(df$y)), + length(unique(df$i))), + dimnames=list(unique(df$x), unique(df$y), unique(df$i)) +) + +t2 <- dcast(approach_1, fullnm ~ cellid, value.var = "vals") + +temp <- copy(template) +cellids <- as.numeric(names(t2)[-1]) +t3 <- as.matrix(t2[,!"fullnm"]) +for(i in 1:nrow(t2)){ ##too slow! + values(template) <- t3[i,] + if(i == 1){ + ts_out <- copy(template) + }else{ + add(ts_out) <- template + } +} + +##need x and y - use getXYfromCell, convert to array +gcm_nm <- "GFDL-ESM4" +template_nm <- "access_template.tif" +table_nm <- "access_ts_tbl" + +create_ts_table <- function(gcm_nm, template_nm, tmp_table_nm, table_nm){ + fnames <- list.files("../Common_Files/gcmts_deltas/", pattern = gcm_nm, full.names = T) + for(fname in fnames){ + cat("Processing",fname,"\n") + orig <- rast(fname) + if(fname == fnames[[1]]){ + template <- orig[[1]] + temp <- 1:ncell(template) + values(template) <- temp + writeRaster(template, paste0("../Common_Files/gcmts_templates/",template_nm), overwrite = TRUE) + } + nlays <- nlyr(orig) + splits <- c(seq(1,nlays,by = 10000),nlays+1) + for(i in 1:(length(splits)-1)){ + cat("Start:",splits[i],"Stop:",(splits[i+1]-1),"\n") + t2 <- values(orig[[splits[i]:(splits[i+1]-1)]], mat = TRUE) ##need to process in chunks + t2 <- data.table(t2) + t2[,cellid := 1:nrow(t2)] + t3 <- melt(t2, id.vars = "cellid") + t3[,c("model","var","month","ssp","run","year") := transpose(stri_split_fixed(variable,"_"))] + t3[,c("variable","model") := NULL] + dbWriteTable(conn, tmp_table_nm, t3, row.names = F, append = TRUE) + } + } + q <- + paste0("create table ",table_nm," as ( + select cellid, ssp, run, year, array_agg(value) as vals + from ",tmp_table_nm," + group by cellid, ssp, run, year + )") + cat("array aggregate\n") + dbExecute(conn, q) + +} + +create_ts_table("CNRM-ESM2-1","cnrm_template.tif","cnrm_temp","cnrm_array") +q <- paste0("create table gfdl_array as ( + select cellid, ssp, run, year, array_agg(value) as vals + from (select * from gfdl_temp order by cellid, year, ssp, run, var, month) as a + group by cellid, ssp, run, year + )") + +dbExecute(conn, q) +create_ts_table("GFDL-ESM4","gfdl_template.tif","gfdl_temp","gfdl_array") + +##testing +my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +ts_dat <- plot_timeSeries_input(my_points, gcms = "GFDL-ESM4") +plot_timeSeries(ts_dat) + + +orig <- rast("../Common_Files/gcmts_deltas/gcmts.tmin.ACCESS-ESM1-5.deltas.tif") +template <- orig[[1]] +temp <- 1:ncell(template) +values(template) <- temp +plot(template) + +writeRaster(template, "access_template.tif") + + +dat2 <- dbGetQuery(conn, "select cellid, ssp, year, vals from access_array where cellid in (1998,1999,2000,2001,2002,2003)") +setDT(dat2) +results <- tbl(conn, sql("select cellid, ssp, year, run, vals from access_array where cellid in (1998,1999,2000,2001,2002,2003)")) +approach_1 <- + results %>% + mutate(vals = unnest(vals)) %>% + collect() + +s <- orig[[1:12]] +values(s) <- NA + +setDT(approach_1) +tmp <- approach_1[ssp == "ssp370" & year == "2050" & run == "ensembleMean",] +tmp[, month := rep(1:12, nrow(tmp)/12)] +t2 <- dcast(tmp, cellid ~ month, value.var = "vals") +set.values(s, t2$cellid, as.matrix(t2[,!"cellid"]), layer = 1:12) +plot(s) +##rows are cellids, cols are layer values + +library(terra) +s <- rast(ncols=5, nrows=5, nlyrs=3) +set.values(s, 1:10, runif(10), layer=2) + +set.values(s, 11:20, cbind(1, runif(10)), layer=2:3) diff --git a/man-roxygen/run_nm.R b/man-roxygen/run_nm.R new file mode 100644 index 0000000..a728240 --- /dev/null +++ b/man-roxygen/run_nm.R @@ -0,0 +1 @@ +#' @param run_nm character. `NULL` or length >= 1. Name of specified run(s) to return, instead of using `max_run`. Use the `list_runs_*()` functions to list available runs.Defaults to `NULL`. \ No newline at end of file diff --git a/man/climr-package.Rd b/man/climr-package.Rd index 17668bd..f08ee31 100644 --- a/man/climr-package.Rd +++ b/man/climr-package.Rd @@ -6,7 +6,7 @@ \alias{climr-package} \title{climr: Downscaling climate data in R} \description{ -`climr` is an R package that builds on the downscaling concepts operationalized in the ClimateNA tool (climatena.ca) (Wang et al. 2016). It provides downscaling of observational and simulated climate data using change-factor downscaling, a simple method that adds low-spatial-resolution climate anomalies to a high-spatial-resolution reference climatological map, with additional elevation adjustment for temperature. Elevation-adjusted monthly values of basic climate elements (temperature and precipitation) are then used to estimate derived variables (e.g., degree-days) based on published equations and parameters from Wang et al. 2016. `climr` is designed to be fast and to minimize local data storage requirements. To do so, it uses a remote PostGIS database, and optionally caches data locally. +Builds on the downscaling concepts operationalized in the 'ClimateNA' tool (\url{https://climatena.ca}) (Wang et al. 2016 \doi{10.1371/journal.pone.0156720}). It provides downscaling of observational and simulated climate data using change-factor downscaling, a simple method that adds low-spatial-resolution climate anomalies to a high-spatial-resolution reference climatological map, with additional elevation adjustment for temperature. Elevation-adjusted monthly values of basic climate elements (temperature and precipitation) are then used to estimate derived variables (e.g., degree-days) based on published equations and parameters from Wang et al. 2016. This package is designed to be fast and to minimize local data storage requirements. To do so, it uses a remote 'PostGIS' database, and optionally caches data locally. } \author{ \strong{Maintainer}: Kiri Daust \email{kiri.daust@gov.bc.ca} diff --git a/man/colSelect.Rd b/man/colSelect.Rd new file mode 100644 index 0000000..a61abbb --- /dev/null +++ b/man/colSelect.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_timeSeries.R +\name{colSelect} +\alias{colSelect} +\title{function for specifying the color} +\usage{ +colSelect(scenario, gcm, pal.scenario, scenarios, pal, pal.gcms) +} +\arguments{ +\item{scenario}{TODO} + +\item{gcm}{TODO} + +\item{pal.scenario}{TODO} + +\item{scenarios}{TODO} + +\item{pal}{character. color palette. Options are "scenario", for use when comparing scenarios, +and "gcms", for use when comparing GCMs.} + +\item{pal.gcms}{TODO} +} +\description{ +function for specifying the color +} diff --git a/man/data-option-lists.Rd b/man/data-option-lists.Rd index a4d8169..6e9581b 100644 --- a/man/data-option-lists.Rd +++ b/man/data-option-lists.Rd @@ -4,14 +4,15 @@ \alias{list_gcms} \alias{list_ssps} \alias{list_gcm_periods} -\alias{list_run} +\alias{list_runs_ssp} +\alias{list_runs_historic} \alias{list_refmaps} \alias{list_obs_periods} \alias{list_vars} \alias{list_obs_years} \alias{list_gcm_ssp_years} \alias{list_gcm_hist_years} -\title{List available runs, global circulation models, periods and climate scenarios} +\title{List available runs, global climate models (GCMs), time periods and scenarios (SSPs)} \usage{ list_gcms() @@ -19,7 +20,9 @@ list_ssps() list_gcm_periods() -list_run(dbCon, gcms) +list_runs_ssp(gcm, ssp) + +list_runs_historic(gcm) list_refmaps() @@ -34,9 +37,9 @@ list_gcm_ssp_years() list_gcm_hist_years() } \arguments{ -\item{dbCon}{A db connection object created by \code{data_connect}.} +\item{gcm}{Name of GCM} -\item{gcms}{Character vector to specify requested GCMs} +\item{ssp}{Name of scenario Must be one of the elements in list_ssps().} \item{set}{character. One of All, Monthly, Seasonal, Annual, or any combination thereof. Defaults to "All".} @@ -46,25 +49,27 @@ list_gcm_hist_years() a character vector. } \description{ -\code{list_gcms} lists available global circulation models. +\code{list_gcms} lists available global climate models. + +\code{list_ssps} lists available greenhouse gas concentration scenarios (SSP-RCPs). -\code{list_ssps} lists available shared socioeconomic pathways. +\code{list_gcm_periods} lists available 20-year normal periods for GCM simulations. -\code{list_gcm_periods} lists available periods. +lists available runs for a given GCM/ssp. -\code{list_run} lists available runs for a given GCM. +lists available runs from the historical simulation (1851-2014) for a specified GCM. \code{list_refmaps} lists available reference maps of gridded climate normals -\code{list_obs_periods} lists available observational periods +\code{list_obs_periods} lists available normal periods for observational climate data -\code{list_vars} lists climate variables +\code{list_vars} lists available climate variables -\code{list_obs_years} lists available years for obs time series +\code{list_obs_years} lists available years for time series of observational climate data -\code{list_gcm_ssp_years} lists available years for future projections' time series +\code{list_gcm_ssp_years} lists available years for time series of global climate model future simulations -\code{list_gcm_hist_years} lists available years for obs projections' time series +\code{list_gcm_hist_years} lists available years for time series of global climate model historical simulations } \details{ Currently available reference maps of gridded climate normals (\code{list_refmaps()}) are: diff --git a/man/downscale.Rd b/man/downscale.Rd index 6b27042..15b265c 100644 --- a/man/downscale.Rd +++ b/man/downscale.Rd @@ -16,6 +16,7 @@ downscale( gcm_ssp_years = NULL, gcm_hist_years = NULL, max_run = 0L, + run_nm = NULL, cache = TRUE, ... ) @@ -64,6 +65,8 @@ Defaults to \code{NULL}.} A value of 0 returns the \code{ensembleMean} only. Runs are included in the order they are found in the models data until \code{max_run} is reached. Defaults to 0L.} +\item{run_nm}{character. \code{NULL} or length >= 1. Name of specified run(s) to return, instead of using \code{max_run}. Use the \verb{list_runs_*()} functions to list available runs.Defaults to \code{NULL}.} + \item{cache}{logical. Cache data locally? Default \code{TRUE}} \item{...}{other arguments passed to \code{\link[=downscale_core]{downscale_core()}}. Namely: \code{return_refperiod}, @@ -110,7 +113,7 @@ in_xyz <- data.frame( ## historic observational time series vars <- c("PPT", "CMD", "Tave_07") climate_norms_hist <- downscale( - xyz = in_xyz, + xyz = in_xyz, which_refmap = "auto", return_refperiod = TRUE, obs_periods = "2001_2020", @@ -129,7 +132,7 @@ climate_norms_hist <- downscale( ## future projections for annual variables from three models climate_norms_fut <- downscale( xyz = in_xyz, which_refmap = "auto", - gcms = list_gcms()[c(1,5,6)], + gcms = list_gcms()[c(1, 5, 6)], ssps = list_ssps()[2], gcm_periods = list_gcm_periods()[1:2], # gcm_ssp_years = 2020:2060, diff --git a/man/downscale_core.Rd b/man/downscale_core.Rd index 236c396..3c906ad 100644 --- a/man/downscale_core.Rd +++ b/man/downscale_core.Rd @@ -78,7 +78,7 @@ in the example for this function. We recommend \code{\link[=downscale]{downscale()}} for most purposes. } \examples{ -## +## library(terra) xyz <- data.frame(lon = runif(10, -130, -106), lat = runif(10, 37, 50), elev = runif(10), id = 1:10) @@ -88,10 +88,10 @@ thebb <- get_bb(xyz) ## get database connection dbCon <- data_connect() -# obtain the climatena 1961-1990 normals for the study area. +# obtain the climatena 1961-1990 normals for the study area. refmap <- input_refmap(dbCon, thebb, reference = "refmap_climatena") -# obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. +# obtain the low-resolution climate data for a single gcm, 20-year period, and ssp scenario. gcm_raw <- input_gcms(dbCon, thebb, list_gcms()[3], list_ssps()[1], period = list_gcm_periods()[2]) # downscale the GCM data @@ -100,8 +100,10 @@ gcm_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = gcm_raw, var # create an input of uniform warming of 2 degrees Celsius and no precipitation change, for use as a null comparison to the GCM warming null <- gcm_raw #' use the gcm input object as a template names(null) <- "null_2C" -names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep="_")) -for(var in names(null[[1]])){ values(null[[1]][[var]]) <- if(length(grep("PPT", var)==1)) 1 else 2 } #' repopulate with the null values +names(null[[1]]) <- sapply(strsplit(names(null[[1]]), "_"), function(x) paste("null2C", x[2], x[3], "NA", "NA", "NA", "NA", sep = "_")) +for (var in names(null[[1]])) { + values(null[[1]][[var]]) <- if (length(grep("PPT", var) == 1)) 1 else 2 +} #' repopulate with the null values # downscale the null values for variables of interest null_downscaled <- downscale_core(xyz = xyz, refmap = refmap, gcms = null, vars = c("MAT", "PAS")) diff --git a/man/gcms-input-data.Rd b/man/gcms-input-data.Rd index 9a936cc..1de1469 100644 --- a/man/gcms-input-data.Rd +++ b/man/gcms-input-data.Rd @@ -13,7 +13,8 @@ input_gcms( ssps = list_ssps(), period = list_gcm_periods(), max_run = 0L, - cache = TRUE + cache = TRUE, + run_nm = NULL ) input_gcm_hist( @@ -22,7 +23,8 @@ input_gcm_hist( gcms = list_gcms(), years = 1901:2014, max_run = 0L, - cache = TRUE + cache = TRUE, + run_nm = NULL ) input_gcm_ssp( @@ -32,7 +34,9 @@ input_gcm_ssp( ssps = list_ssps(), years = 2020:2030, max_run = 0L, - cache = TRUE + cache = TRUE, + run_nm = NULL, + fast = TRUE ) } \arguments{ @@ -55,8 +59,12 @@ are found in the models data until \code{max_run} is reached. Defaults to 0L.} \item{cache}{logical. Specifying whether to cache new data locally or no. Defaults to \code{TRUE}.} +\item{run_nm}{character. \code{NULL} or length >= 1. Name of specified run(s) to return, instead of using \code{max_run}. Use the \verb{list_runs_*()} functions to list available runs.Defaults to \code{NULL}.} + \item{years}{Numeric or character vector in \code{2020:2100}. Defaults to \code{2020:2030}. See \code{\link[=list_gcm_ssp_years]{list_gcm_ssp_years()}} for available years.} + +\item{fast}{Logical. Should we use the faster method of downloading data from the database using arrays instead of Postgis rasters?} } \value{ A \code{list} of \code{SpatRasters}, each with possibly multiple layers, that can diff --git a/man/plot.ensemble.Rd b/man/plot.ensemble.Rd new file mode 100644 index 0000000..132dab1 --- /dev/null +++ b/man/plot.ensemble.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_timeSeries.R +\name{plot.ensemble} +\alias{plot.ensemble} +\title{Function for plotting time series for gcms or compiled ensemble} +\usage{ +plot.ensemble( + x, + var, + scenarios.selected, + scenarios, + showrange = TRUE, + simplify = TRUE, + gcm = NULL, + pal, + pal.scenario, + pal.gcms, + refline = FALSE, + showmean = TRUE, + endlabel = "change", + element, + compile = TRUE, + var2 = NULL, + element1, + element2, + yeartime.names, + yeartimes, + yeartime +) +} +\arguments{ +\item{x}{climate data} + +\item{var}{TODO} + +\item{scenarios.selected}{TODO} + +\item{scenarios}{TODO} + +\item{showrange}{logical. Plot a shaded region indicating the minimum and maximum of the +selected ensemble of GCM simulations for each selected scenario.} + +\item{simplify}{logical. Simplify the ensemble range and mean using a smoothing spline.} + +\item{gcm}{TODO} + +\item{pal}{character. color palette. Options are "scenario", for use when comparing scenarios, +and "gcms", for use when comparing GCMs.} + +\item{pal.scenario}{TODO} + +\item{pal.gcms}{TODO} + +\item{refline}{logical. Plot the 1961-1990 reference period mean for the selected var +and extend this line to the year 2100 as a visual reference.} + +\item{showmean}{logical. Plot the ensemble mean time series. Multi-model ensemble means are +calculated from the mean of simulations for each model.} + +\item{endlabel}{character. Add a label to the end of each simulated time series. Options +are "change", to indicate the change in year 2100 relative to the 1961-1990 baseline, or "gcms" +to indicate the global climate model.} + +\item{element}{TODO} + +\item{compile}{logical. Compile multiple global climate models into a multi-model ensemble. +If \code{FALSE} the single-model ensembles are plotted individually.} + +\item{var2}{character. A second climate var to plot in combination with \code{var1}. +options are \code{\link[=list_vars]{list_vars()}}.} + +\item{element1}{TODO} + +\item{element2}{TODO} + +\item{yeartime.names}{TODO} + +\item{yeartimes}{TODO} + +\item{yeartime}{TODO} +} +\description{ +Function for plotting time series for gcms or compiled ensemble +} diff --git a/man/plot_timeSeries.Rd b/man/plot_timeSeries.Rd index 83b4d13..26a02cd 100644 --- a/man/plot_timeSeries.Rd +++ b/man/plot_timeSeries.Rd @@ -91,8 +91,7 @@ to indicate the global climate model.} \item{yearlines}{logical. Add vertical lines on every fifth year as a visual reference} -\item{legend_pos}{character. Position of the legend. Viable options are c("bottomright", -"bottomleft", "topleft", and "topright").} +\item{legend_pos}{character. Position of the legend. Viable options are \code{c("bottomright", "bottomleft", "topleft", "topright")}.} } \value{ NULL. Draws a plot in the active graphics device. @@ -117,44 +116,42 @@ Some combinations of \code{var1} and \code{var2} are not compatible or meaningfu Examples of meaningful combinations are winter vs summer values of the same climate var or minimum vs. maximum temperatures. -Downloads of GCM time series take a long time. The \code{plot_timeSeries_input()} function can take >1hr -to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -we recommend users dedicate some time to run this function in background. Once the time series are cached, they +Downloads of GCM time series take some time. The \code{plot_timeSeries_input()} function can take ~5 minutes +to run for the first time it is called for a location. Once the time series are cached, they don't need to be downloaded again. } \examples{ -if(FALSE){ -# data frame of arbitrary points -my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +if (FALSE) { + # data frame of arbitrary points + my_points <- data.frame(lon = c(-127.7300, -127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) -# generate the input data -my_data <- plot_timeSeries_input(my_points) + # generate the input data + my_data <- plot_timeSeries_input(my_points) -# use the input to create a plot -plot_timeSeries(my_data, var1 = "Tmin_sm") + # use the input to create a plot + plot_timeSeries(my_data, var1 = "Tmin_sm") -# compare observational time series -plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) + # compare observational time series + plot_timeSeries(my_data, var1 = "Tmin_sm", obs_ts_dataset = c("cru.gpcc", "climatena")) -# compare mean daily minimum and maximum temperatures -plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") + # compare mean daily minimum and maximum temperatures + plot_timeSeries(my_data, var1 = "Tmin_sm", var2 = "Tmax_sm") -# compare summer and winter temperatures (without simplifying the ensemble range) -plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) + # compare summer and winter temperatures (without simplifying the ensemble range) + plot_timeSeries(my_data, var1 = "Tmax_sm", var2 = "Tmax_wt", simplify = FALSE) -# compare global climate models -plot_timeSeries(my_data, gcms = list_gcms()[c(7,13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar=c(3,3,0.1,6), showObserved = FALSE) + # compare global climate models + plot_timeSeries(my_data, gcms = list_gcms()[c(7, 13)], pal = "gcms", ssps = list_ssps()[2], showmean = FALSE, compile = FALSE, simplify = FALSE, endlabel = "gcms", mar = c(3, 3, 0.1, 6), showObserved = FALSE) -# export plot to a temporary directory, including a title -figDir <- tempdir() -png( - filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", - width = 6, height = 5, pointsize = 10, res = 300 -) -plot_timeSeries(my_data, var1 = "Tmin_sm", mar=c(3,3,2,4)) -title("Historical and projected summer night-time warming in the Bulkley Valley, BC") -dev.off() + # export plot to a temporary directory, including a title + figDir <- tempdir() + png( + filename = file.path(figDir, "plot_test.png"), type = "cairo", units = "in", + width = 6, height = 5, pointsize = 10, res = 300 + ) + plot_timeSeries(my_data, var1 = "Tmin_sm", mar = c(3, 3, 2, 4)) + title("Historical and projected summer night-time warming in the Bulkley Valley, BC") + dev.off() } - } diff --git a/man/plot_timeSeries_input.Rd b/man/plot_timeSeries_input.Rd index 5af98aa..06fd831 100644 --- a/man/plot_timeSeries_input.Rd +++ b/man/plot_timeSeries_input.Rd @@ -64,21 +64,20 @@ without needing to generate the inputs each time. This function generates standardized inputs for one or multiple locations at any spatial scale. If multiple locations are specified, the output is the average of the climate variables for all locations. -Downloads of GCM time series take a long time. The \code{plot_timeSeries_input()} function can take >1hr -to run for the first time it is called for a location. We are looking into ways to speed this up, but until then -we recommend users dedicate some time to run this function in background. Once the time series are cached, they +Downloads of GCM time series take some time. The \code{plot_timeSeries_input()} function can take ~5 minutes +to run for the first time it is called for a location. Once the time series are cached, they don't need to be downloaded again. } \examples{ -if(FALSE){ -# data frame of arbitrary points -my_points <- data.frame(lon = c(-127.7300,-127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) +if (FALSE) { + # data frame of arbitrary points + my_points <- data.frame(lon = c(-127.7300, -127.7500), lat = c(55.34114, 55.25), elev = c(711, 500), id = 1:2) -# generate the input data -my_data <- plot_timeSeries_input(my_points) + # generate the input data + my_data <- plot_timeSeries_input(my_points) -# use the input to create a plot -plot_timeSeries(my_data, variable1 = "Tmin_sm") + # use the input to create a plot + plot_timeSeries(my_data, var1 = "Tmin_sm") } -#' +#' } diff --git a/tests/manual_tests/EDA_Eref.R b/tests/manual_tests/EDA_Eref.R new file mode 100644 index 0000000..b259700 --- /dev/null +++ b/tests/manual_tests/EDA_Eref.R @@ -0,0 +1,158 @@ +## messy code used to troubleshoot the Eref modification required to the Hargreaves equation. +## colin mahony + + +day_month <- c(31, 28.25, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) +day_julian <- c(15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349) + +calc_Eref <- function(m, tmmin, tmmax, latitude) { + Eref <- numeric(length(tmmax)) + tmean <- (tmmax + tmmin) / 2 + i <- which(tmean >= 0) + day_month <- c(31, 28.25, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) + day_julian <- c(15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349) + # Paper unclear, 1.18 - 0.0065 in Appendix, 1.18 - 0.0067 in paper + # Wangetal2012_ClimateWNA_JAMC-D-11-043.pdf + # Probably missing + Eref[i] <- 0.0023 * day_month[m] * + calc_S0_I(day_julian[m], tmean[i], latitude[i]) * + (tmean[i] + 17.8) * sqrt(tmmax[i] - tmmin[i]) * + (1.18 - 0.0065 * latitude[i]) + + Eref[is.na(tmmax)] <- tmmax[is.na(tmmax)] ## use tmmax[is.na(tmmax)] to respect NA type + + return(Eref) +} + +# Eref by latitude and month at Tmean=13C +x <- seq(40, 80, 0.1) +plot(x, rep(0, length(x)), ylim=c(0,80), col="white", xlab="Latitude", ylab="Eref", main="Eref by latitude and month at Tmean=13C") +latitude=x +tmmin=rep(11, length(x)) +tmmax=rep(15, length(x)) +for(m in 1:12){ + y <- calc_Eref(m, tmmin, tmmax, latitude) + points(x,y, col=m, pch=m, cex=0.1) + text(40, y[1], m, col=m) + text(80, y[length(x)], m, col=m) +} + +# Eref by latitude and Tmean in june +x <- seq(40, 80, 0.1) +plot(x, rep(0, length(x)), ylim=c(0,80), col="white", xlab="Latitude", ylab="Eref", main="Eref by latitude and Tmean in june") +latitude=x +m=6 +for(t in seq(-4, 16, 2)){ + tmmin=rep(t-4, length(x)) + tmmax=rep(t+4, length(x)) + y <- calc_Eref(m, tmmin, tmmax, latitude) + points(x,y, col=m, pch=m, cex=0.1) + text(40, y[1], t, col=m) + text(80, y[length(x)], t, col=m) +} + +calc_S0_I <- function(d, tm, latitude) { + # BASIC COMPUTER PROGRAM FOR ESTIMATING DAILY RA VALUES + # D=JULIAN DAY (JANUARY 1=1) + # DEC=DECLINATION OF THE SUN IN RADIANS + # ES=MEAN MONTHLY DISTANCE OF THE SUN TO THE EARTH DIVIDED BY THE MEAN ANNUAL DISTANCE + # LD=LATITUDES IN DEGREES + # LDM=MINUTES OF LATITUDES + # RA=MEAN MONTHLY EXTRATERRESTRIAL RADIATION IN MM/DAY + # RAL=MEAN MONTHLY EXTRATERRESTRIAL RADIATION IN LANGLEYS/DAY + # TC=MEAN DAILY TEMPERATURE IN DEGREE CELSIUS + Y <- cos(0.0172142 * (d + 192L)) + DEC <- 0.40876 * Y + ES <- 1.0028 + 0.03269 * Y + XLR <- latitude / 57.2958 + Z <- -tan(XLR) * tan(DEC) + OM <- -atan(Z / sqrt(-Z * Z + 1)) + pi / 2 + OM[!is.finite(OM)] <- if(m %in% 4:9) 3.1 else 0 + # CALCULATE THE DAILY EXTRATERRESTRIAL RADIATION IN LANGLEYS/DAY + DL <- OM / 0.1309 + RAL <- 120 * (DL * sin(XLR) * sin(DEC) + 7.639 * cos(XLR) * cos(DEC) * sin(OM)) / ES + # CALCULATE THE EXTRATERRESTRIAL RADIATION IN MM/DAY + RA <- RAL * 10 / (595.9 - 0.55 * tm) + + return(RA) +} + +#plot of OM by month and latitude (the only two variables it depends on ) +lat.min <- 65 +x <- seq(lat.min, 90, 0.01) +latitude=x +plot(x, rep(0, length(x)), ylim=c(0,3.2), col="white", xlab="Latitude", ylab="OM", main="OM value by latitude and month") +for(m in 1:12){ +d=day_julian[m] +Y <- cos(0.0172142 * (d + 192L)) +DEC <- 0.40876 * Y +ES <- 1.0028 + 0.03269 * Y +XLR <- latitude / 57.2958 +Z <- -tan(XLR) * tan(DEC) +OM <- -atan(Z / sqrt(-Z * Z + 1)) + pi / 2 +OM[!is.finite(OM)] <- if(m %in% 4:9) 3.14 else 0 # NB this is a modification of Hargreaves program to provide finite values above the arctic circle. +points(x,OM, col=m, pch=m, cex=0.1) +text(lat.min, OM[1], m) +# print(min(OM, na.rm=T)) +print(max(OM, na.rm=T)) +} + +library(terra) +library(climr) +library(data.table) + +#test for noram dem after loading the changes +dem <- rast("//objectstore2.nrs.bcgov/ffec/DEM/DEM_NorAm/dem_noram_lowres.tif") + +my_grid <- as.data.frame(dem, cells = TRUE, xy = TRUE) +colnames(my_grid) <- c("id", "lon", "lat", "elev") # rename column names to what climr expects + +## A simple climr query. This will return the observed 1961-1990 and 2001-2020 mean annual temperature (MAT) for the raster grid points. +climr <- downscale(xyz = my_grid, which_refmap = "refmap_climatena", vars = list_vars()) + +par(mfrow=c(1,3)) + +eref_climr <- climr$Eref +## populate the raster grid with the downscaled climate values +X <- rast(dem) # use the DEM as a template raster +X[climr[, id]] <- eref_climr # populate the raster cells with the 2001-2020 mean annual temperature (MAT) values, using the `id` field as the link. +plot(X, main="climr Eref") + +climna <- fread("C:/Users/CMAHONY/OneDrive - Government of BC/Data/ClimateNA_v750/noram_Normal_1961_1990MSY.csv") +eref_climna <- climna$Eref +eref_climna[eref_climna<0] <- 0 +## populate the raster grid with the downscaled climate values +X <- rast(dem) # use the DEM as a template raster +X[climna[, id1]] <- eref_climna # populate the raster cells with the 2001-2020 mean annual temperature (MAT) values, using the `id` field as the link. +plot(X, main="ClimateNA Eref") + +eref_climna[!is.finite(eref_climr)] <- NA + +X[climna[, id1]] <- eref_climna-eref_climr +plot(X, main="Eref Difference (ClimateNA - climr)") + +X[climr[, id]] <- climr$Eref_06 # populate the raster cells with the 2001-2020 mean annual temperature (MAT) values, using the `id` field as the link. +plot(X, main="climr Eref") + +monthcodes <- c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12") +par(mfrow=c(3,4)) +for(m in 1:12){ + X[climr[, id]] <- climr[[grep("Eref", names(climr))[m+1]]] # populate the raster cells with the 2001-2020 mean annual temperature (MAT) values, using the `id` field as the link. + plot(X, main=paste(month.name[m], "Eref")) +} + +slice <- climr[which(my_grid$lon==my_grid$lon[which.min(my_grid$lon > -120)]),] +grid_slice <- my_grid[which(my_grid$lon==my_grid$lon[which.min(my_grid$lon > -120)]),] +climr_slice <- downscale(xyz = grid_slice, which_refmap = "refmap_climatena", vars = list_vars()) + +latitude=my_grid$lat[which(my_grid$lon==my_grid$lon[which.min(my_grid$lon > -120)])] +tmmin=slice$Tmin_06 +tmmax=slice$Tmax_06 +m=6 +x <- latitude +y <- calc_Eref(m, tmmin, tmmax, latitude) +plot(latitude, slice$Eref_06, xlab="Latitude", ylab="Eref", main=paste(month.name[m], "Eref from slice of climr query on raster")) +plot(latitude, y, xlab="Latitude", ylab="Eref", main=paste(month.name[m], "Eref from climr function using climr temperatures")) +plot(latitude, climr_slice$Eref_06, xlab="Latitude", ylab="Eref", main=paste(month.name[m], "Eref from climr query of slice using climr temperatures")) +plot(x, tmmax, xlab="Latitude", ylab="Tmax", main=paste(month.name[m], "climr tmax")) + diff --git a/tests/manual_tests/check_CompareClimateNA.R b/tests/manual_tests/check_CompareClimateNA.R new file mode 100644 index 0000000..3c8ee04 --- /dev/null +++ b/tests/manual_tests/check_CompareClimateNA.R @@ -0,0 +1,67 @@ +## comparison of climr to ClimateNA for north america +## Colin Mahony + +library(terra) +library(climr) +library(data.table) + +#low res dem for north america +dem <- rast("//objectstore2.nrs.bcgov/ffec/DEM/DEM_NorAm/dem_noram_lowres.tif") + +my_grid <- as.data.frame(dem, cells = TRUE, xy = TRUE) +colnames(my_grid) <- c("id", "lon", "lat", "elev") # rename column names to what climr expects + +## climr data +climr <- downscale(xyz = my_grid, which_refmap = "refmap_climatena", vars = list_vars()) + +## climateNA data +climna_grid <- my_grid[,c(1,1,2,3,4)] +colnames(climna_grid) <- c("id1", "id2", "lon", "lat", "el") # rename column names to what climr expects +write.csv(climna_grid, "C:/Users/CMAHONY/OneDrive - Government of BC/Data/ClimateNA_v750/noram.csv") +## Run climateNA and then return to script. +climna <- fread("C:/Users/CMAHONY/OneDrive - Government of BC/Data/ClimateNA_v750/noram_Normal_1961_1990MSY.csv") + +## comparison plot +X <- rast(dem) # use the DEM as a template raster +data("variables") +figDir <- tempdir() + +vars <- list_vars() +vars <- vars[-which(vars%in%c("Tave", "Tmin", "Tmax"))] + +# for(var in vars[grep("NFFD", vars)]){ +for(var in vars){ + png(filename=paste(figDir, "/climrVclimna", var, "png",sep="."), type="cairo", units="in", width=6.5, height=2, pointsize=10, res=600) + + par(mfrow=c(1,3), mar=c(0,0,2,2)) + values(X) <- NA + + data_climr <- climr[,get(var)] + X[climr[, id]] <- data_climr + plot(X, main=paste("climr", var), axes=F) + + var_climna <- variables$Code_ClimateNA[which(variables$Code==var)] + data_climna <- climna[,get(var_climna)] + data_climna[data_climna == -9999] <- NA + X[climr[, id]] <- data_climna + plot(X, main=paste("ClimateNA", var_climna), axes=F) + + # var_type <- if(length(grep("Eref|CMD|PPT|PAS|DD|AHM|MSP|MAP", var))>0) "ratio" else "interval" + var_type <- "interval" + diff <- if(var_type=="interval") data_climna-data_climr else data_climna/data_climr + plotrange <- quantile(diff, c(0.001, 0.999), na.rm=T) + ceiling <- max(abs(plotrange)) + pal <- colorRampPalette(c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#F7F7F7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC", "#053061"))(99) + X[climna[, id1]] <- diff + plot(X, axes=F, range= if(is.finite(ceiling)) c(-ceiling, ceiling) else NULL, main= if(var_type=="interval") "Difference (ClimateNA - climr)" else "Difference (ClimateNA/climr)", col=pal) + + # library(scales) + # par(mar=c(3,3,2,2)) + # hist_climr <- hist(data_climr, xlab=var, main="") + # hist(data_climna, add=T, col=alpha("dodgerblue", 0.5), breaks=hist_climr$breaks) + # legend("topright", legend=c("climr", "ClimateNA"), fill=c("gray", alpha("dodgerblue", 0.5)), bty="n") + + print(var) + dev.off() +} + diff --git a/tests/testthat/test-climr_downscale.R b/tests/testthat/test-climr_downscale.R index 7914af7..9260219 100644 --- a/tests/testthat/test-climr_downscale.R +++ b/tests/testthat/test-climr_downscale.R @@ -76,14 +76,14 @@ test_that("test dowscale basic and spatial", { # vars = c("PPT", "CMD", "CMI", "Tave_01", "Tave_07"), # out_spatial = TRUE, plot = "CMD" # ) - # + # # ds_hist_spatial2 <- downscale( # xyz = xyz, which_refmap = "auto", # obs_periods = "2001_2020", # vars = c("PPT", "CMD", "CMI", "Tave_01", "Tave_07"), # out_spatial = TRUE, plot = "CMD" # ) - # + # # ds_hist_spatial2 <- downscale( # xyz = xyz, which_refmap = "auto", # obs_periods = "2001_2020", @@ -94,7 +94,7 @@ test_that("test dowscale basic and spatial", { # vars = c("PPT", "CMD", "CMI", "Tave_01", "Tave_07"), # out_spatial = TRUE, plot = "CMD" # ) - # + # # ds_hist_spatial2 <- downscale( # xyz = xyz, which_refmap = "auto", # obs_periods = "2001_2020", @@ -110,6 +110,7 @@ test_that("test dowscale basic and spatial", { test_that("test downscale with different argument combinations", { testInit("data.table") + testInit("terra") ## a small no. of points xyz <- data.frame( @@ -162,6 +163,7 @@ test_that("test downscale with different argument combinations", { argsCombos <- unique(argsCombos) argsCombos <- argsCombos[89:108,] + out <- apply(argsCombos, 1, function(args, xyz) { args <- args[!is.na(args)] suppressWarnings(args$xyz <- xyz) # coerces to list. @@ -186,7 +188,6 @@ test_that("test downscale with different argument combinations", { args$gcm_hist_years <- eval(parse(text = args$gcm_hist_years)) } - #browser() out <- try(do.call(downscale, args)) test <- is(out, "data.table") @@ -202,11 +203,12 @@ test_that("test downscale with different argument combinations", { test4[i] <- all(args$obs_periods %in% unique(out$PERIOD)) i <- i + 1 } - if ("obs_years" %in% checkArgs) { - test4[i] <- all(args$obs_years %in% unique(out$PERIOD)) - i <- i + 1 - } + # if ("obs_years" %in% checkArgs) { + # test4[i] <- all(args$obs_years %in% unique(out$PERIOD)) + # i <- i + 1 + # } if ("gcms" %in% checkArgs) { + message(paste(c(names(out), checkArgs, "\n"),sep = " ")) test4[i] <- all(args$gcms %in% unique(out$GCM)) i <- i + 1 } diff --git a/vignettes/climr_workflow_int.Rmd b/vignettes/climr_workflow_int.Rmd index 21d9ca1..5aa23e6 100644 --- a/vignettes/climr_workflow_int.Rmd +++ b/vignettes/climr_workflow_int.Rmd @@ -276,8 +276,8 @@ ds_out_ts <- downscale( xyz = xyzDT, which_refmap = "auto", obs_years = 2001:2015, ## currently up to 2015 - gcm_hist_years = 2001:2020, ## currently up to 2010 - gcm_ssp_years = 2021:2040, ## currently starting at 2021 + gcm_hist_years = 2001:2014, ## currently up to 2010 + gcm_ssp_years = 2015:2040, ## currently starting at 2021 gcms = "EC-Earth3", ssps = "ssp245", max_run = 1, diff --git a/vignettes/methods_ensembleSelection.Rmd b/vignettes/methods_ensembleSelection.Rmd index 52b0e71..0681fe4 100644 --- a/vignettes/methods_ensembleSelection.Rmd +++ b/vignettes/methods_ensembleSelection.Rmd @@ -47,20 +47,29 @@ SSP2-4.5, SSP3-7.0, and SSP5-8.5). eight models failed this criterion. **Criterion 4. One model per institution.** This criterion is a widely applied best practice in ensemble selection [(Leduc et al. 2016)](https://journals.ametsoc.org/view/journals/clim/29/23/jcli-d-15-0761.1.xml). -The rationale for each selection is provided below. \* CNRM-ESM2-1 (ECS -4.8^o^C) was chosen over CNRM-CM6-1 (ECS 4.9^o^C) to give preference to -the ESM over the AOGCM configuration; \* EC-Earth3 (ECS 4.3^o^C) was -arbitrarily chosen over EC-Earth3-Veg (ECS 4.3^o^C); \* INM-CM5-0 (ECS -1.9^o^C) was arbitrarily chosen over INM-CM4-8 (ECS 1.8^o^C); \* -MPI-ESM1-2-HR (ECS 3.0^o^C) was chosen over MPI-ESM1-2-LR (ECS 3.0^o^C) -due to its high resolution and availability of \>1 run for all but -SSP5-8.5. \* MIROC6 (ECS 2.6^o^C) was chosen over MIROC-ES2L (ECS -2.7^o^C) because it has more runs/scenario, higher resolution, and -because the latter has very high temperature bias over British Columbia. -\* For the purposes of this ensemble, different physics or forcing -schemes were considered different models. We used only the r\*i1p3f1 -variants of the GISS-E2-1-G model, as these had the most complete set of -scenario simulations. +The rationale for each selection is provided below. + +- CNRM-ESM2-1 (ECS 4.8^o^C) was chosen over CNRM-CM6-1 (ECS 4.9^o^C) + to give preference to the ESM over the AOGCM configuration; + +- EC-Earth3 (ECS 4.3^o^C) was arbitrarily chosen over EC-Earth3-Veg + (ECS 4.3^o^C); + +- INM-CM5-0 (ECS 1.9^o^C) was arbitrarily chosen over INM-CM4-8 (ECS + 1.8^o^C); + +- MPI-ESM1-2-HR (ECS 3.0^o^C) was chosen over MPI-ESM1-2-LR (ECS + 3.0^o^C) due to its high resolution and availability of \>1 run for + all but SSP5-8.5. + +- MIROC6 (ECS 2.6^o^C) was chosen over MIROC-ES2L (ECS 2.7^o^C) + because it has more runs/scenario, higher resolution, and because + the latter has very high temperature bias over British Columbia. + +- For the purposes of this ensemble, different physics or forcing + schemes were considered different models. We used only the r\*i1p3f1 + variants of the GISS-E2-1-G model, as these had the most complete + set of scenario simulations. **Criterion 5. No closely related models.** Models that share components were excluded, following Figure 5 of [Brunner et al. @@ -185,14 +194,14 @@ variability of individual model runs does not affect most users of climr illustrates the importance of careful ensemble selection in analysis of climate time series. -![](cmip6NA_NAM_Tmax_sm1.png) ![](cmip6NA_NAM_Tmax_sm2.png) -***Figure 1: Comparison of the ClimateNA 13-model ensemble (dashed -lines) and 8-model ensemble (solid lines)***. *The variable is mean -summer daily maximum temperature (tasmax/Tmax_sm) averaged over North -America. The 8-model ensemble excludes models with ECS outside the -IPCC-assessed 2-5^o^C range. (top) Ensemble mean projections for the -four main CMIP6 marker scenarios. (bottom) Ensemble mean and full range -(min/max for all simulations of all models) for the SSP2-4.5 scenario.* +![](cmip6NA_NAM_Tmax_sm1.png) ![](cmip6NA_NAM_Tmax_sm2.png) ***Figure 1: +Comparison of the ClimateNA 13-model ensemble (dashed lines) and 8-model +ensemble (solid lines)***. *The variable is mean summer daily maximum +temperature (tasmax/Tmax_sm) averaged over North America. The 8-model +ensemble excludes models with ECS outside the IPCC-assessed 2-5^o^C +range. (top) Ensemble mean projections for the four main CMIP6 marker +scenarios. (bottom) Ensemble mean and full range (min/max for all +simulations of all models) for the SSP2-4.5 scenario.* ### Reconciling the equilibrium climate sensitivity of the CMIP6 ensemble with observational constraints (Excerpted from Mahony et al. 2022)