class: title-slide, left, middle <h1> Obtaining, cleaning and visualizing <br> hydrological data with <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 581 512"><path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"/></svg> </h1> <br> <h3> Alex Hurley <br> <em>University of Birmingham</em> </h3> <br> <br> .small[<svg style="height:0.8em;top:.04em;position:relative;fill:lightgrey;" viewBox="0 0 576 512"><path d="M488 312.7V456c0 13.3-10.7 24-24 24H348c-6.6 0-12-5.4-12-12V356c0-6.6-5.4-12-12-12h-72c-6.6 0-12 5.4-12 12v112c0 6.6-5.4 12-12 12H112c-13.3 0-24-10.7-24-24V312.7c0-3.6 1.6-7 4.4-9.3l188-154.8c4.4-3.6 10.8-3.6 15.3 0l188 154.8c2.7 2.3 4.3 5.7 4.3 9.3zm83.6-60.9L488 182.9V44.4c0-6.6-5.4-12-12-12h-56c-6.6 0-12 5.4-12 12V117l-89.5-73.7c-17.7-14.6-43.3-14.6-61 0L4.4 251.8c-5.1 4.2-5.8 11.8-1.6 16.9l25.5 31c4.2 5.1 11.8 5.8 16.9 1.6l235.2-193.7c4.4-3.6 10.8-3.6 15.3 0l235.2 193.7c5.1 4.2 12.7 3.5 16.9-1.6l25.5-31c4.2-5.2 3.4-12.7-1.7-16.9z"/></svg> [aglhurley.rbind.io](https://aglhurley.rbind.io) <svg style="height:0.8em;top:.04em;position:relative;fill:lightgrey;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg> [aglhurley](https://twitter.com/aglhurley)] .title-logo-box[![](./static/img/rhydro_logo_alt.png)] --- # Goals -- .middle[ - <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 512 512"><path d="M239.1 6.3l-208 78c-18.7 7-31.1 25-31.1 45v225.1c0 18.2 10.3 34.8 26.5 42.9l208 104c13.5 6.8 29.4 6.8 42.9 0l208-104c16.3-8.1 26.5-24.8 26.5-42.9V129.3c0-20-12.4-37.9-31.1-44.9l-208-78C262 2.2 250 2.2 239.1 6.3zM256 68.4l192 72v1.1l-192 78-192-78v-1.1l192-72zm32 356V275.5l160-65v133.9l-160 80z"/></svg> Introduce **useful packages** - <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 640 512"><path d="M434.7 64h-85.9c-8 0-15.7 3-21.6 8.4l-98.3 90c-.1.1-.2.3-.3.4-16.6 15.6-16.3 40.5-2.1 56 12.7 13.9 39.4 17.6 56.1 2.7.1-.1.3-.1.4-.2l79.9-73.2c6.5-5.9 16.7-5.5 22.6 1 6 6.5 5.5 16.6-1 22.6l-26.1 23.9L504 313.8c2.9 2.4 5.5 5 7.9 7.7V128l-54.6-54.6c-5.9-6-14.1-9.4-22.6-9.4zM544 128.2v223.9c0 17.7 14.3 32 32 32h64V128.2h-96zm48 223.9c-8.8 0-16-7.2-16-16s7.2-16 16-16 16 7.2 16 16-7.2 16-16 16zM0 384h64c17.7 0 32-14.3 32-32V128.2H0V384zm48-63.9c8.8 0 16 7.2 16 16s-7.2 16-16 16-16-7.2-16-16c0-8.9 7.2-16 16-16zm435.9 18.6L334.6 217.5l-30 27.5c-29.7 27.1-75.2 24.5-101.7-4.4-26.9-29.4-24.8-74.9 4.4-101.7L289.1 64h-83.8c-8.5 0-16.6 3.4-22.6 9.4L128 128v223.9h18.3l90.5 81.9c27.4 22.3 67.7 18.1 90-9.3l.2-.2 17.9 15.5c15.9 13 39.4 10.5 52.3-5.4l31.4-38.6 5.4 4.4c13.7 11.1 33.9 9.1 45-4.7l9.5-11.7c11.2-13.8 9.1-33.9-4.6-45.1z"/></svg> Highlight synergies - <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 640 512"><path d="M629.657 343.598L528.971 444.284c-9.373 9.372-24.568 9.372-33.941 0L394.343 343.598c-9.373-9.373-9.373-24.569 0-33.941l10.823-10.823c9.562-9.562 25.133-9.34 34.419.492L480 342.118V160H292.451a24.005 24.005 0 0 1-16.971-7.029l-16-16C244.361 121.851 255.069 96 276.451 96H520c13.255 0 24 10.745 24 24v222.118l40.416-42.792c9.285-9.831 24.856-10.054 34.419-.492l10.823 10.823c9.372 9.372 9.372 24.569-.001 33.941zm-265.138 15.431A23.999 23.999 0 0 0 347.548 352H160V169.881l40.416 42.792c9.286 9.831 24.856 10.054 34.419.491l10.822-10.822c9.373-9.373 9.373-24.569 0-33.941L144.971 67.716c-9.373-9.373-24.569-9.373-33.941 0L10.343 168.402c-9.373 9.373-9.373 24.569 0 33.941l10.822 10.822c9.562 9.562 25.133 9.34 34.419-.491L96 169.881V392c0 13.255 10.745 24 24 24h243.549c21.382 0 32.09-25.851 16.971-40.971l-16.001-16z"/></svg> Showcase capabilities (**processing**, **stats**) - <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 512 512"><path d="M500 384c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12H12c-6.6 0-12-5.4-12-12V76c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v308h436zM372.7 159.5L288 216l-85.3-113.7c-5.1-6.8-15.5-6.3-19.9 1L96 248v104h384l-89.9-187.8c-3.2-6.5-11.4-8.7-17.4-4.7z"/></svg> Visualize (**static** and **interacive** graphs)] --- # Approach - Download **hydrometric** and **ancillary** (gridded) data for stations defined with `AOI` in Glacier National Park (BC, Canada) - Calculate *P-pET* (`raster`), *Runoff*, and flow statistics (`fasstr`) - Visualize results (P-pET vs. R, interactive flow stats) .center[<img src="./03_get-clean-viz_files/ext-figure/flow-chart-pckgs.png" alt="drawing" width="450"/>] --- layout: false class: inverse, center, middle # Get Started --- # Set-up: packages and utils ```r install.packages("devtools") # hydro and met devtools::install_github("bcgov/fasstr") install.packages("daymetr") install.packages("tidyhydat") # spatial devtools::install_github("mikejohnson51/AOI") install.packages("rgeos") install.packages("raster") install.packages("leaflet") # general purpose and viz install.packages("dplyr") install.packages("purrr") install.packages("ggplot2") install.packages("plotly") install.packages("DT") ``` -- - Relevant packages loaded when necessary - Typically use `package::function()` for clarity --- class: inverse, center, middle # A: Define area of interest --- # Define area - use `AOI` to define bounding box around Glacier National Park (100 by 100 km) ```r # Specify Regions and have a peak rockies <- AOI::getAOI(clip = list("Glacier National Park Canada", 100, 100), km = TRUE) class(rockies) ``` ``` ## [1] "SpatialPolygons" ## attr(,"package") ## [1] "sp" ``` - `rockies@bbox` useful for other functions! --- layout: true # Define area --- ```r library(dplyr) library(leaflet) rockies %>% AOI::check() ``` ---
--- layout:false class: inverse, center, middle # B: Climate Data --- # Download climate data - `daymetr`: gridded climate data for North America as NetCDF - downloads files to defined path - provided in LCC projection ```r # for rockies path_rockies <- "./03_get-clean-viz_files/data/rockies" # download data from daymet service *params <- c("dayl", "tmin", "tmax", "prcp") params %>% purrr::walk( * ~daymetr::download_daymet_ncss(location = c(rockies@bbox[[4]], rockies@bbox[[1]], rockies@bbox[[2]], rockies@bbox[[3]]), # top left to bottom right start = 2010, end = 2011, param = .x, frequency = "daily", path = path_rockies)) ``` --- # Load climate rasters - Define projections ```r # Projected CS proj4.Lambert <- "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +a=6378137 +b=6356752.314706705 +units=km +no_defs" # Geographic CS proj4.WGS <- "+init=epsg:4326" ``` --- # Load climate rasters - read individual NetCDFs, stack as rasters - reproject ```r path_rockies <- "./03_get-clean-viz_files/data/rockies" # load data, reproject *params <- c("dayl", "tmin", "tmax", "prcp") rockies_stacks <- params %>% purrr::map( function(x){ list.files(path_rockies, pattern = x, full.names = TRUE) %>% * raster::stack() %>% * raster::`projection<-`(., proj4.Lambert) %>% * raster::projectRaster(crs = proj4.WGS) } ) %>% setNames(params) ``` --- # Load climate rasters ```r # Jan 1, 2010 + 399 days raster::plot(rockies_stacks$prcp[[400]]) ``` <img src="03_get-clean-viz_files/figure-html/raster-plot-init-1.png" width="360" style="display: block; margin: auto;" /> --- layout: true # Calculate P - pET --- - based on `ET.Hamon()` from `Ecohydrology` <svg style="height:0.8em;top:.04em;position:relative;fill:#FFDD00;" viewBox="0 0 512 512"><path d="M239.1 6.3l-208 78c-18.7 7-31.1 25-31.1 45v225.1c0 18.2 10.3 34.8 26.5 42.9l208 104c13.5 6.8 29.4 6.8 42.9 0l208-104c16.3-8.1 26.5-24.8 26.5-42.9V129.3c0-20-12.4-37.9-31.1-44.9l-208-78C262 2.2 250 2.2 239.1 6.3zM256 68.4l192 72v1.1l-192 78-192-78v-1.1l192-72zm32 356V275.5l160-65v133.9l-160 80z"/></svg> - applied to basic raster math ```r # custom function to calculate Hamon's PET et.ham <- function(tmin,tmax,dayl){ # modified from Evapotranspiration package Ta <- (tmax + tmin)/2 vs_Tmax <- 0.6108 * exp(17.27 * tmax/(tmax + 237.3)) vs_Tmin <- 0.6108 * exp(17.27 * tmin/(tmin + 237.3)) vas <- (vs_Tmax + vs_Tmin)/2 ET_Hamon.Daily <- 0.55 * 25.4 * (dayl/12)^2 * (216.7 * vas * 10/(Ta + 273.3))/100 return(ET_Hamon.Daily) } # calculate PET over all days (2 years total) *et_rockies <- et.ham(tmin = rockies_stacks$tmin, tmax = rockies_stacks$tmax, dayl = rockies_stacks$dayl / 3600) # remove calc. artefacts (set to NA from previous raster) et_rockies <- raster::mask(et_rockies, rockies_stacks$dayl[[1]]) ``` --- - aggregate to monthly values ```r # set up indices for aggregating over months year_mon <- seq(as.Date("2010-01-01"), as.Date("2011-12-31"), by = "1 day") %>% format("%Y-%m") months <- as.numeric(as.factor(year_mon)) # monthly totals of et *et_monthly <- raster::stackApply(et_rockies, months, fun = sum) # monthly totals of p *p_monthly <- raster::stackApply(rockies_stacks$prcp, months, fun = sum) # rough balance: p - pet *pet_monthly <- p_monthly - et_monthly # mask NA values for clean plotting pet_monthly <- raster::mask(pet_monthly, rockies_stacks$dayl[[1]]) ``` --- layout: true # Calculate P - pET --- ```r # interactive for June leaflet() %>% addPolygons(data = rockies, fillColor = "transparent", color = col, weight = 4) %>% addProviderTiles("Esri.NatGeoWorldMap", group = "Terrain") %>% addRasterImage(x = pet_monthly[[6]], opacity = .8) ``` ---
--- layout: false class: inverse, middle, center # B: Flow data --- layout: true # Identify viable stations --- - use bounding box from `AOI` to filter for stations in `tidyhydat` data base (lat, lon). ```r # requires download *tidyhydat::hy_set_default_db(hydat_path = "D:/ext_data_R/Hydat.sqlite3") *rockies_stns <- tidyhydat::allstations %>% filter(between(LONGITUDE, rockies@bbox[1,1], rockies@bbox[1,2]), between(LATITUDE, rockies@bbox[2,1], rockies@bbox[2,2]), HYD_STATUS == "ACTIVE") ``` -- <table> <thead> <tr> <th style="text-align:left;"> STATION_NUMBER </th> <th style="text-align:left;"> STATION_NAME </th> <th style="text-align:left;"> HYD_STATUS </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> 08NA002 </td> <td style="text-align:left;"> COLUMBIA RIVER AT NICHOLSON </td> <td style="text-align:left;"> ACTIVE </td> </tr> <tr> <td style="text-align:left;"> 08NA006 </td> <td style="text-align:left;"> KICKING HORSE RIVER AT GOLDEN </td> <td style="text-align:left;"> ACTIVE </td> </tr> <tr> <td style="text-align:left;"> 08NB005 </td> <td style="text-align:left;"> COLUMBIA RIVER AT DONALD </td> <td style="text-align:left;"> ACTIVE </td> </tr> <tr> <td style="text-align:left;"> 08NB012 </td> <td style="text-align:left;"> BLAEBERRY RIVER ABOVE WILLOWBANK CREEK </td> <td style="text-align:left;"> ACTIVE </td> </tr> </tbody> </table> --- layout: true # Download hydrometric data --- - `fasstr::screen_flow_data()` for missing data checks ```r *rockies_q <- tidyhydat::hy_daily_flows(station_number = rockies_stns$STATION_NUMBER, start_date = "2000-01-01", end_date = "2015-12-31") check_me <- rockies_q %>% * fasstr::screen_flow_data() ``` --- - aggregate flow over months ```r rockies_q_month <- rockies_q %>% filter(Date >= "2010-01-01", Date <= "2011-12-31") %>% # cleaning * fasstr::fill_missing_dates() %>% # add info (table joining) * fasstr::add_basin_area() %>% # calcs mutate(r_mm_day = Value / Basin_Area_sqkm * 86400 / 1e6 * 1e3) %>% # aggregate group_by(STATION_NUMBER, year_mon = format(Date, "%Y-%m")) %>% * summarise(r_mm_month = sum(r_mm_day, na.rm = TRUE)) %>% mutate(date_time = lubridate::ymd(paste0(year_mon, "-01"), tz = "MST")) ``` --- layout: false # Extract P-pET for stations ```r # extract data from raster at station locations # Typically would use catchment shape files *pet_stns <- raster::extract(pet_monthly, * y = rockies_stns[ ,c("LONGITUDE","LATITUDE")]) %>% as.data.frame() %>% bind_cols(rockies_stns[ ,"STATION_NUMBER"]) %>% setNames(c(year_mon %>% unique(), "stn")) %>% tidyr::gather(-stn, key = "year_mon", value = "p_pet_mm") %>% mutate(date_time = lubridate::ymd(paste0(year_mon, "-01"), tz = "MST")) ``` --- layout: true # Runoff vs. P-pET --- ```r library(ggplot2) pet_r_plot <- pet_stns %>% filter(stn %in% rockies_q_month$STATION_NUMBER) %>% ggplot(aes(x = date_time, y = p_pet_mm)) + # add climate data geom_bar(stat = "identity", aes(fill = "P-pET"), alpha = 0.5) + # add R data geom_bar(inherit.aes = FALSE, data = rockies_q_month %>% rename(stn = STATION_NUMBER), aes(x = date_time, y = r_mm_month, fill = "R"), stat = "identity", alpha = 0.5) + geom_hline(yintercept = 0, linetype = 2) + # misc theme_pub(base_size = 18) + labs(x = "Date", y = "Flux (mm/month)", fill = "Measure") + coord_flip() + facet_wrap(~stn) ``` --- <img src="03_get-clean-viz_files/figure-html/pet-r-exec-1.png" width="720" style="display: block; margin: auto;" /> --- layout: false class: inverse, middle, center # C: Flow stats in <br> # interactive visualization --- # Outcome: Interactive visualization with two panels showing: .center[<img src="./03_get-clean-viz_files/ext-figure/interactive-outcome.png" alt="drawing" width="450"/>] for Beaver River station --- layout: true # Prepare flow data --- - `fasstr::calc_daily_stats()` for a range of flow statistics over chosen interval, e.g. (min, max, Q5, Q25, Q75, Q95, mean/median) ```r # choose a station beaver_river <- rockies_stns %>% filter(STATION_NUMBER == "08NB019") # calculate flow stats *beaver_q_stats <- fasstr::calc_daily_stats(rockies_q %>% * filter(STATION_NUMBER == beaver_river$STATION_NUMBER)) # add 2011 flow data for interactive viz beaver_q_stats <- beaver_q_stats %>% left_join(rockies_q %>% filter(STATION_NUMBER == beaver_river$STATION_NUMBER, Date >= "2010-01-01", Date < "2011-01-01") %>% mutate(DayofYear = lubridate::yday(Date)), by = "DayofYear") ``` ---
--- layout: false # Prepare precipitation data - Equivalent to earlier raster extraction ```r # extract P data and make df *p_beaver <- raster::extract(rockies_stacks$prcp, * y = beaver_river[ ,c("LONGITUDE","LATITUDE")]) %>% as.data.frame() %>% bind_cols(beaver_river %>% select(STATION_NUMBER, STATION_NAME)) %>% setNames(c(seq(as.Date("2010-01-01"), as.Date("2011-12-31"), by = "1 day") %>% as.character(), "stn", "name")) %>% tidyr::gather(-stn, -name, key = "date", value = "p_mm_day") %>% mutate(date_time = lubridate::ymd(date, tz = "MST")) ``` --- layout: true # Set-up interactive plots --- ## Flow stats ```r library(plotly) *q_beaver <- plot_ly(data = beaver_q_stats, x = ~DayofYear) %>% # Add max range add_ribbons(ymin = ~Minimum, ymax = ~Maximum, color = I("gray80"), name = "Max. Range (2000 - 2015)") %>% # add percentiles add_ribbons(ymin = ~P5, ymax = ~P95, color = I("steelblue1"), name = "Q5-Q95 (2000 - 2015)") %>% # add 20 add_lines(y = ~Value, color = I("darkorange"), name = "Beaver River (2010)") %>% layout(yaxis = list(title = "mean daily Q (Cm/s)")) ``` --- ## Precip ```r precip_beaver <- p_beaver %>% filter(date_time >= "2010-01-01", date_time < "2011-01-01") %>% mutate(DayofYear = lubridate::yday(date_time)) %>% # set up plot * plot_ly(data = ., type = "bar", name = "Precip (2010)") %>% * add_bars(x = ~DayofYear, * y = ~p_mm_day, color = I("steelblue1")) %>% layout(yaxis = list(title = "P (mm/day)"), width = 800) ``` --- layout: true # Generate final visualization --- ```r subplot(precip_beaver, q_beaver, nrows = 2, shareX = TRUE, titleY = TRUE) ``` ---
--- layout: false ## Summary .center[<img src="./03_get-clean-viz_files/ext-figure/flow-chart-pckgs.png" alt="drawing" width="450"/>] -- Used range of general-purpose and hydrologic packages in concert to: - obtain spatial and time series data (climate, hydrometric) - wrangle, clean and pre-process data - make a range of static and interactive visualization (maps, figures)