Modelling gridMET climate for CAN PIAL WPBR transects from Daymet data
Stephen Huysman
2023-09-14
Introduction
Jamison et al. (In Preparation) developed a range-wide model of probability of White Pine Blister Rust infection in whitebark pine based on August-September Relative Humidity and Temperature. This model was developed using Daymet historical climate data Thornton et al. (2014). We wish to develop an analog of this model based on historical gridMET historical climate data (Abatzoglou 2013), so that it can be used with cliamte projections based on gridMET and MACA (Abatzoglou and Brown 2012) such as those developed by Tercek et al. (2021).
The dataset analyzed by Jamison et al. (In Preparation) contains transects located in parks in Canada, which are covered by Daymet historical climate data but not by gridMET. Here, we attempt to see if Daymet climate data can reasonably be modeled from gridMET data. The approach is to fit a linear regression to Daymet temperature data with gridMET temperature data as the response variable. We select Glacier NP transects to use to develop the model as they are the nearest transects to the Canadian Park transects and thus assumed to be most climatically similar.
Setup
library(daymetr)
library(here)
library(lubridate)
library(tidyverse)
library(modelr)
## Transect location data provided by EJ
parks <- read_csv("SITE_LOCATIONS.csv")## Rows: 745 Columns: 5
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): network, park, transect_id
## dbl (2): lat, long
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.## All Canadian transects
trans_can <- parks %>% filter(park == "CAN")
## All Glacier NP transects
trans_glac <- parks %>% filter(park == "GLAC")
##trans_glac <- parks %>% filter(park != "CAN")
gridmet_dir <- file.path("out")
## Load gridMET historical data for GLAC transects. These were downloaded previously using gridmet-historical.R script
gridmet_glac_files <- file.path(gridmet_dir, paste0(trans_glac$transect_id, ".csv"))
gridmet_glac_data <- do.call(rbind, lapply(gridmet_glac_files, function(x)
transform(read.csv(x), transect_id = basename(x))))
head(gridmet_glac_data); tail(gridmet_glac_data)## X Date P TmaxK TminK RHmax RHmin VPD RH Tmax Tmin T transect_id
## 1 1 1979-01-01 0.0 252.8 242.6 91.2 56.8 0.02 74.00 -20.35 -30.55 -25.45 WB0401.csv
## 2 2 1979-01-02 0.5 258.2 247.4 90.7 57.0 0.03 73.85 -14.95 -25.75 -20.35 WB0401.csv
## 3 3 1979-01-03 0.0 259.2 250.0 99.4 53.5 0.05 76.45 -13.95 -23.15 -18.55 WB0401.csv
## 4 4 1979-01-04 0.0 254.2 241.8 91.7 46.2 0.04 68.95 -18.95 -31.35 -25.15 WB0401.csv
## 5 5 1979-01-05 0.0 256.1 242.0 83.6 42.3 0.05 62.95 -17.05 -31.15 -24.10 WB0401.csv
## 6 6 1979-01-06 0.0 259.1 244.8 86.9 51.4 0.06 69.15 -14.05 -28.35 -21.20 WB0401.csv## X Date P TmaxK TminK RHmax RHmin VPD RH Tmax Tmin T transect_id
## 1074805 16280 2023-07-28 0 289.6 278.5 64.1 33.5 0.75 48.80 16.45 5.35 10.90 WC-002.csv
## 1074806 16281 2023-07-29 0 294.3 279.6 72.1 30.1 0.97 51.10 21.15 6.45 13.80 WC-002.csv
## 1074807 16282 2023-07-30 0 297.7 282.8 60.1 25.1 1.35 42.60 24.55 9.65 17.10 WC-002.csv
## 1074808 16283 2023-07-31 0 297.7 283.9 40.7 19.6 1.58 30.15 24.55 10.75 17.65 WC-002.csv
## 1074809 16284 2023-08-01 0 296.2 284.7 30.2 17.1 1.67 23.65 23.05 11.55 17.30 WC-002.csv
## 1074810 16285 2023-08-02 0 299.6 284.5 41.3 14.9 1.82 28.10 26.45 11.35 18.90 WC-002.csvgridmet_glac_data$transect_id = substr(gridmet_glac_data$transect_id,1,nchar(gridmet_glac_data$transect_id)-4)
head(gridmet_glac_data); tail(gridmet_glac_data)## X Date P TmaxK TminK RHmax RHmin VPD RH Tmax Tmin T transect_id
## 1 1 1979-01-01 0.0 252.8 242.6 91.2 56.8 0.02 74.00 -20.35 -30.55 -25.45 WB0401
## 2 2 1979-01-02 0.5 258.2 247.4 90.7 57.0 0.03 73.85 -14.95 -25.75 -20.35 WB0401
## 3 3 1979-01-03 0.0 259.2 250.0 99.4 53.5 0.05 76.45 -13.95 -23.15 -18.55 WB0401
## 4 4 1979-01-04 0.0 254.2 241.8 91.7 46.2 0.04 68.95 -18.95 -31.35 -25.15 WB0401
## 5 5 1979-01-05 0.0 256.1 242.0 83.6 42.3 0.05 62.95 -17.05 -31.15 -24.10 WB0401
## 6 6 1979-01-06 0.0 259.1 244.8 86.9 51.4 0.06 69.15 -14.05 -28.35 -21.20 WB0401## X Date P TmaxK TminK RHmax RHmin VPD RH Tmax Tmin T transect_id
## 1074805 16280 2023-07-28 0 289.6 278.5 64.1 33.5 0.75 48.80 16.45 5.35 10.90 WC-002
## 1074806 16281 2023-07-29 0 294.3 279.6 72.1 30.1 0.97 51.10 21.15 6.45 13.80 WC-002
## 1074807 16282 2023-07-30 0 297.7 282.8 60.1 25.1 1.35 42.60 24.55 9.65 17.10 WC-002
## 1074808 16283 2023-07-31 0 297.7 283.9 40.7 19.6 1.58 30.15 24.55 10.75 17.65 WC-002
## 1074809 16284 2023-08-01 0 296.2 284.7 30.2 17.1 1.67 23.65 23.05 11.55 17.30 WC-002
## 1074810 16285 2023-08-02 0 299.6 284.5 41.3 14.9 1.82 28.10 26.45 11.35 18.90 WC-002gridmet_glac_data$Date = ymd(gridmet_glac_data$Date)
## Daymet data dir. Used for subsequent routine to download Daymet data or to load in previously downloaded data
daymet_dir <- file.path("daymet")Download Daymet data
## Generate properly formatted csv of locations to run download_daymet_batch
## Only needs to be ran once
locations <- rbind(
data.frame(site = trans_can$transect_id, lat = trans_can$lat, lon = trans_can$long),
data.frame(site = trans_glac$transect_id, lat = trans_glac$lat, lon = trans_glac$long)
)
locations_file <- file.path(daymet_dir, "locations.csv")
write.table(locations, locations_file, sep = ",", col.names = TRUE, row.names = FALSE, quote = FALSE)
download_daymet_batch(
file_location = locations_file,
start = 1980,
end = 2022,
internal = FALSE,
silent = FALSE,
path = daymet_dir
)Read Daymet data
daymet_glac_files <- file.path(daymet_dir, paste0(trans_glac$transect_id, "_1980_2022.csv"))
daymet_glac_data <- bind_rows(lapply(daymet_glac_files, function(x) {
read_daymet(site = basename(x), x) %>%
pivot_wider(names_from = measurement, values_from = value) %>%
drop_na("tmax..deg.c.", "tmin..deg.c.", "vp..Pa.") %>%
mutate(
site = str_split_i(site, pattern = "_1980_2022.csv", i = 1),
tavg = (tmax..deg.c. + tmin..deg.c.) / 2,
svp = 610.8 * exp((17.27 * tavg / (237.3 + tavg))),
rh = vp..Pa. / svp,
date = ymd(parse_date_time(paste(year, yday), orders = "yj"))
)
}))
daymet_can_files <- file.path(daymet_dir, paste0(trans_can$transect_id, "_1980_2022.csv"))
daymet_can_files <- daymet_can_files[file.exists(daymet_can_files)] ## Some transects had NA for lat/lon and were not able to be downloaded from Daymet
daymet_can_data <- bind_rows(lapply(daymet_can_files, function(x) {
read_daymet(site = basename(x), x) %>%
pivot_wider(names_from = measurement, values_from = value) %>%
drop_na("tmax..deg.c.", "tmin..deg.c.", "vp..Pa.") %>%
mutate(
site = str_split_i(site, pattern = "_1980_2022.csv", i = 1),
tavg = (tmax..deg.c. + tmin..deg.c.) / 2,
svp = 610.8 * exp((17.27 * tavg / (237.3 + tavg))),
rh = vp..Pa. / svp,
date = ymd(parse_date_time(paste(year, yday), orders = "yj"))
) })) %>%
rename(
tavg.daymet = tavg,
tmin.daymet = tmin..deg.c.,
tmax.daymet = tmax..deg.c.,
rh.daymet = rh
) %>%
mutate(year = year(date), month = month(date)) %>%
select(site, latitude, longitude, date, year, month, tavg.daymet, tmin.daymet, tmax.daymet, rh.daymet) %>%
drop_na(tavg.daymet, tmin.daymet, tmax.daymet, rh.daymet)
daymet_can_data$month <- factor(daymet_can_data$month)
## Combine gridMET and Daymet data for GLAC
all_glac_data <- left_join(daymet_glac_data, gridmet_glac_data, by = c("site" = "transect_id", "date" = "Date")) %>%
rename(
tmax.daymet = tmax..deg.c.,
tmin.daymet = tmin..deg.c.,
tavg.daymet = tavg,
rh.daymet = rh,
tmax.gridmet = Tmax,
tmin.gridmet = Tmin,
tavg.gridmet = T,
rh.gridmet = RH
) %>%
mutate(year = year(date), month = month(date)) %>%
select(site, latitude, longitude, date, year, month, tmax.daymet, tmin.daymet, tavg.daymet, rh.daymet, tmax.gridmet, tmin.gridmet, tavg.gridmet, rh.gridmet)
all_glac_data$month <- factor(all_glac_data$month)Modeling
par(mfrow = c(2, 2))
lm_min_temp <- lm(tmin.gridmet ~ tmin.daymet, data = all_glac_data)
summary(lm_min_temp)##
## Call:
## lm(formula = tmin.gridmet ~ tmin.daymet, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.0043 -1.8562 -0.1239 1.6496 26.4365
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6978149 0.0029937 -233.1 <2e-16 ***
## tmin.daymet 0.8589173 0.0003152 2724.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.853 on 1035868 degrees of freedom
## Multiple R-squared: 0.8775, Adjusted R-squared: 0.8775
## F-statistic: 7.423e+06 on 1 and 1035868 DF, p-value: < 2.2e-16plot(lm_min_temp)
lm_max_temp <- lm(tmax.gridmet ~ tmax.daymet, data = all_glac_data)
summary(lm_max_temp)##
## Call:
## lm(formula = tmax.gridmet ~ tmax.daymet, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7162 -1.7412 -0.0149 1.7015 26.3350
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2207224 0.0033361 66.16 <2e-16 ***
## tmax.daymet 0.9557934 0.0002818 3391.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.853 on 1035868 degrees of freedom
## Multiple R-squared: 0.9174, Adjusted R-squared: 0.9174
## F-statistic: 1.15e+07 on 1 and 1035868 DF, p-value: < 2.2e-16plot(lm_max_temp)
## lm_avg_temp <- lm(tavg.gridmet ~ tavg.daymet, data = all_glac_data)
## summary(lm_avg_temp)
## plot(lm_avg_temp)
lm_min_temp_month <- lm(tmin.gridmet ~ tmin.daymet * month, data = all_glac_data)
summary(lm_min_temp_month)##
## Call:
## lm(formula = tmin.gridmet ~ tmin.daymet * month, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7225 -1.6024 -0.0915 1.4724 22.4302
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.2993994 0.0137666 -239.668 < 2e-16 ***
## tmin.daymet 0.6839564 0.0009677 706.769 < 2e-16 ***
## month2 0.3250893 0.0201766 16.112 < 2e-16 ***
## month3 0.4343387 0.0191807 22.645 < 2e-16 ***
## month4 1.5118811 0.0181128 83.470 < 2e-16 ***
## month5 2.9899284 0.0162559 183.929 < 2e-16 ***
## month6 4.0362128 0.0177660 227.187 < 2e-16 ***
## month7 5.1949485 0.0218819 237.408 < 2e-16 ***
## month8 5.0407824 0.0208189 242.126 < 2e-16 ***
## month9 3.4088815 0.0168814 201.931 < 2e-16 ***
## month10 1.8625405 0.0165258 112.705 < 2e-16 ***
## month11 0.7477385 0.0183900 40.660 < 2e-16 ***
## month12 -0.3038716 0.0197175 -15.411 < 2e-16 ***
## tmin.daymet:month2 0.0503458 0.0014030 35.885 < 2e-16 ***
## tmin.daymet:month3 -0.0130553 0.0015533 -8.405 < 2e-16 ***
## tmin.daymet:month4 0.0048138 0.0019194 2.508 0.01214 *
## tmin.daymet:month5 0.0658253 0.0025548 25.765 < 2e-16 ***
## tmin.daymet:month6 0.0211843 0.0027521 7.697 1.39e-14 ***
## tmin.daymet:month7 -0.0066317 0.0025872 -2.563 0.01037 *
## tmin.daymet:month8 0.0080779 0.0025078 3.221 0.00128 **
## tmin.daymet:month9 0.0932186 0.0022252 41.891 < 2e-16 ***
## tmin.daymet:month10 0.0903253 0.0018157 49.748 < 2e-16 ***
## tmin.daymet:month11 0.0692842 0.0015276 45.354 < 2e-16 ***
## tmin.daymet:month12 0.0219103 0.0013725 15.964 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 1035846 degrees of freedom
## Multiple R-squared: 0.9033, Adjusted R-squared: 0.9033
## F-statistic: 4.208e+05 on 23 and 1035846 DF, p-value: < 2.2e-16plot(lm_min_temp_month)
lm_max_temp_month <- lm(tmax.gridmet ~ tmax.daymet * month, data = all_glac_data)
summary(lm_max_temp_month)##
## Call:
## lm(formula = tmax.gridmet ~ tmax.daymet * month, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.9459 -1.5625 0.0389 1.5955 23.1310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.366897 0.009456 -144.551 < 2e-16 ***
## tmax.daymet 0.766908 0.001408 544.868 < 2e-16 ***
## month2 0.490220 0.013369 36.669 < 2e-16 ***
## month3 1.703120 0.012894 132.089 < 2e-16 ***
## month4 2.129542 0.015515 137.258 < 2e-16 ***
## month5 3.493237 0.019382 180.233 < 2e-16 ***
## month6 5.400492 0.024368 221.619 < 2e-16 ***
## month7 6.581468 0.032846 200.375 < 2e-16 ***
## month8 6.639652 0.030715 216.169 < 2e-16 ***
## month9 4.465006 0.020605 216.691 < 2e-16 ***
## month10 2.336446 0.014927 156.520 < 2e-16 ***
## month11 0.532526 0.012841 41.472 < 2e-16 ***
## month12 -0.400633 0.013814 -29.002 < 2e-16 ***
## tmax.daymet:month2 0.027744 0.001956 14.187 < 2e-16 ***
## tmax.daymet:month3 -0.014213 0.002103 -6.758 1.4e-11 ***
## tmax.daymet:month4 0.076424 0.002153 35.505 < 2e-16 ***
## tmax.daymet:month5 0.052115 0.002049 25.429 < 2e-16 ***
## tmax.daymet:month6 -0.024925 0.002061 -12.096 < 2e-16 ***
## tmax.daymet:month7 -0.017652 0.002136 -8.263 < 2e-16 ***
## tmax.daymet:month8 -0.036003 0.002066 -17.423 < 2e-16 ***
## tmax.daymet:month9 0.028926 0.001899 15.233 < 2e-16 ***
## tmax.daymet:month10 0.032347 0.001923 16.821 < 2e-16 ***
## tmax.daymet:month11 0.061167 0.002013 30.390 < 2e-16 ***
## tmax.daymet:month12 -0.004926 0.001908 -2.581 0.00984 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.524 on 1035846 degrees of freedom
## Multiple R-squared: 0.9354, Adjusted R-squared: 0.9354
## F-statistic: 6.517e+05 on 23 and 1035846 DF, p-value: < 2.2e-16plot(lm_max_temp_month)
## lm_avg_temp_month <- lm(tavg.gridmet ~ tavg.daymet * month, data = all_glac_data)
## summary(lm_avg_temp_month)
## plot(lm_avg_temp_month)
lm_min_temp_month_lat <- lm(tmin.gridmet ~ tmin.daymet * month + latitude, data = all_glac_data)
summary(lm_min_temp_month_lat)##
## Call:
## lm(formula = tmin.gridmet ~ tmin.daymet * month + latitude, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.9290 -1.6025 -0.1015 1.4641 22.5613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.1532742 0.7032172 54.255 < 2e-16 ***
## tmin.daymet 0.6830351 0.0009662 706.907 < 2e-16 ***
## month2 0.3235017 0.0201429 16.060 < 2e-16 ***
## month3 0.4295731 0.0191488 22.433 < 2e-16 ***
## month4 1.5077159 0.0180826 83.379 < 2e-16 ***
## month5 2.9966756 0.0162291 184.648 < 2e-16 ***
## month6 4.0707330 0.0177460 229.389 < 2e-16 ***
## month7 5.2579999 0.0218715 240.405 < 2e-16 ***
## month8 5.0914141 0.0208018 244.759 < 2e-16 ***
## month9 3.4275406 0.0168561 203.341 < 2e-16 ***
## month10 1.8695061 0.0164986 113.313 < 2e-16 ***
## month11 0.7485102 0.0183593 40.770 < 2e-16 ***
## month12 -0.3047672 0.0196845 -15.483 < 2e-16 ***
## latitude -0.8523055 0.0144560 -58.958 < 2e-16 ***
## tmin.daymet:month2 0.0502393 0.0014006 35.869 < 2e-16 ***
## tmin.daymet:month3 -0.0139145 0.0015508 -8.973 < 2e-16 ***
## tmin.daymet:month4 0.0027212 0.0019166 1.420 0.156
## tmin.daymet:month5 0.0603215 0.0025523 23.634 < 2e-16 ***
## tmin.daymet:month6 0.0133264 0.0027508 4.845 1.27e-06 ***
## tmin.daymet:month7 -0.0143220 0.0025862 -5.538 3.06e-08 ***
## tmin.daymet:month8 0.0018559 0.0025058 0.741 0.459
## tmin.daymet:month9 0.0903806 0.0022220 40.675 < 2e-16 ***
## tmin.daymet:month10 0.0896805 0.0018126 49.475 < 2e-16 ***
## tmin.daymet:month11 0.0688926 0.0015251 45.173 < 2e-16 ***
## tmin.daymet:month12 0.0218616 0.0013702 15.955 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.531 on 1035845 degrees of freedom
## Multiple R-squared: 0.9037, Adjusted R-squared: 0.9036
## F-statistic: 4.048e+05 on 24 and 1035845 DF, p-value: < 2.2e-16plot(lm_min_temp_month_lat)
lm_max_temp_month_lat <- lm(tmax.gridmet ~ tmax.daymet * month + latitude, data = all_glac_data)
summary(lm_max_temp_month_lat)##
## Call:
## lm(formula = tmax.gridmet ~ tmax.daymet * month + latitude, data = all_glac_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.7418 -1.5528 0.0288 1.5806 23.3243
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.168172 0.695211 82.231 < 2e-16 ***
## tmax.daymet 0.765174 0.001403 545.434 < 2e-16 ***
## month2 0.492179 0.013323 36.942 < 2e-16 ***
## month3 1.711962 0.012850 133.223 < 2e-16 ***
## month4 2.147873 0.015464 138.899 < 2e-16 ***
## month5 3.525544 0.019320 182.485 < 2e-16 ***
## month6 5.455862 0.024294 224.574 < 2e-16 ***
## month7 6.667490 0.032750 203.588 < 2e-16 ***
## month8 6.700238 0.030619 218.826 < 2e-16 ***
## month9 4.494927 0.020538 218.855 < 2e-16 ***
## month10 2.351215 0.014878 158.037 < 2e-16 ***
## month11 0.536441 0.012797 41.919 < 2e-16 ***
## month12 -0.402464 0.013767 -29.234 < 2e-16 ***
## latitude -1.203342 0.014291 -84.205 < 2e-16 ***
## tmax.daymet:month2 0.028128 0.001949 14.433 < 2e-16 ***
## tmax.daymet:month3 -0.015274 0.002096 -7.287 3.16e-13 ***
## tmax.daymet:month4 0.075691 0.002145 35.284 < 2e-16 ***
## tmax.daymet:month5 0.051076 0.002043 25.007 < 2e-16 ***
## tmax.daymet:month6 -0.026842 0.002054 -13.070 < 2e-16 ***
## tmax.daymet:month7 -0.020213 0.002129 -9.493 < 2e-16 ***
## tmax.daymet:month8 -0.037272 0.002059 -18.098 < 2e-16 ***
## tmax.daymet:month9 0.028698 0.001892 15.165 < 2e-16 ***
## tmax.daymet:month10 0.032456 0.001916 16.935 < 2e-16 ***
## tmax.daymet:month11 0.060794 0.002006 30.309 < 2e-16 ***
## tmax.daymet:month12 -0.004863 0.001902 -2.557 0.0106 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.515 on 1035845 degrees of freedom
## Multiple R-squared: 0.9358, Adjusted R-squared: 0.9358
## F-statistic: 6.291e+05 on 24 and 1035845 DF, p-value: < 2.2e-16plot(lm_max_temp_month_lat)
## lm_avg_temp_month_lat <- lm(tavg.gridmet ~ tavg.daymet * month + latitude, data = all_glac_data)
## summary(lm_avg_temp_month_lat)
## plot(lm_avg_temp_month_lat)
anova(lm_min_temp_month, lm_min_temp)## Analysis of Variance Table
##
## Model 1: tmin.gridmet ~ tmin.daymet * month
## Model 2: tmin.gridmet ~ tmin.daymet
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1035846 6657236
## 2 1035868 8432644 -22 -1775408 12557 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm_min_temp_month_lat, lm_min_temp_month)## Analysis of Variance Table
##
## Model 1: tmin.gridmet ~ tmin.daymet * month + latitude
## Model 2: tmin.gridmet ~ tmin.daymet * month
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1035845 6634970
## 2 1035846 6657236 -1 -22266 3476.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm_max_temp_month, lm_max_temp)## Analysis of Variance Table
##
## Model 1: tmax.gridmet ~ tmax.daymet * month
## Model 2: tmax.gridmet ~ tmax.daymet
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1035846 6599163
## 2 1035868 8434053 -22 -1834891 13092 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lm_max_temp_month_lat, lm_max_temp_month)## Analysis of Variance Table
##
## Model 1: tmax.gridmet ~ tmax.daymet * month + latitude
## Model 2: tmax.gridmet ~ tmax.daymet * month
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1035845 6554297
## 2 1035846 6599163 -1 -44865 7090.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Models are improved by adding month and lat, from results of ESS F-test
## Add predictions to dataframes
all_glac_data <- all_glac_data %>% add_predictions(lm_min_temp_month_lat, var = "tmin.gridmet.pred") %>%
add_predictions(lm_max_temp_month_lat, var = "tmax.gridmet.pred")
daymet_can_data <- daymet_can_data %>% add_predictions(lm_min_temp_month_lat, var = "tmin.gridmet.pred") %>%
add_predictions(lm_max_temp_month_lat, var = "tmax.gridmet.pred")Plots
Visualize historical and predicted data for one site for each GLAC and CAN, for months of August and September, 2000.
all_glac_data %>%
pivot_longer(cols = c(tmax.daymet, tmin.daymet, tmax.gridmet, tmin.gridmet, tmin.gridmet.pred, tmax.gridmet.pred)) %>%
filter(year(date) == 2000, month(date) == 8 | month(date) == 9, site == "WB0401") %>%
ggplot() +
geom_line(mapping = aes(x = date, y = value, color = name)) +
ggtitle("GLAC historical Daymet and gridMET data with predicted gridMET")
daymet_can_data %>%
pivot_longer(cols = c(tmax.daymet, tmin.daymet, tmin.gridmet.pred, tmax.gridmet.pred)) %>%
filter(year(date) == 2020, month(date) == 8 | month(date) == 9, site == "JA-03-050") %>%
ggplot() +
geom_line(mapping = aes(x = date, y = value, color = name)) +
ggtitle("CAN historical Daymet data with predicted gridMET")
Reconstruct Tavg data from modelled Tmin and Tmax
all_glac_data <- all_glac_data %>%
mutate(tavg.gridmet.pred = (tmin.gridmet.pred + tmax.gridmet.pred) / 2)
daymet_can_data <- daymet_can_data %>%
mutate(tavg.gridmet.pred = (tmin.gridmet.pred + tmax.gridmet.pred) / 2)
all_glac_data %>%
pivot_longer(cols = c(tavg.daymet, tavg.gridmet, tavg.gridmet.pred)) %>%
filter(year(date) == 2000, month(date) == 8 | month(date) == 9, site == "WB0401") %>%
ggplot() +
geom_line(mapping = aes(x = date, y = value, color = name)) +
ggtitle("Reconstructed gridMET Tavg data vs historical gridMET for GLAC")
daymet_can_data %>%
pivot_longer(cols = c(tavg.daymet, tavg.gridmet.pred)) %>%
filter(year(date) == 2000, month(date) == 8 | month(date) == 9, site == "JA-03-050") %>%
ggplot() +
geom_line(mapping = aes(x = date, y = value, color = name)) +
ggtitle("Reconstructed gridMET Tavg data vs historical gridMET for GLAC")
R Version
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 12 (bookworm)
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.21.so; LAPACK version 3.11.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
## [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/Boise
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] modelr_0.1.11 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
## [9] ggplot2_3.4.3 tidyverse_2.0.0 lubridate_1.9.2 here_1.0.1 daymetr_1.7
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 utf8_1.2.3 generics_0.1.3 stringi_1.7.12 hms_1.1.3 digest_0.6.33 magrittr_2.0.3
## [8] evaluate_0.21 grid_4.3.1 timechange_0.2.0 fastmap_1.1.1 jsonlite_1.8.7 rprojroot_2.0.3 backports_1.4.1
## [15] fansi_1.0.4 scales_1.2.1 codetools_0.2-19 jquerylib_0.1.4 cli_3.6.1 rlang_1.1.1 crayon_1.5.2
## [22] bit64_4.0.5 munsell_0.5.0 cachem_1.0.8 yaml_2.3.7 withr_2.5.0 tools_4.3.1 parallel_4.3.1
## [29] tzdb_0.4.0 colorspace_2.1-0 ncdf4_1.21 broom_1.0.5 vctrs_0.6.3 R6_2.5.1 lifecycle_1.0.3
## [36] fs_1.6.3 bit_4.0.5 vroom_1.6.3 pkgconfig_2.0.3 bslib_0.5.1 pillar_1.9.0 gtable_0.3.4
## [43] glue_1.6.2 xfun_0.40 tidyselect_1.2.0 knitr_1.44 farver_2.1.1 htmltools_0.5.6 rmarkdown_2.24
## [50] labeling_0.4.3 compiler_4.3.1References
Abatzoglou, John T. 2013. “Development of Gridded Surface
Meteorological Data for Ecological Applications and Modelling.”
International Journal of Climatology 33 (1): 121–31. https://doi.org/10.1002/joc.3413.
Abatzoglou, John T., and Timothy J. Brown. 2012. “A Comparison of
Statistical Downscaling Methods Suited for Wildfire
Applications.” International Journal of Climatology 32
(5): 772–80. https://doi.org/10.1002/joc.2312.
Jamison, Elizabeth-Ann, David Thoma, Erin Shanahan, and Kristin Legg. In
Preparation. “Climatic Vulnerability of Whitebark Pine to White
Pine Blister Rust Infection.”
Tercek, Michael T., David Thoma, John E. Gross, Kirk Sherrill, Stefanie
Kagone, and Gabriel Senay. 2021. “Historical Changes in Plant
Water Use and Need in the Continental United
States.” PLOS ONE 16 (9): e0256586. https://doi.org/10.1371/journal.pone.0256586.
Thornton, P. E., M. M. Thornton, B. W. Mayer, N. Wilhelmi, Y. Wei, R.
Devarakonda, and R. B. Cook. 2014. “Daymet: Daily Surface
Weather Data on a 1-Km Grid for North
America, Version 2.” ORNL DAAC, May.
https://doi.org/10.3334/ORNLDAAC/1219.