In mt-climate-office/mcor: The Core Montana Climate Office R Package
data_out <- "../vignettes/data/raw-data/"
huc <- 8
options(knitr.table.format = "html")
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
cache = FALSE,
out.width = "100%"
)
library(mcor)
library(magrittr)
library(raster)
library(tidyverse)
library(sf)
library(mapview)
library(kableExtra)
mapviewOptions(basemaps = c("CartoDB.Positron"))
The goal of this vignette is to document how the MCO makes comparisons between historic, "normal" (1981--2010 and 1971-2000) SWE data from the NRCS SNOTEL network and the VIC hydroclimatology models derived from monthly BCSD downscaled climate models. We compare between 1981--2010 and 1971--2000 normal instrumental April 1 SWE values (from SNOTEL) and modeled mid-century (2040--2069) values from the BCSD VIC model. This is a point reconstruction that may be aggregated to hydrological basins in the style of NRCS SWE maps.
Summary {-}
In this script, we will:
- Download the 2040--2069 mid-century monthly BCSD VIC data for Montana and the HUC
r huc hydrological basins it overlaps.
- Download the 1981--2010 and 1971--2000 monthly VIC data for that same area.
- Download the 1981--2010 and 1971--2000 normal (median) instrumental SNOTEL data from within that same area.
- Extract the RCP 4.5 BCSD VIC data from the locations of the instrumental SNOTEL data for comparison.
- Generate maps of percent of normal April 1 SWE.
Downloading and loading monthly BCSD VIC data for Montana {-}
The first step is to download the monthly RCP 4.5 BCSD VIC data for Montana. These data are stored in several places; here, we access them through a THREDDS server maintained by the USGS Water Resources Mission Area (formerly the USGS Center for Integrated Data Analytics, CIDA). In this section, we:
- Download the data covering the area of all HUC
r huc watershed that intersect with the state of Montana.
- Extract out only April 1 SWE values.
- Convert the data from millimeters to inches.
- Calculate the quartiles of the 31 models for each year.
- Calculate the median quartile values over the mid-century (2040--2069), 1981--2010, and 1971--2000 periods.
dir.create(data_out,
showWarnings = FALSE,
recursive = TRUE)
mm_to_in <- function(x){
x %>%
magrittr::multiply_by(0.0393701) # mm to inches
}
bcsd_swe_midcentury_vars <- thredds::tds_ncss_list_vars("https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html") %$%
name %>%
stringr::str_subset("_swe") %>%
stringr::str_subset("_rcp45_") %>%
sort()
bcsd_swe_midcentury <-
bcsd_swe_midcentury_vars %>%
purrr::map_chr(function(var){
thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html",
out_file = stringr::str_c(data_out,"/",var,"_2040-2069.nc"),
bbox = mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_transform(4326) %>%
# magrittr::add(c(360,0)) %>%
sf::st_bbox(),
vars = var,
ncss_args = list(time_start = "2040-01-01",
time_end = "2069-12-31"),
overwrite = FALSE)
}) %>%
purrr::compact() %>%
magrittr::set_names(bcsd_swe_midcentury_vars) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(function(x){
x %>%
raster::subset(
x %>%
names() %>%
gsub("X","",.) %>%
as.Date(format = "%Y.%m.%d") %>%
lubridate::month() %in%
c(4) %>%
which()
) %>%
mm_to_in() %>%
raster::as.list()
}) %>%
# The next eight lines first calculates the quantiles of the 31 models for each year
# then the median for each percentile across the mid-century
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(raster::calc, quantile, na.rm = T) %>%
purrr::map(raster::as.list) %>%
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>%
purrr::map(raster::calc, median, na.rm = TRUE) %>%
raster::stack(quick = TRUE)
bcsd_swe_1971_2000 <-
bcsd_swe_midcentury_vars %>%
purrr::map_chr(function(var){
thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html",
out_file = stringr::str_c(data_out,"/",var,"_1971-2000.nc"),
bbox = mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_transform(4326) %>%
# magrittr::add(c(360,0)) %>%
sf::st_bbox(),
vars = var,
ncss_args = list(time_start = "1971-01-01",
time_end = "2000-12-31"),
overwrite = FALSE)
}) %>%
purrr::compact() %>%
magrittr::set_names(bcsd_swe_midcentury_vars) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(function(x){
x %>%
raster::subset(
x %>%
names() %>%
gsub("X","",.) %>%
as.Date(format = "%Y.%m.%d") %>%
lubridate::month() %in%
c(4) %>%
which()
) %>%
mm_to_in() %>%
raster::as.list()
}) %>%
# The next eight lines first calculates the quantiles of the 31 models for each year
# then the median for each percentile across the mid-century
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(raster::calc, quantile, na.rm = T) %>%
purrr::map(raster::as.list) %>%
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>%
purrr::map(raster::calc, median, na.rm = TRUE) %>%
raster::stack(quick = TRUE)
bcsd_swe_1981_2010 <-
bcsd_swe_midcentury_vars %>%
purrr::map_chr(function(var){
thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html",
out_file = stringr::str_c(data_out,"/",var,"_1981-2010.nc"),
bbox = mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_transform(4326) %>%
# magrittr::add(c(360,0)) %>%
sf::st_bbox(),
vars = var,
ncss_args = list(time_start = "1981-01-01",
time_end = "2010-12-31"),
overwrite = FALSE)
}) %>%
purrr::compact() %>%
magrittr::set_names(bcsd_swe_midcentury_vars) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(function(x){
x %>%
raster::subset(
x %>%
names() %>%
gsub("X","",.) %>%
as.Date(format = "%Y.%m.%d") %>%
lubridate::month() %in%
c(4) %>%
which()
) %>%
mm_to_in() %>%
raster::as.list()
}) %>%
# The next eight lines first calculates the quantiles of the 31 models for each year
# then the median for each percentile across the mid-century
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(raster::calc, quantile, na.rm = T) %>%
purrr::map(raster::as.list) %>%
purrr::transpose() %>%
purrr::map(raster::stack, quick = TRUE) %>%
magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>%
purrr::map(raster::calc, median, na.rm = TRUE) %>%
raster::stack(quick = TRUE)
Here is a panel of maps of the 1981--2010 median April 1 SWE (in.) compared to the mid-century median April 1 SWE (in.), with HUC r huc watersheds overlaid.
basins <- mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_transform(4326)
# pal = mapview::mapviewPalette("mapviewTopoColors")
pal <- colorRampPalette(rev(RColorBrewer::brewer.pal(9,"Blues")), alpha = TRUE)
m1 <- mapview::mapview(bcsd_swe_1981_2010[["X0.5."]],
col.regions = pal(100),
at = seq(0, 60, 5),
legend = TRUE,
alpha = 0.5,
na.color = NA,
layer.name = "1981-2010 SWE (in.)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "white")
m2 <- mapview::mapview(bcsd_swe_midcentury[["X0.5."]],
col.regions = pal(100),
at = seq(0, 60, 5),
legend = TRUE,
alpha = 0.5,
na.color = NA,
layer.name = "Mid-century SWE (in.)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "white")
mapview::sync(list(m1,
m2),
ncol = 1)
Downloading the 1971--2000 and 1981--2010 normal instrumental SNOTEL data {-}
Our instrumental data comes from the NRCS SNOTEL network, available through the NRCS National Water and Climate Center. The NRCS publishes two datasets: an inventory of SNOTEL sites, including their associated HUC basins, and the data recorded from the sites. It is important to remember that one SNOTEL site can be associated with multiple basins; the SNOTEL sites are often near basin boundaries, and when aggregating SNOTEL data to great basin reports, the NRCS uses both SNOTEL stations within the basin, and ones that are nearby. We also download both 1971--2000 and 1981--2010 normal instrumental SNOTEL data. The MCO has been using 1981--2010 as its current standard for defining normal climatology, but the water chapter of the MCA used 1971--2000. This section of code:
- Downloads the SNOTEL inventory for SNOTEL sites within Montana HUC
r huc watersheds.
- Downloads the normal data for those SNOTEL sites.
- Joins those two datasets, resulting in spatially-referenced SNOTEL normals data.
date <- "2018-04-01" %>%
lubridate::as_date()
snotel_inventory <- mcor::mco_get_snotel_inventory() %>%
dplyr::left_join(mcor::mco_get_snotel_huc() %>%
dplyr::mutate(`WBD code` = stringr::str_sub(`WBD code`,1,huc)) %>%
dplyr::distinct()) %>%
dplyr::filter(`WBD code` %in% (mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` %in% c(huc)) %$%
`WBD code`)) %>%
sf::st_as_sf() %>%
sf::st_transform(mcor::mt_state_plane) %>%
dplyr::select(-`Station Id`:-`End Date`)
snotel_data <- mcor::mco_get_snotel_data(stations = snotel_inventory$Station %>%
unique(),
variables = c('WTEQ::normal_1981',
'WTEQ::normal_1971'),
start_date = date,
end_date = date)
snotel_data <- snotel_inventory %>%
dplyr::left_join(snotel_data,
by = c("Station")) %>%
dplyr::arrange(Station) %>%
stats::na.omit() %>%
dplyr::select(Station,
`WBD code`,
`Normal Snow Water Equivalent (1981-2010) (in) Start of Day Values`,
`Normal Snow Water Equivalent (1971-2000) (in) Start of Day Values`) %>%
dplyr::rename(`SNOTEL 1981-2010 normal (in.)` = `Normal Snow Water Equivalent (1981-2010) (in) Start of Day Values`,
`SNOTEL 1971-2000 normal (in.)` = `Normal Snow Water Equivalent (1971-2000) (in) Start of Day Values`)
Here is a panel of maps of the 1981--2010 and 1971--2000 April 1 SWE values.
snotel_data_slim <-
snotel_data %>%
dplyr::select(-`WBD code`) %>%
dplyr::distinct() %>%
dplyr::rename(`Normal_1981` = "SNOTEL 1981-2010 normal (in.)",
`Normal_1971` = "SNOTEL 1971-2000 normal (in.)")
m1 <- mapview::mapview(snotel_data_slim,
zcol = "Normal_1981",
col.regions = pal(100),
at = seq(0, 60, 5),
legend = TRUE,
alpha = 0.5,
na.color = NA,
layer.name = "1981-2010 Norm. (in.)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "black")
m2 <- mapview::mapview(snotel_data_slim,
zcol = "Normal_1971",
col.regions = pal(100),
at = seq(0, 60, 5),
legend = TRUE,
alpha = 0.5,
na.color = NA,
layer.name = "1971-2000 Norm. (in.)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "black")
mapview::sync(list(m1,
m2),
ncol = 1)
Extract the BCSD VIC data from the SNOTEL locations {-}
This section of code extracts the BCSD mid-century data for the SNOTEL locations, and attaches those data to the SNOTEL data downloaded above. Notice that the BCSD values for the extracted SNOTEL locations are very different from their measured values.
snotel_bcsd_midcentury <- raster::extract(bcsd_swe_midcentury,
snotel_data %>%
sf::st_cast() %>%
sf::st_transform(4326),
df = TRUE) %>%
dplyr::select(-ID) %>%
magrittr::set_names(c("BCSD Midcentury min.",
"BCSD Midcentury 25%",
"BCSD 2049-2060 median (in.)",
"BCSD Midcentury 75%",
"BCSD Midcentury max.")) %>%
dplyr::select(`BCSD 2049-2060 median (in.)`)
snotel_bcsd_normal_1971_2000 <- raster::extract(bcsd_swe_1971_2000,
snotel_data %>%
sf::st_cast() %>%
sf::st_transform(4326),
df = TRUE) %>%
dplyr::select(-ID) %>%
magrittr::set_names(c("BCSD 1971-2000 min.",
"BCSD 1971-2000 25%",
"BCSD 1971-2000 normal (in.)",
"BCSD 1971-2000 75%",
"BCSD 1971-2000 max.")) %>%
dplyr::select(`BCSD 1971-2000 normal (in.)`)
snotel_bcsd_normal_1981_2010 <- raster::extract(bcsd_swe_1981_2010,
snotel_data %>%
sf::st_cast() %>%
sf::st_transform(4326),
df = TRUE) %>%
dplyr::select(-ID) %>%
magrittr::set_names(c("BCSD 1981-2010 min.",
"BCSD 1981-2010 25%",
"BCSD 1981-2010 normal (in.)",
"BCSD 1981-2010 75%",
"BCSD 1981-2010 max.")) %>%
dplyr::select(`BCSD 1981-2010 normal (in.)`)
snotel_data %<>%
dplyr::bind_cols(snotel_bcsd_midcentury) %>%
dplyr::bind_cols(snotel_bcsd_normal_1971_2000) %>%
dplyr::bind_cols(snotel_bcsd_normal_1981_2010) %>%
dplyr::select(Station,
`WBD code`,
`SNOTEL 1981-2010 normal (in.)`,
`BCSD 1981-2010 normal (in.)`,
`SNOTEL 1971-2000 normal (in.)`,
`BCSD 1971-2000 normal (in.)`,
`BCSD 2049-2060 median (in.)`
)
snotel_data %>%
sf::st_set_geometry(NULL) %>%
dplyr::select(-`WBD code`) %>%
dplyr::distinct() %>%
knitr::kable(digits = 2) %>%
kable_styling() %>%
scroll_box(width = "100%", height = "500px")
Mapping percent of normal April 1 SWE {-}
Here, we calculate the % of normal SWE, and make some maps relating future projections to normal conditions.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9]))
snotel_data_1981 <-
snotel_data %>%
dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1981-2010 normal (in.)`) %>%
dplyr::distinct()
snotel_data_1971 <-
snotel_data %>%
dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1971-2000 normal (in.)`) %>%
dplyr::distinct()
m1 <- mapview::mapview(snotel_data_1981,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1981-2010)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "black")
m2 <- mapview::mapview(snotel_data_1971,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1971-2000)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "black")
mapview::sync(list(m1,
m2),
ncol = 1)
These maps make the same comparison among the gridded modeled data.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9]))
# bcsd_swe_midcentury_nozero <- bcsd_swe_midcentury
# bcsd_swe_1981_2010_nozero <- bcsd_swe_1981_2010
# bcsd_swe_1971_2000_nozero <- bcsd_swe_1971_2000
#
# bcsd_swe_midcentury_nozero[bcsd_swe_midcentury == 0] <- 0.0001
# bcsd_swe_1981_2010_nozero[bcsd_swe_1981_2010 == 0] <- 0.0001
# bcsd_swe_1971_2000_nozero[bcsd_swe_1971_2000 == 0] <- 0.0001
bcsd_swe_percent_normal_1981 <- 100 * bcsd_swe_midcentury / bcsd_swe_1981_2010
bcsd_swe_percent_normal_1971 <- 100 * bcsd_swe_midcentury / bcsd_swe_1971_2000
# bcsd_swe_percent_normal_1981[bcsd_swe_percent_normal_1981 > 200] <- NA
# bcsd_swe_percent_normal_1971[bcsd_swe_percent_normal_1971 > 200] <- NA
m1 <- mapview::mapview(bcsd_swe_percent_normal_1981[["X0.5."]],
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1981-2010)") +
mapview::mapview(basins,
alpha.regions = 0,
label = basins$Watershed,
color = "black")
m2 <- mapview::mapview(bcsd_swe_percent_normal_1971[["X0.5."]] %>%
round(),
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1971-2000)") +
mapview::mapview(basins,
alpha.regions = 0,
label = FALSE,
color = "black")
mapview::sync(list(m1,
m2),
ncol = 1)
Finally, we to emulate the NRCS SWE maps. In the first set of maps, we aggregate modeled values over each HUC basin. This roughly matches what was done during the Montana Climate Assessment. In the second set, we aggregate stations within basins by taking their mean, following the procedures used by the NRCS to make their standard SWE maps.
Midcentury % of normal April 1st SWE, modeled to modeled, gridded values {-}
These maps compare modeled normal April 1 SWE values to modeled midcentury values. Gridded reconstructions are averaged over entire basins.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9]))
watersheds_cropped <-
mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc)# %>%
# sf::st_intersection(mt_state) %>%
# dplyr::filter(sf::st_area(.) > units::set_units(1000000000,m^2)) %>%
# dplyr::rename(`WBD code` = `WBD.code`,
# `Hydrologic Unit` = Hydrologic.Unit)
bcsd_swe_percent_normal_1981_basins <-
watersheds_cropped %>%
dplyr::bind_cols(bcsd_swe_percent_normal_1981 %>%
raster::extract(watersheds_cropped %>%
sf::st_transform(4326),
fun = mean,
na.rm=TRUE,
df = TRUE) %>%
dplyr::select(-ID) %>%
dplyr::rename(`Midcentury min.` = `X0.`,
`Midcentury 25%` = `X0.25.`,
`Midcentury median` = `X0.5.`,
`Midcentury 75%` = `X0.75.`,
`Midcentury max.` = `X1.`)) %>%
mutate_all(funs(replace(., is.na(.), NA))) %>%
dplyr::rename("Midcentury_median" = "Midcentury median") %>%
dplyr::select(Watershed, Midcentury_median) %>%
na.omit()
bcsd_swe_percent_normal_1971_basins <-
watersheds_cropped %>%
dplyr::bind_cols(bcsd_swe_percent_normal_1971 %>%
raster::extract(watersheds_cropped %>%
sf::st_transform(4326),
fun = mean,
na.rm=TRUE,
df = TRUE) %>%
dplyr::select(-ID) %>%
dplyr::rename(`Midcentury min.` = `X0.`,
`Midcentury 25%` = `X0.25.`,
`Midcentury median` = `X0.5.`,
`Midcentury 75%` = `X0.75.`,
`Midcentury max.` = `X1.`)) %>%
mutate_all(funs(replace(., is.na(.), NA))) %>%
dplyr::rename("Midcentury_median" = "Midcentury median") %>%
dplyr::select(Watershed, Midcentury_median) %>%
na.omit()
m1 <- mapview::mapview(bcsd_swe_percent_normal_1981_basins,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1981-2010), gridded",
label = paste0(bcsd_swe_percent_normal_1981_basins$Watershed,": ",bcsd_swe_percent_normal_1981_basins$Midcentury_median %>% round(),"%"))
m2 <- mapview::mapview(bcsd_swe_percent_normal_1971_basins,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1971-2000), gridded",
label = paste0(bcsd_swe_percent_normal_1971_basins$Watershed,": ",bcsd_swe_percent_normal_1971_basins$Midcentury_median %>% round(),"%"))
mapview::sync(list(m1,
m2),
ncol = 1)
We can also look at the three basins in figure 3-10 of the Montana Climate Assessment, in an effort to reproduce that figure. Here, we calculate the basin-wide boxplots for three basins roughly conforming to the basins in the MCA, using the 1971--2000 normals; these don't conform precisely with figure 3-10, as the pour points for defining the basins are underspecified in the documentation for for Montana Climate Assessment. Here, we define the Clark Fork above St. Regis as including the Middle Clark Fork, Bitterroot, Flint-Rock, Blackfoot, and Upper Clark Fork HUC 8 basins; the Missouri above Toston as including the Gallatin, Madison, Jefferson, Ruby, Beaverhead, and Big Hole basins; and the Yellowstone above Billings as including the Upper Yellowstone, Stillwater, Clarks Fork Yellowstone, Shields, and Yellowstone Headwaters basins. This analysis includes points for each of 31 climate models available in the BCSD VIC dataset.
bcsd_swe_1971_2000_all_models <-
bcsd_swe_midcentury_vars %>%
purrr::map_chr(function(var){
thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html",
out_file = stringr::str_c(data_out,"/",var,"_1971-2000.nc"),
bbox = mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_transform(4326) %>%
# magrittr::add(c(360,0)) %>%
sf::st_bbox(),
vars = var,
ncss_args = list(time_start = "1971-01-01",
time_end = "2000-12-31"),
overwrite = FALSE)
}) %>%
purrr::compact() %>%
magrittr::set_names(bcsd_swe_midcentury_vars) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(function(x){
x %>%
raster::subset(
x %>%
names() %>%
gsub("X","",.) %>%
as.Date(format = "%Y.%m.%d") %>%
lubridate::month() %in%
c(4) %>%
which()
) %>%
mm_to_in() %>%
raster::as.list()
}) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(mean) %>%
raster::stack(quick = TRUE)
bcsd_swe_midcentury_all_models <-
bcsd_swe_midcentury_vars %>%
purrr::map_chr(function(var){
thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html",
out_file = stringr::str_c(data_out,"/",var,"_2040-2069.nc"),
bbox = mcor::mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_bbox() %>%
sf::st_as_sfc() %>%
sf::st_transform(4326) %>%
# magrittr::add(c(360,0)) %>%
sf::st_bbox(),
vars = var,
ncss_args = list(time_start = "2040-01-01",
time_end = "2069-12-31"),
overwrite = FALSE)
}) %>%
purrr::compact() %>%
magrittr::set_names(bcsd_swe_midcentury_vars) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(function(x){
x %>%
raster::subset(
x %>%
names() %>%
gsub("X","",.) %>%
as.Date(format = "%Y.%m.%d") %>%
lubridate::month() %in%
c(4) %>%
which()
) %>%
mm_to_in() %>%
raster::as.list()
}) %>%
purrr::map(raster::stack, quick = TRUE) %>%
purrr::map(mean) %>%
raster::stack(quick = TRUE)
watersheds_cropped <-
mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc) %>%
sf::st_transform(4326) %>%
dplyr::mutate(.,ID = 1:nrow(.))
(bcsd_swe_midcentury_all_models / bcsd_swe_1971_2000_all_models) %>%
raster::extract(watersheds_cropped,
df = TRUE) %>%
tidyr::gather(Model, Value, -ID) %>%
# dplyr::rename(`Midcentury median` = `X0.5.`) %>%
dplyr::select(ID, Value) %>%
dplyr::left_join(x = watersheds_cropped, y = .) %>%
dplyr::filter(!is.nan(Value),
!is.na(Value),
!is.infinite(Value)) %>%
dplyr::mutate(Watershed = factor(Watershed)) %>%
dplyr::mutate(Watershed = forcats::fct_collapse(Watershed,
`Clark Fork above St. Regis` = c("Middle Clark Fork",
"Bitterroot",
"Flint-Rock",
"Blackfoot",
"Upper Clark Fork"),
`Missouri above Toston` = c("Gallatin",
"Madison",
"Jefferson",
"Ruby",
"Beaverhead",
"Big Hole"),
`Yellowstone above Billings` = c("Upper Yellowstone",
"Stillwater",
"Clarks Fork Yellowstone",
"Shields",
"Yellowstone Headwaters"))) %>%
dplyr::filter(Watershed %in% c("Clark Fork above St. Regis",
"Missouri above Toston",
"Yellowstone above Billings"),
Value >= 0,
Value <= 2) %>%
dplyr::mutate(Watershed = forcats::fct_drop(Watershed),
Watershed = forcats::fct_relevel(Watershed,
"Clark Fork above St. Regis",
"Missouri above Toston",
"Yellowstone above Billings")) %>%
# sf::st_set_geometry(NULL) %>%
# dplyr::select(Watershed, `Midcentury median`) %>%
ggplot2::ggplot(ggplot2::aes(x = Watershed,
y = 100 * (Value - 1))) +
ggplot2::geom_boxplot(outlier.shape = NA) +
ylab("% Change in April 1 SWE")
Midcentury % of normal April 1st SWE, station to modeled, point values {-}
These maps compare empirical April 1 SWE values from SNOTEL stations to the modeled midcentury SWE values extracted from those same station locations (i.e., the data from the station table above).
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9]))
watersheds_cropped <-
mt_watersheds_simple %>%
dplyr::filter(`Hydrologic Unit` == huc)
snotel_data_1981 <-
snotel_data %>%
dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1981-2010 normal (in.)`) %>%
group_by(`WBD code`) %>%
dplyr::summarise(`Stations Count` = n(),
`Midcentury_median` = mean(Midcentury_median) %>%
round()) %>%
dplyr::filter(`Stations Count` >= 3) %>%
sf::st_set_geometry(NULL) %>%
dplyr::left_join(watersheds_cropped) %>%
sf::st_sf() %>%
na.omit() %>%
dplyr::select(-`Hydrologic Unit`) %>%
dplyr::select(`WBD code`, Watershed, dplyr::everything())
snotel_data_1971 <-
snotel_data %>%
dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1971-2000 normal (in.)`) %>%
group_by(`WBD code`) %>%
dplyr::summarise(`Stations Count` = n(),
`Midcentury_median` = mean(Midcentury_median) %>%
round()) %>%
dplyr::filter(`Stations Count` >= 3) %>%
sf::st_set_geometry(NULL) %>%
dplyr::left_join(watersheds_cropped) %>%
sf::st_sf() %>%
na.omit() %>%
dplyr::select(-`Hydrologic Unit`) %>%
dplyr::select(`WBD code`, Watershed, dplyr::everything())
m1 <- mapview::mapview(snotel_data_1981,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1981-2010)",
label = paste0(snotel_data_1981$Watershed,": ",snotel_data_1981$Midcentury_median,"%"))
m2 <- mapview::mapview(snotel_data_1971,
zcol = "Midcentury_median",
col.regions = pal(100),
at = seq(0, 200, 20),
legend = TRUE,
na.color = NA,
layer.name = "% Norm. (1971-2000)",
label = paste0(snotel_data_1971$Watershed,": ",snotel_data_1971$Midcentury_median,"%"))
mapview::sync(list(m1,
m2),
ncol = 1)
mt-climate-office/mcor documentation built on March 27, 2024, 6:30 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
data_out <- "../vignettes/data/raw-data/" huc <- 8 options(knitr.table.format = "html") knitr::opts_chunk$set( collapse = TRUE, comment = "#>", cache = FALSE, out.width = "100%" ) library(mcor) library(magrittr) library(raster) library(tidyverse) library(sf) library(mapview) library(kableExtra) mapviewOptions(basemaps = c("CartoDB.Positron"))
The goal of this vignette is to document how the MCO makes comparisons between historic, "normal" (1981--2010 and 1971-2000) SWE data from the NRCS SNOTEL network and the VIC hydroclimatology models derived from monthly BCSD downscaled climate models. We compare between 1981--2010 and 1971--2000 normal instrumental April 1 SWE values (from SNOTEL) and modeled mid-century (2040--2069) values from the BCSD VIC model. This is a point reconstruction that may be aggregated to hydrological basins in the style of NRCS SWE maps.
Summary {-}
In this script, we will:
- Download the 2040--2069 mid-century monthly BCSD VIC data for Montana and the HUC
r huchydrological basins it overlaps. - Download the 1981--2010 and 1971--2000 monthly VIC data for that same area.
- Download the 1981--2010 and 1971--2000 normal (median) instrumental SNOTEL data from within that same area.
- Extract the RCP 4.5 BCSD VIC data from the locations of the instrumental SNOTEL data for comparison.
- Generate maps of percent of normal April 1 SWE.
Downloading and loading monthly BCSD VIC data for Montana {-}
The first step is to download the monthly RCP 4.5 BCSD VIC data for Montana. These data are stored in several places; here, we access them through a THREDDS server maintained by the USGS Water Resources Mission Area (formerly the USGS Center for Integrated Data Analytics, CIDA). In this section, we:
- Download the data covering the area of all HUC
r hucwatershed that intersect with the state of Montana. - Extract out only April 1 SWE values.
- Convert the data from millimeters to inches.
- Calculate the quartiles of the 31 models for each year.
- Calculate the median quartile values over the mid-century (2040--2069), 1981--2010, and 1971--2000 periods.
dir.create(data_out, showWarnings = FALSE, recursive = TRUE) mm_to_in <- function(x){ x %>% magrittr::multiply_by(0.0393701) # mm to inches } bcsd_swe_midcentury_vars <- thredds::tds_ncss_list_vars("https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html") %$% name %>% stringr::str_subset("_swe") %>% stringr::str_subset("_rcp45_") %>% sort() bcsd_swe_midcentury <- bcsd_swe_midcentury_vars %>% purrr::map_chr(function(var){ thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html", out_file = stringr::str_c(data_out,"/",var,"_2040-2069.nc"), bbox = mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_transform(4326) %>% # magrittr::add(c(360,0)) %>% sf::st_bbox(), vars = var, ncss_args = list(time_start = "2040-01-01", time_end = "2069-12-31"), overwrite = FALSE) }) %>% purrr::compact() %>% magrittr::set_names(bcsd_swe_midcentury_vars) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(function(x){ x %>% raster::subset( x %>% names() %>% gsub("X","",.) %>% as.Date(format = "%Y.%m.%d") %>% lubridate::month() %in% c(4) %>% which() ) %>% mm_to_in() %>% raster::as.list() }) %>% # The next eight lines first calculates the quantiles of the 31 models for each year # then the median for each percentile across the mid-century purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(raster::calc, quantile, na.rm = T) %>% purrr::map(raster::as.list) %>% purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>% purrr::map(raster::calc, median, na.rm = TRUE) %>% raster::stack(quick = TRUE) bcsd_swe_1971_2000 <- bcsd_swe_midcentury_vars %>% purrr::map_chr(function(var){ thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html", out_file = stringr::str_c(data_out,"/",var,"_1971-2000.nc"), bbox = mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_transform(4326) %>% # magrittr::add(c(360,0)) %>% sf::st_bbox(), vars = var, ncss_args = list(time_start = "1971-01-01", time_end = "2000-12-31"), overwrite = FALSE) }) %>% purrr::compact() %>% magrittr::set_names(bcsd_swe_midcentury_vars) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(function(x){ x %>% raster::subset( x %>% names() %>% gsub("X","",.) %>% as.Date(format = "%Y.%m.%d") %>% lubridate::month() %in% c(4) %>% which() ) %>% mm_to_in() %>% raster::as.list() }) %>% # The next eight lines first calculates the quantiles of the 31 models for each year # then the median for each percentile across the mid-century purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(raster::calc, quantile, na.rm = T) %>% purrr::map(raster::as.list) %>% purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>% purrr::map(raster::calc, median, na.rm = TRUE) %>% raster::stack(quick = TRUE) bcsd_swe_1981_2010 <- bcsd_swe_midcentury_vars %>% purrr::map_chr(function(var){ thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html", out_file = stringr::str_c(data_out,"/",var,"_1981-2010.nc"), bbox = mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_transform(4326) %>% # magrittr::add(c(360,0)) %>% sf::st_bbox(), vars = var, ncss_args = list(time_start = "1981-01-01", time_end = "2010-12-31"), overwrite = FALSE) }) %>% purrr::compact() %>% magrittr::set_names(bcsd_swe_midcentury_vars) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(function(x){ x %>% raster::subset( x %>% names() %>% gsub("X","",.) %>% as.Date(format = "%Y.%m.%d") %>% lubridate::month() %in% c(4) %>% which() ) %>% mm_to_in() %>% raster::as.list() }) %>% # The next eight lines first calculates the quantiles of the 31 models for each year # then the median for each percentile across the mid-century purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(raster::calc, quantile, na.rm = T) %>% purrr::map(raster::as.list) %>% purrr::transpose() %>% purrr::map(raster::stack, quick = TRUE) %>% magrittr::set_names(paste0(seq(0,1,0.25),"%")) %>% purrr::map(raster::calc, median, na.rm = TRUE) %>% raster::stack(quick = TRUE)
Here is a panel of maps of the 1981--2010 median April 1 SWE (in.) compared to the mid-century median April 1 SWE (in.), with HUC r huc watersheds overlaid.
basins <- mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_transform(4326) # pal = mapview::mapviewPalette("mapviewTopoColors") pal <- colorRampPalette(rev(RColorBrewer::brewer.pal(9,"Blues")), alpha = TRUE) m1 <- mapview::mapview(bcsd_swe_1981_2010[["X0.5."]], col.regions = pal(100), at = seq(0, 60, 5), legend = TRUE, alpha = 0.5, na.color = NA, layer.name = "1981-2010 SWE (in.)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "white") m2 <- mapview::mapview(bcsd_swe_midcentury[["X0.5."]], col.regions = pal(100), at = seq(0, 60, 5), legend = TRUE, alpha = 0.5, na.color = NA, layer.name = "Mid-century SWE (in.)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "white") mapview::sync(list(m1, m2), ncol = 1)
Downloading the 1971--2000 and 1981--2010 normal instrumental SNOTEL data {-}
Our instrumental data comes from the NRCS SNOTEL network, available through the NRCS National Water and Climate Center. The NRCS publishes two datasets: an inventory of SNOTEL sites, including their associated HUC basins, and the data recorded from the sites. It is important to remember that one SNOTEL site can be associated with multiple basins; the SNOTEL sites are often near basin boundaries, and when aggregating SNOTEL data to great basin reports, the NRCS uses both SNOTEL stations within the basin, and ones that are nearby. We also download both 1971--2000 and 1981--2010 normal instrumental SNOTEL data. The MCO has been using 1981--2010 as its current standard for defining normal climatology, but the water chapter of the MCA used 1971--2000. This section of code:
- Downloads the SNOTEL inventory for SNOTEL sites within Montana HUC
r hucwatersheds. - Downloads the normal data for those SNOTEL sites.
- Joins those two datasets, resulting in spatially-referenced SNOTEL normals data.
date <- "2018-04-01" %>% lubridate::as_date() snotel_inventory <- mcor::mco_get_snotel_inventory() %>% dplyr::left_join(mcor::mco_get_snotel_huc() %>% dplyr::mutate(`WBD code` = stringr::str_sub(`WBD code`,1,huc)) %>% dplyr::distinct()) %>% dplyr::filter(`WBD code` %in% (mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` %in% c(huc)) %$% `WBD code`)) %>% sf::st_as_sf() %>% sf::st_transform(mcor::mt_state_plane) %>% dplyr::select(-`Station Id`:-`End Date`) snotel_data <- mcor::mco_get_snotel_data(stations = snotel_inventory$Station %>% unique(), variables = c('WTEQ::normal_1981', 'WTEQ::normal_1971'), start_date = date, end_date = date) snotel_data <- snotel_inventory %>% dplyr::left_join(snotel_data, by = c("Station")) %>% dplyr::arrange(Station) %>% stats::na.omit() %>% dplyr::select(Station, `WBD code`, `Normal Snow Water Equivalent (1981-2010) (in) Start of Day Values`, `Normal Snow Water Equivalent (1971-2000) (in) Start of Day Values`) %>% dplyr::rename(`SNOTEL 1981-2010 normal (in.)` = `Normal Snow Water Equivalent (1981-2010) (in) Start of Day Values`, `SNOTEL 1971-2000 normal (in.)` = `Normal Snow Water Equivalent (1971-2000) (in) Start of Day Values`)
Here is a panel of maps of the 1981--2010 and 1971--2000 April 1 SWE values.
snotel_data_slim <- snotel_data %>% dplyr::select(-`WBD code`) %>% dplyr::distinct() %>% dplyr::rename(`Normal_1981` = "SNOTEL 1981-2010 normal (in.)", `Normal_1971` = "SNOTEL 1971-2000 normal (in.)") m1 <- mapview::mapview(snotel_data_slim, zcol = "Normal_1981", col.regions = pal(100), at = seq(0, 60, 5), legend = TRUE, alpha = 0.5, na.color = NA, layer.name = "1981-2010 Norm. (in.)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "black") m2 <- mapview::mapview(snotel_data_slim, zcol = "Normal_1971", col.regions = pal(100), at = seq(0, 60, 5), legend = TRUE, alpha = 0.5, na.color = NA, layer.name = "1971-2000 Norm. (in.)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "black") mapview::sync(list(m1, m2), ncol = 1)
Extract the BCSD VIC data from the SNOTEL locations {-}
This section of code extracts the BCSD mid-century data for the SNOTEL locations, and attaches those data to the SNOTEL data downloaded above. Notice that the BCSD values for the extracted SNOTEL locations are very different from their measured values.
snotel_bcsd_midcentury <- raster::extract(bcsd_swe_midcentury, snotel_data %>% sf::st_cast() %>% sf::st_transform(4326), df = TRUE) %>% dplyr::select(-ID) %>% magrittr::set_names(c("BCSD Midcentury min.", "BCSD Midcentury 25%", "BCSD 2049-2060 median (in.)", "BCSD Midcentury 75%", "BCSD Midcentury max.")) %>% dplyr::select(`BCSD 2049-2060 median (in.)`) snotel_bcsd_normal_1971_2000 <- raster::extract(bcsd_swe_1971_2000, snotel_data %>% sf::st_cast() %>% sf::st_transform(4326), df = TRUE) %>% dplyr::select(-ID) %>% magrittr::set_names(c("BCSD 1971-2000 min.", "BCSD 1971-2000 25%", "BCSD 1971-2000 normal (in.)", "BCSD 1971-2000 75%", "BCSD 1971-2000 max.")) %>% dplyr::select(`BCSD 1971-2000 normal (in.)`) snotel_bcsd_normal_1981_2010 <- raster::extract(bcsd_swe_1981_2010, snotel_data %>% sf::st_cast() %>% sf::st_transform(4326), df = TRUE) %>% dplyr::select(-ID) %>% magrittr::set_names(c("BCSD 1981-2010 min.", "BCSD 1981-2010 25%", "BCSD 1981-2010 normal (in.)", "BCSD 1981-2010 75%", "BCSD 1981-2010 max.")) %>% dplyr::select(`BCSD 1981-2010 normal (in.)`) snotel_data %<>% dplyr::bind_cols(snotel_bcsd_midcentury) %>% dplyr::bind_cols(snotel_bcsd_normal_1971_2000) %>% dplyr::bind_cols(snotel_bcsd_normal_1981_2010) %>% dplyr::select(Station, `WBD code`, `SNOTEL 1981-2010 normal (in.)`, `BCSD 1981-2010 normal (in.)`, `SNOTEL 1971-2000 normal (in.)`, `BCSD 1971-2000 normal (in.)`, `BCSD 2049-2060 median (in.)` ) snotel_data %>% sf::st_set_geometry(NULL) %>% dplyr::select(-`WBD code`) %>% dplyr::distinct() %>% knitr::kable(digits = 2) %>% kable_styling() %>% scroll_box(width = "100%", height = "500px")
Mapping percent of normal April 1 SWE {-}
Here, we calculate the % of normal SWE, and make some maps relating future projections to normal conditions.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9])) snotel_data_1981 <- snotel_data %>% dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1981-2010 normal (in.)`) %>% dplyr::distinct() snotel_data_1971 <- snotel_data %>% dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1971-2000 normal (in.)`) %>% dplyr::distinct() m1 <- mapview::mapview(snotel_data_1981, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1981-2010)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "black") m2 <- mapview::mapview(snotel_data_1971, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1971-2000)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "black") mapview::sync(list(m1, m2), ncol = 1)
These maps make the same comparison among the gridded modeled data.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9])) # bcsd_swe_midcentury_nozero <- bcsd_swe_midcentury # bcsd_swe_1981_2010_nozero <- bcsd_swe_1981_2010 # bcsd_swe_1971_2000_nozero <- bcsd_swe_1971_2000 # # bcsd_swe_midcentury_nozero[bcsd_swe_midcentury == 0] <- 0.0001 # bcsd_swe_1981_2010_nozero[bcsd_swe_1981_2010 == 0] <- 0.0001 # bcsd_swe_1971_2000_nozero[bcsd_swe_1971_2000 == 0] <- 0.0001 bcsd_swe_percent_normal_1981 <- 100 * bcsd_swe_midcentury / bcsd_swe_1981_2010 bcsd_swe_percent_normal_1971 <- 100 * bcsd_swe_midcentury / bcsd_swe_1971_2000 # bcsd_swe_percent_normal_1981[bcsd_swe_percent_normal_1981 > 200] <- NA # bcsd_swe_percent_normal_1971[bcsd_swe_percent_normal_1971 > 200] <- NA m1 <- mapview::mapview(bcsd_swe_percent_normal_1981[["X0.5."]], col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1981-2010)") + mapview::mapview(basins, alpha.regions = 0, label = basins$Watershed, color = "black") m2 <- mapview::mapview(bcsd_swe_percent_normal_1971[["X0.5."]] %>% round(), col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1971-2000)") + mapview::mapview(basins, alpha.regions = 0, label = FALSE, color = "black") mapview::sync(list(m1, m2), ncol = 1)
Finally, we to emulate the NRCS SWE maps. In the first set of maps, we aggregate modeled values over each HUC basin. This roughly matches what was done during the Montana Climate Assessment. In the second set, we aggregate stations within basins by taking their mean, following the procedures used by the NRCS to make their standard SWE maps.
Midcentury % of normal April 1st SWE, modeled to modeled, gridded values {-}
These maps compare modeled normal April 1 SWE values to modeled midcentury values. Gridded reconstructions are averaged over entire basins.
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9])) watersheds_cropped <- mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc)# %>% # sf::st_intersection(mt_state) %>% # dplyr::filter(sf::st_area(.) > units::set_units(1000000000,m^2)) %>% # dplyr::rename(`WBD code` = `WBD.code`, # `Hydrologic Unit` = Hydrologic.Unit) bcsd_swe_percent_normal_1981_basins <- watersheds_cropped %>% dplyr::bind_cols(bcsd_swe_percent_normal_1981 %>% raster::extract(watersheds_cropped %>% sf::st_transform(4326), fun = mean, na.rm=TRUE, df = TRUE) %>% dplyr::select(-ID) %>% dplyr::rename(`Midcentury min.` = `X0.`, `Midcentury 25%` = `X0.25.`, `Midcentury median` = `X0.5.`, `Midcentury 75%` = `X0.75.`, `Midcentury max.` = `X1.`)) %>% mutate_all(funs(replace(., is.na(.), NA))) %>% dplyr::rename("Midcentury_median" = "Midcentury median") %>% dplyr::select(Watershed, Midcentury_median) %>% na.omit() bcsd_swe_percent_normal_1971_basins <- watersheds_cropped %>% dplyr::bind_cols(bcsd_swe_percent_normal_1971 %>% raster::extract(watersheds_cropped %>% sf::st_transform(4326), fun = mean, na.rm=TRUE, df = TRUE) %>% dplyr::select(-ID) %>% dplyr::rename(`Midcentury min.` = `X0.`, `Midcentury 25%` = `X0.25.`, `Midcentury median` = `X0.5.`, `Midcentury 75%` = `X0.75.`, `Midcentury max.` = `X1.`)) %>% mutate_all(funs(replace(., is.na(.), NA))) %>% dplyr::rename("Midcentury_median" = "Midcentury median") %>% dplyr::select(Watershed, Midcentury_median) %>% na.omit() m1 <- mapview::mapview(bcsd_swe_percent_normal_1981_basins, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1981-2010), gridded", label = paste0(bcsd_swe_percent_normal_1981_basins$Watershed,": ",bcsd_swe_percent_normal_1981_basins$Midcentury_median %>% round(),"%")) m2 <- mapview::mapview(bcsd_swe_percent_normal_1971_basins, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1971-2000), gridded", label = paste0(bcsd_swe_percent_normal_1971_basins$Watershed,": ",bcsd_swe_percent_normal_1971_basins$Midcentury_median %>% round(),"%")) mapview::sync(list(m1, m2), ncol = 1)
We can also look at the three basins in figure 3-10 of the Montana Climate Assessment, in an effort to reproduce that figure. Here, we calculate the basin-wide boxplots for three basins roughly conforming to the basins in the MCA, using the 1971--2000 normals; these don't conform precisely with figure 3-10, as the pour points for defining the basins are underspecified in the documentation for for Montana Climate Assessment. Here, we define the Clark Fork above St. Regis as including the Middle Clark Fork, Bitterroot, Flint-Rock, Blackfoot, and Upper Clark Fork HUC 8 basins; the Missouri above Toston as including the Gallatin, Madison, Jefferson, Ruby, Beaverhead, and Big Hole basins; and the Yellowstone above Billings as including the Upper Yellowstone, Stillwater, Clarks Fork Yellowstone, Shields, and Yellowstone Headwaters basins. This analysis includes points for each of 31 climate models available in the BCSD VIC dataset.
bcsd_swe_1971_2000_all_models <- bcsd_swe_midcentury_vars %>% purrr::map_chr(function(var){ thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html", out_file = stringr::str_c(data_out,"/",var,"_1971-2000.nc"), bbox = mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_transform(4326) %>% # magrittr::add(c(360,0)) %>% sf::st_bbox(), vars = var, ncss_args = list(time_start = "1971-01-01", time_end = "2000-12-31"), overwrite = FALSE) }) %>% purrr::compact() %>% magrittr::set_names(bcsd_swe_midcentury_vars) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(function(x){ x %>% raster::subset( x %>% names() %>% gsub("X","",.) %>% as.Date(format = "%Y.%m.%d") %>% lubridate::month() %in% c(4) %>% which() ) %>% mm_to_in() %>% raster::as.list() }) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(mean) %>% raster::stack(quick = TRUE) bcsd_swe_midcentury_all_models <- bcsd_swe_midcentury_vars %>% purrr::map_chr(function(var){ thredds::tds_ncss_download(ncss_url = "https://cida.usgs.gov/thredds/ncss/BCSD_mon_VIC/dataset.html", out_file = stringr::str_c(data_out,"/",var,"_2040-2069.nc"), bbox = mcor::mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_transform(4326) %>% # magrittr::add(c(360,0)) %>% sf::st_bbox(), vars = var, ncss_args = list(time_start = "2040-01-01", time_end = "2069-12-31"), overwrite = FALSE) }) %>% purrr::compact() %>% magrittr::set_names(bcsd_swe_midcentury_vars) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(function(x){ x %>% raster::subset( x %>% names() %>% gsub("X","",.) %>% as.Date(format = "%Y.%m.%d") %>% lubridate::month() %in% c(4) %>% which() ) %>% mm_to_in() %>% raster::as.list() }) %>% purrr::map(raster::stack, quick = TRUE) %>% purrr::map(mean) %>% raster::stack(quick = TRUE) watersheds_cropped <- mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) %>% sf::st_transform(4326) %>% dplyr::mutate(.,ID = 1:nrow(.)) (bcsd_swe_midcentury_all_models / bcsd_swe_1971_2000_all_models) %>% raster::extract(watersheds_cropped, df = TRUE) %>% tidyr::gather(Model, Value, -ID) %>% # dplyr::rename(`Midcentury median` = `X0.5.`) %>% dplyr::select(ID, Value) %>% dplyr::left_join(x = watersheds_cropped, y = .) %>% dplyr::filter(!is.nan(Value), !is.na(Value), !is.infinite(Value)) %>% dplyr::mutate(Watershed = factor(Watershed)) %>% dplyr::mutate(Watershed = forcats::fct_collapse(Watershed, `Clark Fork above St. Regis` = c("Middle Clark Fork", "Bitterroot", "Flint-Rock", "Blackfoot", "Upper Clark Fork"), `Missouri above Toston` = c("Gallatin", "Madison", "Jefferson", "Ruby", "Beaverhead", "Big Hole"), `Yellowstone above Billings` = c("Upper Yellowstone", "Stillwater", "Clarks Fork Yellowstone", "Shields", "Yellowstone Headwaters"))) %>% dplyr::filter(Watershed %in% c("Clark Fork above St. Regis", "Missouri above Toston", "Yellowstone above Billings"), Value >= 0, Value <= 2) %>% dplyr::mutate(Watershed = forcats::fct_drop(Watershed), Watershed = forcats::fct_relevel(Watershed, "Clark Fork above St. Regis", "Missouri above Toston", "Yellowstone above Billings")) %>% # sf::st_set_geometry(NULL) %>% # dplyr::select(Watershed, `Midcentury median`) %>% ggplot2::ggplot(ggplot2::aes(x = Watershed, y = 100 * (Value - 1))) + ggplot2::geom_boxplot(outlier.shape = NA) + ylab("% Change in April 1 SWE")
Midcentury % of normal April 1st SWE, station to modeled, point values {-}
These maps compare empirical April 1 SWE values from SNOTEL stations to the modeled midcentury SWE values extracted from those same station locations (i.e., the data from the station table above).
pal <- colorRampPalette(c(RColorBrewer::brewer.pal(9,"RdBu")[1:3],RColorBrewer::brewer.pal(9,"RdBu")[7:9])) watersheds_cropped <- mt_watersheds_simple %>% dplyr::filter(`Hydrologic Unit` == huc) snotel_data_1981 <- snotel_data %>% dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1981-2010 normal (in.)`) %>% group_by(`WBD code`) %>% dplyr::summarise(`Stations Count` = n(), `Midcentury_median` = mean(Midcentury_median) %>% round()) %>% dplyr::filter(`Stations Count` >= 3) %>% sf::st_set_geometry(NULL) %>% dplyr::left_join(watersheds_cropped) %>% sf::st_sf() %>% na.omit() %>% dplyr::select(-`Hydrologic Unit`) %>% dplyr::select(`WBD code`, Watershed, dplyr::everything()) snotel_data_1971 <- snotel_data %>% dplyr::mutate(Midcentury_median = 100 * `BCSD 2049-2060 median (in.)` / `SNOTEL 1971-2000 normal (in.)`) %>% group_by(`WBD code`) %>% dplyr::summarise(`Stations Count` = n(), `Midcentury_median` = mean(Midcentury_median) %>% round()) %>% dplyr::filter(`Stations Count` >= 3) %>% sf::st_set_geometry(NULL) %>% dplyr::left_join(watersheds_cropped) %>% sf::st_sf() %>% na.omit() %>% dplyr::select(-`Hydrologic Unit`) %>% dplyr::select(`WBD code`, Watershed, dplyr::everything()) m1 <- mapview::mapview(snotel_data_1981, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1981-2010)", label = paste0(snotel_data_1981$Watershed,": ",snotel_data_1981$Midcentury_median,"%")) m2 <- mapview::mapview(snotel_data_1971, zcol = "Midcentury_median", col.regions = pal(100), at = seq(0, 200, 20), legend = TRUE, na.color = NA, layer.name = "% Norm. (1971-2000)", label = paste0(snotel_data_1971$Watershed,": ",snotel_data_1971$Midcentury_median,"%")) mapview::sync(list(m1, m2), ncol = 1)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.