In bcgov/forestDroughtTool: Predict Forest Drought Hazard In British Columbia

knitr::opts_chunk$set(echo = TRUE)

Introduction

The purpose of this vignette is to document our process for selecting climate data as input to the asmrCalc() function in the forestDroughtTool R package.

# Load libraries
library(bcmaps)
library(weathercan)
library(dplyr)
library(magrittr)
library(ggplot2)
library(cowplot)
library(forestDroughtTool)

Climate station coverage for BGC units in BC

Environment Canada (EC) daily climate stations seem to be the best so far because they have relatively long records and they record winter precipitation.

We can use a few R packages to identify EC stations with daily data and their corresponding BGC unit.

# Step 1.1 - Select stations in BC with daily data and assign to 'stn'
stn<-
  stations %>% 
  filter(prov=="BC" & interval=="day") %>%  

# Step 1.2 - Convert to a spatial file and merge with BGC units (this will take awhile!)
  st_as_sf(coords=c("lon","lat")) %>% # convert to spatial file
  st_set_crs(4326) %>% # set to WGS1984 datum
  transform_bc_albers() %>% # set to BC albers projection
  st_join(bec()[,"MAP_LABEL"]) %>%  # merge with BEC
  st_drop_geometry() # drop geometry

# Step 1.3 Summarize station coverage
  stn %>% 
    group_by(MAP_LABEL) %>%
    summarise(Num.Stn=n())

Subset stations for BGC units of interest

We can now download data for BGC unit(s) of interest.

# List the BGC units of interest
bgc<-c("SBSdk","SBSwk1")

# Filter stations for those units
stnBGC<-
  stn %>% 
  filter(MAP_LABEL%in%bgc) %>% # filter for BGC units of interest
  mutate(length=end-start+1) %>% # create a column to show record length
  arrange(desc(length)) %>% 
  select(station_name,station_id,MAP_LABEL,start,end,length)

stnBGC

Note there aren't many records for SBSwk1. Might want to drop the variant from the BGC unit, and see if there any any stations in that subzone but in other variants.

Downloading of daily climate data

For this vignette, let's download the Smithers A climate station data. From the above table, we see that the station ID is 487. We'll need this number for the weather download. We'll also limit the download to 1955 to 1995 (somewhat around the 1961-1990 climate normal period.)

# Download climate station data (this can take a long time!)
sm.dat<-
  weather_dl(station_ids=487,interval="day",start="1955-01-01",end="1995-12-31") %>%

  # format climate data (see below for more information)
  select(stn=station_name,date,tmn=min_temp,tmx=max_temp,ppt=total_precip,year,month,day)

You can download data for multiple stations at once.

Formatting of daily climate data

Daily climate data have to be formatted correctly for the asmrCalc() function to work properly. There have to seven columns in the data, and they must be named exactly as follows: - date - tmn - tmx - ppt - year - month - day

As an example, refer to the Prince George daily climate data in the forestDroughtTool package:

PrinceGeorge

Plot data to confirm data coverage

# Create plots 
  p1<-ggplot(sm.dat,aes(x=date,y=ppt))+geom_line()
  p2<-ggplot(sm.dat,aes(x=date,y=tmx))+geom_line()
  p3<-ggplot(sm.dat,aes(x=date,y=tmn))+geom_line()
  plot_grid(p1,p2,p3)

Summarise NA values in data

We need to process missing values in the daily climate data. For this project, we: - omitted any years with >10 consecutive missing data in any climate variable;and - imputed missing data using adjacent values (closest data before and after missing value)

sm.dat_filled<-cleanECData(sm.dat)

If you are cleaning data for multiple station files downloaded into one data frame, use the following code:

# Not run
# climData is multi-station data downloaded using the weathercan package

# X<-select(climData,stn=station_name,date,tmn=min_temp,tmx=max_temp,ppt=total_precip,year,month,day)
# x1<-by(INDICES=X$stn,function(x) cleanECData(x),data=X)
# x2<-dplyr::bind_rows(x1,.id="id")
# x2$stn_name=names(x1)[as.numeric(x2$id)]

If you want to know what years were omitted from analysis, and what daily values were imputed, use the following code:

bcgov/forestDroughtTool documentation built on March 20, 2021, 4:15 p.m.