Hi @cmahony,

Can you add @meztez and @FrankBornais as contributor to this repo please.

Thanks

@meztez the variables are out of order in the output of downscale(). e.g., [81] "Eref07" "CMD07" "Eref08" "CMD08" "CMD_sm" "Eref09" "CMD09" "Eref10"

vars.all <- list_variables()
climr <- data.frame(downscale(
  xyz = xyz,
  normal = normal,
  gcm = NULL,
  var = vars.all
))
names(climr)

> names(climr)
  [1] "ID"       "PPT01"    "PPT02"    "PPT03"    "PPT04"    "PPT05"    "PPT06"    "PPT07"   
  [9] "PPT08"    "PPT09"    "PPT10"    "PPT11"    "PPT12"    "Tmax01"   "Tmax02"   "Tmax03"  
 [17] "Tmax04"   "Tmax05"   "Tmax06"   "Tmax07"   "Tmax08"   "Tmax09"   "Tmax10"   "Tmax11"  
 [25] "Tmax12"   "Tmin01"   "Tmin02"   "Tmin03"   "Tmin04"   "Tmin05"   "Tmin06"   "Tmin07"  
 [33] "Tmin08"   "Tmin09"   "Tmin10"   "Tmin11"   "Tmin12"   "MWMT"     "MCMT"     "MAT"     
 [41] "MAP"      "TD"       "NFFD12"   "NFFD01"   "NFFD02"   "NFFD_wt"  "NFFD03"   "NFFD04"  
 [49] "NFFD05"   "NFFD_sp"  "NFFD06"   "NFFD07"   "NFFD08"   "NFFD_sm"  "NFFD09"   "NFFD10"  
 [57] "NFFD11"   "NFFD_at"  "NFFD"     "bFFP"     "Eref12"   "CMD12"    "Eref01"   "CMD01"   
 [65] "Eref02"   "CMD02"    "CMD_wt"   "Eref03"   "CMD03"    "Eref04"   "CMD04"    "Eref05"  
 [73] "CMD05"    "CMD_sp"   "Eref06"   "CMD06"    "Eref07"   "CMD07"    "Eref08"   "CMD08"   
 [81] "CMD_sm"   "Eref09"   "CMD09"    "Eref10"   "CMD10"    "Eref11"   "CMD11"    "CMD_at"  
 [89] "CMD"      "Tmin_wt"  "Tmin_sp"  "Tmin_sm"  "Tmin_at"  "Tmin"     "Tmax_wt"  "Tmax_sp" 
 [97] "Tmax_sm"  "Tmax_at"  "Tmax"     "CMI"      "CMI_at"   "CMI_sm"   "CMI_sp"   "CMI_wt"  
[105] "CMI01"    "CMI02"    "CMI03"    "CMI04"    "CMI05"    "CMI06"    "CMI07"    "CMI08"   
[113] "CMI09"    "CMI10"    "CMI11"    "CMI12"    "Tave12"   "DD_0_12"  "Tave01"   "DD_0_01" 
[121] "Tave02"   "DD_0_02"  "DD_0_wt"  "Tave03"   "DD_0_03"  "Tave04"   "DD_0_04"  "Tave05"  
[129] "DD_0_05"  "DD_0_sp"  "Tave06"   "DD_0_06"  "Tave07"   "DD_0_07"  "Tave08"   "DD_0_08" 
[137] "DD_0_sm"  "Tave09"   "DD_0_09"  "Tave10"   "DD_0_10"  "Tave11"   "DD_0_11"  "DD_0_at" 
[145] "DD_0"     "DD_18_12" "DD_18_01" "DD_18_02" "DD_18_wt" "DD_18_03" "DD_18_04" "DD_18_05"
[153] "DD_18_sp" "DD_18_06" "DD_18_07" "DD_18_08" "DD_18_sm" "DD_18_09" "DD_18_10" "DD_18_11"
[161] "DD_18_at" "DD_18"    "DD18_12"  "DD18_01"  "DD18_02"  "DD18_wt"  "DD18_03"  "DD18_04" 
[169] "DD18_05"  "DD18_sp"  "DD18_06"  "DD18_07"  "DD18_08"  "DD18_sm"  "DD18_09"  "DD18_10" 
[177] "DD18_11"  "DD18_at"  "DD18"     "DD5_12"   "DD5_01"   "DD5_02"   "DD5_wt"   "DD5_03"  
[185] "DD5_04"   "DD5_05"   "DD5_sp"   "DD5_06"   "DD5_07"   "DD5_08"   "DD5_sm"   "DD5_09"  
[193] "DD5_10"   "DD5_11"   "DD5_at"   "DD5"      "eFFP"     "EMT"      "Eref_wt"  "Eref_sp" 
[201] "Eref_sm"  "Eref_at"  "Eref"     "EXT"      "FFP"      "PAS12"    "PAS01"    "PAS02"   
[209] "PAS_wt"   "PAS03"    "PAS04"    "PAS05"    "PAS_sp"   "PAS06"    "PAS07"    "PAS08"   
[217] "PAS_sm"   "PAS09"    "PAS10"    "PAS11"    "PAS_at"   "PAS"      "PPT_wt"   "PPT_sp"  
[225] "PPT_sm"   "PPT_at"   "PPT"      "SHM"      "Tave_wt"  "Tave_sp"  "Tave_sm"  "Tave_at" 
[233] "Tave"    
> 

@meztez the shiny app wasn't one of the required deliverables but i'm curious about its status and potential for further development.

• Document individual R functions
• Document individual C++
• Produce beautiful package vignettes that explain how to use package

• Understand climate data format provided
• Create maintenance functions to structure climate data
• Create climate data pipeline from github repo to local storage
• Cache lapse rates in memory for R session

References:
https://github.com/bcgov/bcmaps/blob/master/R/cache-utils.R

dem, historical, future will be retrieved from repo

functions for each of the following methods : list available, retrieve, build compatible object to feed back into downscale / lapse rate functions.

I received this message from Sarah Endicott @see24. this is out of scope for Boostao, its an issue for future work.

I noticed that you have the derived variables from Wang 2016 in your repo and I wanted to give you a heads up that I found an issue in that paper when I implemented those functions in R a couple years ago.
I emailed Tongli about it and he agreed the paper was wrong.

First of all in the paper the DD5 equation looks like this:

I believe the directions of the > and < are reversed because it should be a straight line at temperatures above k and a curve heading to zero below k.

Secondly, when I calculated DD5 using the parameters in the paper the DD5 was higher than it should be for many locations. When I investigated the issue I found that a graph made using the parameters in the paper did not match the graph for DD5 vs Jan mean temperature in the paper (see graph below left, blue line). I discovered that if I exchanged the parameters b and T0 from the S2 Table the graph appeared to match the paper (see graphs below red lines, and table below).

Also if you haven't coded the equations yet I might be able to share my code but would have to check with my employer so let me know if you are interested!

@meztez This is a running list of feedback on the test version from Will, Kiri, and Colin. made with the understanding that the package is still in dev and some of these items are already on your mind.

• some of the derived variables don't show up when requested in downscale() e.g., MAT, MCMT. seems to be the annual ones.
• temperatures seem to be returned as degrees*10. any reason for this?
• all months are returned for the basic variables, plus additional requested variables. is there a way to specify only the desired variables? e.g., just Tmax01 and DD5_01. if not, we need this.

Some of the derived variable equations aren't quite right. e.g., MWMT should be mean warmest month temperature (i.e., Tave of the month with the highest Tave for each location). currently it is calculated as annual mean Tmax. same issue for MCMT. this affects some of the downstream variables.

"MWMT" = function() {set(dt, j = "MWMT", value = (v("Tmax01")+v("Tmax02")+v("Tmax03")+v("Tmax04")+v("Tmax05")+v("Tmax06")+v("Tmax07")+v("Tmax08")+v("Tmax09")+v("Tmax10")+v("Tmax11")+v("Tmax12"))/12)},

@meztez I'm not sure if this is straightforward to calculate since in theory the warmest month could differ among cells, and would do so at continental scales. this is perhaps something i could fix as future work, but if you have bandwidth for it that would be great.

something wrong with precipitation as snow. it should be zero in summer but lots of snow happening in climR. e.g. in july below:

winter PAS looks fine.

It does not seem to be supported.

rspatial/terra#298.

I would keep the csv index for now.

climR is structured nicely to scale up to North America. this would be done by:

1. subbing in GCM files for the north american domain
2. adding a 2-4km grid of climate normals for north america, and 800m (native PRISM) normals for BC.

This approach would allow the user to select the grid they want, or supply an 800m PRISM grid of their own in the R environment, e.g. for study areas in the US.

@meztez let us know if you see any problems with this idea or additional considerations/steps.

there are typos returned in list_gcm() that create some duplicates (e.g., "MPI-ESM1-2-HR" "MPI-ESMI-2-HR").

> list_gcm()
 [1] "ACCESS-ESM1-5" "BCC-CSM2-MR"   "CanESM5"       "CNRM-ESM2-1"   "EC-Earth3"     "GFDL-ESM4"    
 [7] "GISS-E2-1-G"   "INM-CM5-0"     "IPSL-CM6A-LR"  "MIROC6"        "MPI-ESM1-2-HR" "MPI-ESMI-2-HR"
[13] "MRI-ESM2-0"    "UKESM1-0-LL"   "UKESMI-0-LL"  
> unique(results$GCM)
 [1] "ACCESS-ESM1-5" "BCC-CSM2-MR"   "CanESM5"       "CNRM-ESM2-1"   "EC-Earth3"     "GFDL-ESM4"    
 [7] "GISS-E2-1-G"   "INM-CM5-0"     "IPSL-CM6A-LR"  "MIROC6"        "MPI-ESM1-2-HR" "MRI-ESM2-0"   
[13] "UKESM1-0-LL"  
> 

Input

• Points (x, y, z)
• GCM (vector)
• Scenarios (vector)
• Maximum number of runs to include (pick in order)
• Base reference layer (historical and lapse rates)
• Climate variables to compute

Output

• Wide on the climate variables
• Row per gcm / scenario / run / period ...

In the parameters table for DD<0 sometime the value of k is missing. The computation of DD_m requires k to determine which formula to use:

Which formula should we use when k is missing?

ref : DD<0 table

Using January Tmax as an example, we can see that the gradients in the mountains (eastern portion) are a bit different in ClimateNA and ClimateBC. Specifically, ClimateBC is indicating some inversions (colder valley bottom than mountain side) in the Quesnel Mountains and some other areas, which is realistic. seems that the inversions in the 800m PRISM are unresolved in the 4km ClimateNA grid. we may want to reduce our domain and use the 800m PRISM grid, or in the longer term assemble a composite grid of western North America.

TL;DR

Topics greatly improve the discoverability of repos; please add the short code from the table below to the topics of your repo so that ministries can use GitHub's search to find out what repos belong to them and other visitors can find useful content (and reuse it!).

Why Topic

In short order we'll add our 800th repo. This large number clearly demonstrates the success of using GitHub and our Open Source initiative. This huge success means it's critical that we work to make our content as discoverable as possible. Through discoverability, we promote code reuse across a large decentralized organization like the Government of British Columbia as well as allow ministries to find the repos they own.

What to do

Below is a table of abbreviation a.k.a short codes for each ministry; they're the ones used in all @gov.bc.ca email addresses. Please add the short codes of the ministry or organization that "owns" this repo as a topic.

That's it, you're done!!!

How to use

Once topics are added, you can use them in GitHub's search. For example, enter something like org:bcgov topic:citz to find all the repos that belong to Citizens' Services. You can refine this search by adding key words specific to a subject you're interested in. To learn more about searching through repos check out GitHub's doc on searching.

Pro Tip 🤓

• If your org is not in the list below, or the table contains errors, please create an issue here.

• While you're doing this, add additional topics that would help someone searching for "something". These can be the language used javascript or R; something like opendata or data for data only repos; or any other key words that are useful.

• Add a meaningful description to your repo. This is hugely valuable to people looking through our repositories.

• If your application is live, add the production URL.

Ministry Short Codes

NOTE See an error or omission? Please create an issue here to get it remedied.

@meztez Precipitation as snow doesn't seem to be working. it is zero for all months.

Precipitation seems to be working as expected with a lapse rate. However, in this case we shouldn't be using a lapse rate. in reality, precipitation doesn't vary with elevation at scales less than a couple of kilometers. ClimateBC doesn't apply a lapse rate to precipitation for this reason. the fix for this is to just to use bilinear interpolation from the baseline precipitation rasters, rather than applying a lapse rate to precipitation. @meztez, it would be great to make this modification. ideally, we would keep the option of applying a lapse rate to precipitation, perhaps with a command in the downscale() function, defaulted off.

PAS is higher in climR than in climateBC. i think this is probably the parameters; the parameters in the ClimateNA publication had some errors and tongli provided some WNA parameters that may not match these. Bruno corrected some errors in the published equations, but the parameters may still be an issue.

TL;DR 🏎️

Teams are encouraged to favour modern inclusive phrasing both in their communication as well as in any source checked into their repositories. You'll find a table at the end of this text with preferred phrasing to socialize with your team.

Words Matter

We're aligning our development community to favour inclusive phrasing for common technical expressions. There is a table below that outlines the phrases that are being retired along with the preferred alternatives.

During your team scrum, technical meetings, documentation, the code you write, etc. use the inclusive phrasing from the table below. That's it - it really is that easy.

For the curious mind, the Public Service Agency (PSA) has published a guide describing how Words Matter in our daily communication. Its an insightful read and a good reminder to be curious and open minded.

What about the master branch?

The word "master" is not inherently bad or non-inclusive. For example people get a masters degree; become a master of their craft; or master a skill. It's generally when the word "master" is used along side the word "slave" that it becomes non-inclusive.

Some teams choose to use the word main for the default branch of a repo as opposed to the more commonly used master branch. While it's not required or recommended, your team is empowered to do what works for them. If you do rename the master branch consider using main so that we have consistency among the repos within our organization.

Preferred Phrasing

Pro Tip 🤓

This list is not comprehensive. If you're aware of other outdated nomenclature please create an issue (PR preferred) with your suggestion.

@meztez There are a few variables still missing from the variable set:

• Rad/MAR - ignore, we're not doing radiation
• RH - i recall you asking about this. any issue with calculating it?
• MSP - may through september precipitation. just need to sum across those months.
• AHM - annual heat-moisture index - (MAT+10)/(MAP/1000)
• DD1040 - degree-days above 10°C and below 40°C. ignore, as we don't yet have the equation for this one.

# climR data
climr <- data.frame(downscale(
  xyz = xyz,
  normal = normal,
  gcm = NULL,
  var = vars.all
))

# ClimateBC data
library(data.table)
fplot=paste("C:\\Users\\CMAHONY\\OneDrive - Government of BC\\Projects\\2022_FLPapp\\inputs\\", studyarea, "_", scenario, "_", proj.year, "_MSY.csv", sep="")
Y0 <- fread(fplot, stringsAsFactors = FALSE, data.table = FALSE) #fread is faster than read.csv
Y0 <- separate(Y0, Year, into = c("GCM","scenario","proj.year"), sep = "_", remove = T)
Y0$proj.year <- gsub(".gcm","",Y0$proj.year)
Y0$proj.year <- substr(Y0$proj.year, 1,4)

climbc <- Y0[which(Y0$GCM==gcm),]

vars.missing <- names(climbc)[-c(1:8)][-which(names(climbc)[-c(1:8)]%in%names(climr)[-1])]

> vars.missing
 [1] "Rad01"  "Rad02"  "Rad03"  "Rad04"  "Rad05"  "Rad06"  "Rad07"  "Rad08"  "Rad09"  "Rad10" 
[11] "Rad11"  "Rad12"  "RH01"   "RH02"   "RH03"   "RH04"   "RH05"   "RH06"   "RH07"   "RH08"  
[21] "RH09"   "RH10"   "RH11"   "RH12"   "Rad_wt" "Rad_sp" "Rad_sm" "Rad_at" "RH_wt"  "RH_sp" 
[31] "RH_sm"  "RH_at"  "MSP"    "AHM"    "MAR"    "RH"     "DD1040"
>

Inputs can be compressed further to increase data update speed.

Normals can be stacked in a single compressed NetCDF file instead of multiple ascii for better handling.

• Create reformat script
• Rename inputs dir raw + pkgm add to propre .RBuildignore
• Refactor pkg to handle optimized format

I'm wondering if data size and speed can be improved by limiting decimal precision. don't need more than 2 decimals on temperature and 1 decimal on the other variables. just a thought. no action needed at this point

ClimateBC/NA v7.3 includes an R interface. it would be good to do a speed test of climR against this interface.

@cmahony The provided digital elevation model has the same resolution as the target climate variables (tasmin, tasmax and pr). In order to apply the lapse rate adjustment, a better resolution DEM would be best. If I use the provided one, the adjustment would be zero all the time. I could use the same one as the CCISS tool, but if you have access to a better one, even beter. What do you think?

CCISS DEM :

Digital Elevation Model (DEM) for British Columbia produced by GeoBC. This data is the TRIM DEM converted to the Canadian Digital Elevation Data (CDED) format. The data consists of an ordered array of ground or reflective surface elevations, recorded in metres, at regularly spaced intervals. The spacing of the grid points is .75 arc seconds north/south. The data was converted into 1:50,000 grids for distribution. The scale of this modified data is 1:250,000 which was captured from the original source data which was at a scale of 1:20,000.

when i do a comparison of ClimateBC future values (here ssp245 and 2041-2060) the terrain effects match perfectly, but there are coarse-scale differences.

These biases are different for each gcm

need to resolve this

https://github.com/bcgov/climR-pnw/blob/138d82cc53043bfdf3cc86cbf8e620bae0db75ed/scripts/lapserate.R#L126

Some lapse rate computations will produce NaN as the surrounding cell deltas are all 0s. Instead of NaN, I'm replacing them with 0s, so they won't impact vectorized downscaling. Otherwise, I would have to check for NaN within the loop.

@cmahony Are if you are ok with this hypothesis.

the values of the climate variables aren't smooth across a smooth elevation gradient. I'm guessing there is a small problem with the lapse rates.

e.g. Tave over a test dem:

this is the dem itself:

here is code to reproduce the example (the dem is in the base code folder):

dem <- raster(paste("dem_test.tif", sep=""))
X <- dem
values(X) <- NA

dem.pts <- rasterToPoints(dem, spatial=T) 
dem.pts <- as.data.frame(dem.pts) 
land <- which(!is.na(values(dem)))
xyz <- data.frame(lon=dem.pts$x, lat=dem.pts$y, elev=dem.pts[,1])

# Select GCM
gcm <- gcm_input(
  gcm = gcms[5],
  ssp = c("ssp245"),
  period = "2041_2060",
  max_run = 0
)

# projected normals
projnormals <- data.frame(downscale(
  xyz = xyz,
  normal = normal,
  gcm = gcm,
  var = vars.all
))

var="Tave"
par(mar=c(0,0,0,1))
for(var in names(projnormals)[-c(1:5)]){
  run.0 <- projnormals[which(projnormals$GCM=="EC.Earth3" 
                          & projnormals$SSP=="ssp245" 
                          & projnormals$RUN=="ensembleMean" 
                          & projnormals$PERIOD=="2041_2060"), 
                    which(names(projnormals)==var)]
  
  values(X)[land] <- run.0
  plot(X)
  mtext(var, cex=1.5, line=-1.5, adj=0.05)
  
}

the MAP for climR and ClimateBC is different. my guess is that precipitation is being averaged over the year in climR, but it should be summed over the year. [we are going with climateBC naming conventions for now, but MAP is a terrible name.]

]

• Reimplement logic in C++ using Eigen
• Build performance benchmarks tools
• In case performance are not satisfying, evaluate alternative C++ libraries (unnecessary)
• Make sure performance are similar on Windows

Rename baseline, future, some arguments for consistency.
Pass all parameters to function within function instead of relying on scoping.

Main computation R interface

• Define user-functions interface
• Possibility utilize class dispatch to adapt to multiple input types (sf points, raster, data.frame, matrix)
• Document logic steps
• Quickly code a non optimized algorithm to do steps in R / Julia, anything as long as it is quick and it works.
• Validate logic
• Create tests to make sure logic is respected

Move lapse_rate to normal / historical object creation as they are really tied together.

Refactor to account for it.

Cover

• Lapse rate
• Package usage
• Extending support for other data sources (moving from github)
• Any other subject that might require more details

@meztez i just noticed that i had put the updated gcm files in the \inputs folder, rather than replacing the old ones in the \inputs\gcmData folder. not sure where you are at with this, so i'll leave it to you to move the new files into the \inputs\gcmData and delete the old ones. sorry about that.

This issue is a kind reminder that your repository has been inactive for 181 days. Some repositories are maintained in accordance with business requirements that infrequently change thus appearing inactive, and some repositories are inactive because they are unmaintained.

To help differentiate products that are unmaintained from products that do not require frequent maintenance, repomountie will open an issue whenever a repository has not been updated in 180 days.

• If this product is being actively maintained, please close this issue.
• If this repository isn't being actively maintained anymore, please archive this repository. Also, for bonus points, please add a dormant or retired life cycle badge.

Thank you for your help ensuring effective governance of our open-source ecosystem!

@meztez Something is wrong with Tmin. it doesn't have the topographical detail that we would expect from applying a lapse rate to the DEM grid. note that this blurriness is not there for Tmax and PPT. I've shown january here, but the problem is there for all months.

@meztez We need a function to downscale the 1961-1990 normals. i see that this isn't specifically stated in the CWU opportunity but it is a core function of ClimateNA and is needed to do analysis on the data. perhaps this is already done but i don't know the usage.

It seems that Tmin doesn’t have a lapse rate applied, while PPT does have a lapse rate applied. You can see this in the January comparisons of ClimateBC and ClimateNA for each variable, attached, but this pattern is the same in all months. I’m guessing that the ClimateNA code to skip elevation adjustment in PPT is being applied to Tmin instead. I have raised this issue with Tongli and waiting for a response. looks like we'll have to wait for a corrected version of ClimateNA and sub in the new baseline.

Tmin:

PPT

Finish append calculations methods.

It receives a data.table with PPT01-PPT12, Tmax01-Tmax12, Tmin01-Tmin12, a vector of extra climate variables code, and grouping?

Grouping could probably be skipped and variables would match Codes in provided file.

The observed 2001-2020 climate is used in many of our analyses. to get this into climR we would provide gridded anomalies in whatever grid resolution, and then modify downscale() so that it can accept these non-GCM anomalies.
@meztez this wasn't specified in the CWU opportunity so is out of scope for this phase.

@meztez the hyphens in the gcm name get turned into periods in downscale(). this will be inconvenient for querying the data.

> gcm <- gcm_input(
+   gcm = c("BCC-CSM2-MR"),
+   ssp = c("ssp245"),
+   period = "2041_2060",
+   max_run = 0
+ )
> results <- downscale(
+   xyz = xyz,
+   normal = normal,
+   gcm = gcm,
+   var = vars.all
+ )
> unique(results$GCM)
[1] "BCC.CSM2.MR"

@meztez if you want to visualize the climate model data you are working with, see an app i made cmip-NA .

when using the BC normals, the lapse rates at the southern, eastern, and northern borders of BC might have some artefacts.

The current baseline is a 2km grid based on a ClimateNA query. so we are essentially elevation-adjusting from a 4km grid. we are mainly taking this approach because we don't have an 800m composite of the PRISM + WorldClim for the western north america study area. However, there are some climatic details (e.g. inversions) that get washed out in the 4km grid. it would be good in the short term to provide the option to use the 800m BC PRISM, for analyses that are within BC. this would involve back-casting the 1971-2000 PRISM normals to 1961-1990.

@meztez we need a list_runs() function that will show the runs (e.g. r1i1p1f1, r2i2p2f2...) available for each GCM.

Huge discrepancies with ClimateBC in NFFD. perhaps this is related to the problem with Tmin. or it might be an error in the published equations.

• Create package skeleton
• Configure .gitignore, .Rbuildignore, special folders + vignette

No Project Lifecycle Badge found in your readme!

Hello! I scanned your readme and could not find a project lifecycle badge. A project lifecycle badge will provide contributors to your project as well as other stakeholders (platform services, executive) insight into the lifecycle of your repository.

What is a Project Lifecycle Badge?

It is a simple image that neatly describes your project's stage in its lifecycle. More information can be found in the project lifecycle badges documentation.

What do I need to do?

I suggest you make a PR into your README.md and add a project lifecycle badge near the top where it is easy for your users to pick it up :). Once it is merged feel free to close this issue. I will not open up a new one :)