Problem with CMI values

[ElizabethKleynhans]

The CMI values generated through climR seem wrong. See reproducible example below:

library(climr)
library(data.table)
library(ggplot2)

Get lat, long, and elevation for 15 locations in BC

#test_pts<-read.csv("C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\escape_clim_2020.csv")

test_pts<-data.frame(id = seq(1,15,by=1),
                      lon = c(-120.1879,-120.4258,-121.9251,-120.3030,-127.5062,-127.6785, -120.7541,-128.8055,-120.7693,-123.2985,-127.4969,-128.3576,-128.6876,-121.6685,-117.7658),
                      lat = c(59.3396, 57.4644, 59.9900, 55.2420, 54.0191, 54.1638, 51.4921, 54.8385, 54.9365, 54.0320, 54.2036, 54.5116, 54.5913, 52.4060,49.7312),
                      elev = c(441.9092,901.2709,461.7851,926.7590,1098.2932,1022.2858,1179.9183, 164.3891, 1073.8667, 740.3359,  867.5170, 309.3564,  221.9352,  715.9024, 1215.3536))

Extract climR data for the 15 locations and plot it

ds_out <- climr_downscale(xyz = test_pts, which_normal = "auto", 
                          gcm_models = c("ACCESS-ESM1-5", "MPI-ESM1-2-HR"), 
                          ssp = c("ssp245"), 
                          #gcm_period = c("2021_2040", "2041_2060","2061_2080"),
                          gcm_ts_years = 2020,
                          max_run = 3, # we want 3 individual runs for each model
                          vars = c("CMI05", "CMI06", "CMI07", "CMI08" ))

ds_out

ds_out2<-ds_out[GCM=="ACCESS-ESM1-5" & SSP=="ssp245",]
data_long_1 <- tidyr::gather(ds_out2, condition, measurement, CMI05:CMI08, factor_key=TRUE)

p <- ggplot(data_long_1, aes(x=condition, y=measurement, fill=RUN)) + 
  labs(title="climR data",x="CMI", y = "")+
  geom_boxplot()
p

Extract climateBC data for the 15 locations and plot it

test_pts<-test_pts[,c("id","lat", "lon", "elev")]

write.csv(test_pts,'C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\test_pts.csv')

setwd("D:/Climatebc_v742"); # set the ClimateBC root directory as the working directory
exe <- "ClimateBC_v7.42.exe"

## 2020
inputFile = '/C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\test_pts.csv' 
outputFile = '/C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\output\\test_pts_2020.csv'
yearPeriod = '/Year_2020.ann'
system2(exe,args= c('/M', yearPeriod, inputFile, outputFile))

test_pts3<-read.csv('C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\output\\test_pts_2020.csv')

test_pts3<-as.data.table(test_pts3)
test_pts3<-test_pts3[,c("id", "Latitude", "Longitude", "Elevation","CMI05", "CMI06", "CMI07", "CMI08")]
summary(test_pts3)

data_long_2 <- tidyr::gather(test_pts3, condition, measurement, CMI05:CMI08, factor_key=TRUE)

# Basic box plot
p <- ggplot(data_long_2, aes(x=condition, y=measurement)) + 
  labs(title="ClimateBC data",x="CMI", y = "")+
  geom_boxplot()
p

Plot the climateBC data and the climR data together in one figure

data_long_1a<-data_long_1[,c("id", "RUN","condition","measurement")]
data_long_2a<-data_long_2[,c("id","condition","measurement")]
data_long_2a$RUN<-"ClimateBC"

data_long<-rbind(data_long_1a,data_long_2a)

p <- ggplot(data_long, aes(x=condition, y=measurement, fill=RUN)) + 
  labs(title="climR data",x="CMI", y = "")+
  geom_boxplot()
p

[cmahony]

@kdaust I modified the reprex to include climr observed values. I'm not sure there actually is a problem. I'll look into what values we would expect for BC.

library(climr)
library(data.table)
library(ggplot2)

# Get lat, long, and elevation for 15 locations in BC

#test_pts<-read.csv("C:\\Work\\caribou\\castor\\R\\fire_sim\\tmp\\escape_clim_2020.csv")

test_pts<-data.frame(id = seq(1,15,by=1),
                     lon = c(-120.1879,-120.4258,-121.9251,-120.3030,-127.5062,-127.6785, -120.7541,-128.8055,-120.7693,-123.2985,-127.4969,-128.3576,-128.6876,-121.6685,-117.7658),
                     lat = c(59.3396, 57.4644, 59.9900, 55.2420, 54.0191, 54.1638, 51.4921, 54.8385, 54.9365, 54.0320, 54.2036, 54.5116, 54.5913, 52.4060,49.7312),
                     elev = c(441.9092,901.2709,461.7851,926.7590,1098.2932,1022.2858,1179.9183, 164.3891, 1073.8667, 740.3359,  867.5170, 309.3564,  221.9352,  715.9024, 1215.3536))

# Extract climR data for the 15 locations and plot it

ds_out <- climr_downscale(xyz = test_pts, which_normal = "auto", 
                          gcm_models = c("ACCESS-ESM1-5"), 
                          ssp = c("ssp245"), 
                          #gcm_period = c("2021_2040", "2041_2060","2061_2080"),
                          gcm_ts_years = 2015,
                          historic_ts = 2015,
                          max_run = 3, # we want 3 individual runs for each model
                          vars = c("CMI05", "CMI06", "CMI07", "CMI08" ))

ds_out

ds_out2<-ds_out[PERIOD=="2015",]
data_long_1 <- tidyr::gather(ds_out2, condition, measurement, CMI05:CMI08, factor_key=TRUE)

p <- ggplot(data_long_1, aes(x=condition, y=measurement, fill=RUN)) + 
  labs(title="climR data",x="CMI", y = "")+
  geom_boxplot()
p

[cmahony]

the points are distributed through BC:

[cmahony]

Wang, Hogg, et al. 2014 suggest that the climr values are reasonable (see "Montane" and "Maritime" in figure below.)

[cmahony]

@ElizabethKleynhans have a look at my comments above. It may be that the problem is actually in ClimateBC. Did you have a priori expectation of different values than what you are getting from climr?

[ElizabethKleynhans]

Hmmm, ok thanks for looking into this. We did an analysis using the climateBC data (assuming it was correct) and then I used climr data to project forward but my projections did not make sense ... which is why I looked into this. How far back do climr's observed values go? We need data from 2009

[cmahony]

@ElizabethKleynhans the climr observed time series is only 1902-2015 but we will be updating to 1901-2023 in the next couple of weeks.

[kdaust]

Thanks both for looking into this. @cmahony I'm not sure how the CMI values in ClimateBC are calculated - I discussed with Will the best approach for calculating CMI, and we decided on the Hogg's 1997 approach (https://doi.org/10.1016/S0168-1923(96)02380-5).

[cmahony]

Things look fine but i'll keep this active until we definitively resolve it.

[ElizabethKleynhans]

Ok sounds good. Ill flag it with Tongli, since I guess that might be where the issue is. ... or at least maybe they calculate CMI in a different way. Either way it would be good to know why the values are so different.

[cmahony]

thanks for taking this on @ElizabethKleynhans! let us know what you hear from tongli.

[cmahony]

@ElizabethKleynhans @kdaust looks like a problem in climr. the values of annual CMI are too high:

here are the values we would expect, based on NRCan data published at https://cfs.nrcan.gc.ca/fc-data-catalogue/read/1

so climr is >10x too high, and missing the negative values in the BG zone that we would expect.

I've looked into the climr equations and it is probably in the calc_PET or calc_SVP functions I'm out of time to take this further today but i'll look into it more later this week.

[ElizabethKleynhans]

Thanks for looking into this again Colin! It will be great to know where the issue is coming from.

[cmahony]

there were two problems:

1. a typo in Hogg (1997) published Equation 3.
2. Hogg's CMI is in cm. we weren't dividing by ten.

The typo is in the calculation of VPD: VPD is saturated vapour pressure minus actual vapor pressure (es-ea). Hogg means to take the average of eTmin and eTmax, but subtracted instead of added.

this is the CMI without correcting the typo (black symbols are climateBC, red are calculations from equations):

here it is with correction of the typo (raw calculations perfectly match ClimateBC values):

[cmahony]

fixed in pull request: https://github.com/bcgov/climr/pull/261. for review by kiri before merge