moose_factor_okanagan.Rmd

---
title: 'Okanagan Timber Supply Review 2020: Moose Analysis'
author: "Tyler Muhly"
date: "25/05/2020"
output:
  word_document: default
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require (sf)
require (data.table)
require (ggplot2)
require (here)
require (ggspatial)
require (dplyr)
require (tidyr)
require (scales)
```

## 1. Importance of Moose in British Columbia and the Okanagan Timber Supply Area
Recent court decisions in British Columbia (BC) have stated that statutory decision makers must use credible information to consider the effects of land management decisions, including allowable annual cut (AAC) determinations, on Indigenous peoples rights to harvest wildlife (e.g., see [William v. British Columbia 2012](https://www.bccourts.ca/jdb-txt/CA/12/02/2012BCCA0285.htm) and [West Moberly First Nations v. British Columbia 2011](https://www.bccourts.ca/jdb-txt/CA/11/02/2011BCCA0247.htm)). Accommodations to First Nations communities for impacts of forest harvest on wildlife may be necessary as part of some AAC decisions. 

In general, moose are a highly valued wildlife species in North America for a variety of cultural and economic reasons ( [Timmermann and Rogers 2005](https://pdfs.semanticscholar.org/55d2/86e83dd7cae6da270c153858b3d1a98d89de.pdf)). Prior to European colonization, moose were used extensively by Indigenous peoples for food, clothing and shelter ( [Moose In British Columbia](http://www.env.gov.bc.ca/wld/documents/moose.pdf)). In the last 100 years, moose have become an important, nutritious, staple food of many interior and coastal First Nations communities in BC ( [First Nations Health Authority fact sheet](https://www.fnha.ca/Documents/Traditional_Food_Fact_Sheets.pdf)). 

Several meetings (Dec. 15, 2016, January 30, 2017 and April 12, 2017) were held with the Secwepemc Nations, including Adam's Lake Indian Band and Splatsin First Nations to discuss key wildlife values in the Okanagan timber supply area (TSA) that should be considered as part of the timber supply review (TSR). Sustainability of moose populations was identfied as a priority value by the Secwepemc Nations. In addition, the Secwempemc Nations previously developed a collaborative [Moose and Watershed stewarship pilot porgram](https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/sustainable-forest-management-practices/moose-watershed-stewardship-pilot) with the Thompson Rivers Natural Resource District to improve moose sustainability in some portions of the Okanagan TSA. 

## 2. Key Relationships Between Moose and Forestry   
Research has shown that forestry activity influences moose density and distribution, both positively and negatively. Forestry cutblocks remove forest canopy, which generally increases the production of deciduous shrub browse on landscapes, which can positively influence moose. Shrub production in forestry cutblocks varies, but appears to peak anywhere from at 5 to 30 years after harvest, and moose appear to use these stand ages the most. Cutblocks less than 5 years old and older than 30 years old (i.e., mid-seral stands), appear to generally be of less forage value to moose and receive less use as foraging habitat. However, older stands, including mature cutblocks, can benefit moose by providing valuable cover habitat. Closed canopy conifer forest stands are important habitat for providing thermal cover in summer and to intercept snow in high snowfall years or areas. 

Forest harvest, when done in moderation and in a way that creates a diversity of forest stand ages and types, can benefit moose. However, the creation of road infrastrucutre to extract timber may negatively affect moose density and distribution overall. Forestry roads can make areas more accessible to moose hunters, increasing hunting mortality, which can limit moose population size. The overall effects of forestry on moose may be negative when roads endure on the landscape and are not actively decommissioned or recovered. 

Moose have an important and complex relationship with forestry development. Based on a review of previous research on the effects of forestry on moose, the relationship is most likely to be positive when:
  - At the scale of a moose home range (i.e., approximately 10 km^2^), forest cutblocks are interspersed with large mature or old forest stands, and cut in a progressive way over a 5 to 10 year period so there is a distribution of cutblock ages 
  - Silvicultural practices on harvested stands allow for growth of some shrubs, particularly along cutblock edges with mature stands   
  - Roads are minimized, blocked, deactivated or restored

## 3. Forestry-related Indicators of Moose Habitat and Population Condition
Based on previous research, and what can be simulated from TSR models, we used the following indicators to assess current and future conditions of moose habitat in the Okanagan TSA:
  - percentage of watershed area that is 5 to 30 years old, with an ideal percentage of 30%
  - percentage of watershed area that is conifer stands greater than 5 ha in size and 15m tall, with an ideal percentage of 40%
  - road density in a watershed area, with an ideal target of less than 1km/km^2^ 

## 4. Current State of Moose Habitat and Populations in the Timber Supply Area
Spatial data were downloaded from [DataBC](https://data.gov.bc.ca/) on April 7, 2020, and saved into a file geodatabase. Data that were downloaded included the [digital road atlas (DRA)](https://catalogue.data.gov.bc.ca/dataset/digital-road-atlas-dra-master-partially-attributed-roads), [forest tenure (FTEN) roads](https://catalogue.data.gov.bc.ca/dataset/forest-tenure-road-segment-lines), [forest vegetation resources inventory (VRI)](https://catalogue.data.gov.bc.ca/dataset/vri-2019-forest-vegetation-composite-rank-1-layer-r1-), [TSA boundaries](https://catalogue.data.gov.bc.ca/dataset/fadm-timber-supply-area-tsa), [freshwater atlas assessment watershed areas (AWAs)](https://catalogue.data.gov.bc.ca/dataset/freshwater-atlas-assessment-watersheds), [landscape unit (LU) boundaries](https://catalogue.data.gov.bc.ca/dataset/landscape-units-of-british-columbia-current) and [wildlife management unit (WMU) boundaries](https://catalogue.data.gov.bc.ca/dataset/wildlife-management-units). 

The DRA, FTEN roads, freshwater atlas and VRI data were 'clipped' to the Okanagan TSA boundary. The linear DRA and FTEN roads data were merged together and then converted to a 20 m resolution raster to remove duplicate roads in both datasets (i.e., roads less than 20 m apart). I then converted data back to linear data using the ArcScan extension, with the following settings: 
- Geometrical intersection 
- Max line width = 20 
- Noise level = 20 
- Compression tolerance = 0.025 
- Smoothing weight = 3 
- Hole size = 0 

Roads data, and data on forest stand crown closure, species, projected age and projected height from VRI were Unioned to TSA, WMU, LU and AWA boundaries in ArcGIS 10.6.

### Moose Population Status and Trends
Moose population data was obtained by searching the BC government's [species inventory web explorer (SIWE)](https://catalogue.data.gov.bc.ca/dataset/species-inventory-web-explorer-siwe) for 'moose' inventory data collected in the Okanagan or Thompson regions. I recorded data on moose density, populations, bull:cow ratios and calf:cow ratios.  

Recent moose density estimates in the Okanagan TSA ranged from 0.85 moose/km^2^ in WMU 8-11 to 0.22 moose/km^2^ in WMU 8-21 (Fig. 1). Most of these estimates were obtained between 2011 to 2019, with the exception of WMU 3-12 (0.26 moose/km^2^) and WMU 8.10 (0.27 moose/km^2^), which were obtained in 1985 and 1999, respectively. 

```{r, moose pop density map, fig.cap = "Figure 1. Most recent moose density estimates by wildlife management unit in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.moose.pop <- read.csv (file = here ("wildlife/data/okanagan/moose_density_20200504.csv"),
                            header = T)
data.wmu <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_wmu", 
                         quiet = T)
data.tsa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "tsa_outer", 
                         quiet = T)

data.wmu.moose <- merge (data.wmu, data.moose.pop, all = TRUE)
data.wmu.moose.dns <- data.wmu.moose %>%
                        filter (!is.na (density)) 
data.wmu.moose.dns <- data.wmu.moose.dns [data.wmu.moose.dns$WILDLIFE_MGMT_UNIT_ID != "4-32" &
                                          data.wmu.moose.dns$WILDLIFE_MGMT_UNIT_ID != "4-39" &
                                          data.wmu.moose.dns$WILDLIFE_MGMT_UNIT_ID != "3-42" &
                                          data.wmu.moose.dns$WILDLIFE_MGMT_UNIT_ID != "8-4" &
                                          data.wmu.moose.dns$WILDLIFE_MGMT_UNIT_ID != "8-15" , ]
data.wmu.moose.dns <- data.wmu.moose.dns [-c (5, 10, 12, 13, 17), ] # 10, 16, 18, 19, 24 
data.wmu.moose.dns$brks_density <- cut (data.wmu.moose.dns$density, 
                                        breaks = c (0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 
                                                    0.80, 0.90, 1.00, Inf), 
                                        labels = c ("0 - 0.09", "0.10 - 0.19", "0.20 - 0.29", "0.30 - 0.39",
                                                    "0.40 - 0.49", "0.50 - 0.59", "0.60 - 0.69", "0.70 - 0.79",
                                                    "0.80 - 0.89", "0.90 - 1.00", ">1.00"))

data.wmu.moose.dns <- cbind (data.wmu.moose.dns, st_coordinates (st_centroid (data.wmu.moose.dns)))

ggplot (data = data.wmu.moose.dns) +
        geom_sf (data = data.wmu.moose.dns, aes (fill = brks_density)) +
        geom_sf (data = data.tsa, fill = NA) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Density (moose/km2)") +
        geom_label(data = data.wmu.moose.dns, 
                   aes(X, Y, label = density), 
                   size = 2.5) +
        theme (axis.title.x = element_blank(),
               axis.title.y = element_blank())
```

The ratio of bull moose to cow moose is often used to indicate hunting pressure on moose populations, and a ratio of greater than 30 bulls to 100 cows is a typical management target, where populations below that indicate a heavily hunted, and potentially unstable population ( [Young and Boertje 2008](https://alcesjournal.org/index.php/alces/article/view/37); [Walker et al. 2017](http://a100.gov.bc.ca/pub/siwe/download.do?docId=34052)). 

Recent bull:cow ratios indicate that moose populations were under relatively high hunting pressure in the south and central portions of the TSA, ranging from 6 to 16 bulls:cows (Fig. 2). Ratios were highest in the eastern and west-central portions of the TSA, reaching up to 85 bulls:cows in WMU 8-14.

```{r, moose bull/cow ratio, fig.cap = "Figure 2. Most recent moose bull/cow ratios by wildlife management unit in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.moose.pop <- read.csv (file = here ("wildlife/data/okanagan/moose_density_20200504.csv"),
                            header = T)
data.wmu <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_wmu", 
                         quiet = T)
data.tsa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "tsa_outer", 
                         quiet = T)
data.wmu.moose <- merge (data.wmu, data.moose.pop, all = TRUE)
data.wmu.moose.bc <- data.wmu.moose %>%
                        filter (!is.na (bulls_100cows))
### FILTER OUT PARTIAL  WMUs ###
#### NOT 2 to 11, 32, 46 to 49, 62
data.wmu.moose.bc.2 <- data.wmu.moose.bc [c (1, 12:31, 33:45, 50:61), ] 
data.wmu.moose.bc.2 <- data.wmu.moose.bc.2 [c (1, 2, 8, 15, 19, 21, 24, 25, 27, 31, 32, 34, 35, 41, 45), ] 
# get the latest number
data.wmu.moose.bc.2$brks_bc <- cut (data.wmu.moose.bc.2$bulls_100cows, 
                                        breaks = c (0, 15, 30, 45, 60, 75, 90, Inf), 
                                        labels = c ("0 - 14", "15 - 29", "30 - 44", "45 - 59",
                                                    "60 - 74", "75 - 89", ">90"))
data.wmu.moose.bc.2 <- cbind (data.wmu.moose.bc.2, st_coordinates (st_centroid (data.wmu.moose.bc.2)))

ggplot (data = data.wmu.moose.bc.2) +
        geom_sf (data = data.wmu.moose.bc.2, aes (fill = brks_bc)) +
        geom_sf (data = data.tsa, fill = NA) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Bulls:100 Cows") +
        geom_label(data = data.wmu.moose.bc.2, 
                   aes(X, Y, label = bulls_100cows), 
                   size = 2.5) +
        theme (axis.title.x = element_blank(),
               axis.title.y = element_blank())

```

The ratio of calves to cow moose is often used as an indicator of moose population trend, where ratios of 25 to 30 indicate a stable population, and ratios greater than 30 indicate an increasing population ( [FLNRORD 2019](http://www.env.gov.bc.ca/fw/wildlife/management-issues/docs/2019%20Moose%20Factsheet.pdf)). 

The most recent calf:cow ratios were lowest in the north-central and south-central portions of the TSA, and highest in the central and far southern portions of the TSA (Fig. 3). The lowest calf:cow ratios (19) were less than 25, indicating a potential decreasing population in those areas. However, the majority of WMUs were close to or above 30 calves:100 cows, indicating a high potential for stable to increasing population trends for much of the TSA.

```{r, calf/cow ratio, fig.cap = "Figure 3. Most recent moose calf/cow ratios by wildlife management unit in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.moose.pop <- read.csv (file = here ("wildlife/data/okanagan/moose_density_20200504.csv"),
                            header = T)
data.wmu <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_wmu", 
                         quiet = T)
data.tsa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "tsa_outer", 
                         quiet = T)
data.wmu.moose <- merge (data.wmu, data.moose.pop, all = TRUE)
data.wmu.moose.bc <- data.wmu.moose %>%
                        filter (!is.na (calves_100cows))
### FILTER OUT PARTIAL  WMUs ###
#### NOT 2 to 11, 32, 46 to 49, 62
data.wmu.moose.bc.2 <- data.wmu.moose.bc [c (1, 12:31, 33:45, 50:61), ] 
data.wmu.moose.bc.2 <- data.wmu.moose.bc.2 [c (1, 2, 8, 15, 19, 21, 24, 25, 27, 31, 32, 34, 35, 41, 45), ] 
data.wmu.moose.bc.2$brks_bc <- cut (data.wmu.moose.bc.2$calves_100cows, 
                                        breaks = c (0, 10, 20, 30, 40, 50, 60, 70, Inf), 
                                        labels = c ("0 - 9", "10 - 19", "20 - 29", "30 - 39",
                                                    "40 - 49", "50 - 59", "60 - 69", ">70"))
data.wmu.moose.bc.2 <- cbind (data.wmu.moose.bc.2, st_coordinates (st_centroid (data.wmu.moose.bc.2)))

ggplot (data = data.wmu.moose.bc.2) +
        geom_sf (data = data.wmu.moose.bc.2, aes (fill = brks_bc)) +
        geom_sf (data = data.tsa, fill = NA) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Calves:100 Cows") +
        geom_label(data = data.wmu.moose.bc.2, 
                   aes(X, Y, label = calves_100cows), 
                   size = 2.5) +
        theme (axis.title.x = element_blank(),
               axis.title.y = element_blank())

```

Recently, the moose population appears to have been increasing across the Okanagan region. In 2007, the population was estimated at 2,200 animals ( [Gyug 2007](http://a100.gov.bc.ca/pub/siwe/details.do?projectId=4279&surveyId=10889&pagerOffset=0)), then it was estimated at 3,913 animals in 2013 and 4,352 in 2017 ( [Walker et al. 2017](http://a100.gov.bc.ca/pub/siwe/details.do?projectId=4279&surveyId=36576&pagerOffset=0)). However, within the Okanagan TSA there were some WMUs where moose populations appeared to be decreasing (Fig. 4). Moose densities decreased from 0.38 moose/km^2^ to 0.31 moose/km^2^ and 0.38 moose/km^2^ to 0.26 moose/km^2^ in WMUs 8-23 and 8-24, respectively, between 2003 and 2014. 

```{r, moose pop density trend, fig.cap = "Figure 4. Trends in moose density in wildlife management units with more than one recent density estimate within the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.moose.pop <- read.csv (file = here ("wildlife/data/okanagan/moose_density_20200504.csv"),
                            header = T)
data.wmu <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_wmu", 
                         quiet = T)
data.wmu.moose <- merge (data.wmu, data.moose.pop, all = TRUE)
data.wmu.moose <- data.wmu.moose [- c(74:77), ]
data.wmu.moose.trend <- data.wmu.moose %>%
                          filter (# WILDLIFE_MGMT_UNIT_ID == "4-39" | WILDLIFE_MGMT_UNIT_ID == "8-1" | 
                                  WILDLIFE_MGMT_UNIT_ID == "8-6" | # WILDLIFE_MGMT_UNIT_ID == "8-8" | 
                                  # WILDLIFE_MGMT_UNIT_ID == "8-9" | 
                                  # WILDLIFE_MGMT_UNIT_ID == "8-10" | 
                                  WILDLIFE_MGMT_UNIT_ID == "8-11" | # WILDLIFE_MGMT_UNIT_ID == "8-12" | 
                                  # WILDLIFE_MGMT_UNIT_ID == "8-14" | 
                                  # WILDLIFE_MGMT_UNIT_ID == "8-21" | 
                                  WILDLIFE_MGMT_UNIT_ID == "8-23" | WILDLIFE_MGMT_UNIT_ID == "8-24"  
                                  # WILDLIFE_MGMT_UNIT_ID == "8-25"
                                  ) # these units have >1 record for density
data.wmu.moose.trend$year <- as.numeric (data.wmu.moose.trend$year)

ggplot (data = data.wmu.moose.trend,
        aes (x = year, y = density)) +
  geom_point () +
  geom_smooth (method = lm) +
  facet_grid (~ WILDLIFE_MGMT_UNIT_ID) +
  theme_bw () +
  theme (axis.text.x = element_text (angle = 90)) +
  scale_y_continuous (limits = c (0, 1),
                      breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  ylab ("Density (moose/km2)") + 
  xlab ("Year")
  
```

### Road Density
Current road density (km/km^2^) estimates were relatively high (greater than 2 km/km^2^) across much of the Okanagan TSA (Fig. 5). Lower road densities, less than 1 km/km^2^, typically occurred in the northeast (Monashee mountains and towards Revelstoke), and southwest (EC Manning and Cathedral provincial parks) portions of the TSA. 

```{r, calc current road density, eval = F, echo = F, message = F, warning = F, include = F}
data.roads <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                           layer = "union_roads_final", 
                           quiet = T)
data.roads <- as.data.frame (data.roads)
data.road.length <- as.data.table (data.roads %>%
                                    group_by (WATERSHED_FEATURE_ID) %>%
                                    summarise (sum_length_km = sum (Length_km)))
data.awa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_fawa", 
                         quiet = T)
data.awa.roads <- merge (data.awa, data.road.length, all = TRUE)
data.awa.roads$road_density <- data.awa.roads$sum_length_km / (data.awa.roads$AREA_HA*0.01) # ha to km2
st_write (data.awa.roads, 
          dsn = here ("wildlife/data/okanagan/data_awa_roads.shp"))
```

```{r, map current road density, fig.cap = "Figure 5. Map of current road density by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.roads <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_roads.shp"),
                               quiet = T)
data.awa.roads$rd_dn [is.na (data.awa.roads$rd_dn)] <- 0
data.awa.roads$brks_rd_dn <- cut (data.awa.roads$rd_dn, 
                                  breaks = c (-0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, Inf), 
                                  labels = c ("0 - 0.19", "0.20 - 0.39", "0.40 - 0.59", "0.60 - 0.79",
                                              "0.80 - 0.99", "1.00 - 1.19", "1.20 - 1.39", "1.40 - 1.59",
                                              "1.60 - 1.79", "1.80 - 1.99", ">2.00"))
ggplot (data = data.awa.roads) +
        geom_sf (data = data.awa.roads, aes (fill = brks_rd_dn)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Road Density (km/km2)") 

```

### Forest Cover for Moose
The AWA's with the highest proportion of conifer forest cover for moose (i.e., stands at least 5 ha in size and 15 m tall) occurred in the west-central portion of the TSA, west of Okanagan Lake, north of Shuswap Lake and north of Osoyoos (Fig. 6). There were a limited number of AWA's with little amounts of conifer forest cover along the western edge of the TSA. AWA's with an approximately 40% proportion of conifer forest cover occurred in the east-central and north-west portions of the TSA. 

```{r, calc area of conifer stands >15m, eval = F, echo = F, message = F, warning = F, include = F}
data.vri <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "union_vri_lean", 
                         quiet = T)
data.vri <- as.data.frame (data.vri)
data.area.height <- as.data.table (data.vri %>%
                                    filter (SPECIES_CD_1 == "B" | SPECIES_CD_1 == "BA" | SPECIES_CD_1 == "BB" |
                                            SPECIES_CD_1 == "BG" | SPECIES_CD_1 == "BL" | SPECIES_CD_1 == "BN" |
                                            SPECIES_CD_1 == "C" | SPECIES_CD_1 == "CW" | SPECIES_CD_1 == "F" | 
                                            SPECIES_CD_1 == "FD" | SPECIES_CD_1 == "FDC" | SPECIES_CD_1 == "FDI" | 
                                            SPECIES_CD_1 == "H" | SPECIES_CD_1 == "HM" | SPECIES_CD_1 == "HW" | 
                                            SPECIES_CD_1 == "HX" | SPECIES_CD_1 == "HXM" | SPECIES_CD_1 == "P" | 
                                            SPECIES_CD_1 == "PA" | SPECIES_CD_1 == "PF" | SPECIES_CD_1 == "PJ" | 
                                            SPECIES_CD_1 == "PL" | SPECIES_CD_1 == "PLC" | SPECIES_CD_1 == "PLI" | 
                                            SPECIES_CD_1 == "PR" | SPECIES_CD_1 == "PW" | SPECIES_CD_1 == "PX" | 
                                            SPECIES_CD_1 == "PXJ" | SPECIES_CD_1 == "PY" | SPECIES_CD_1 == "S" | 
                                            SPECIES_CD_1 == "SA" | SPECIES_CD_1 == "SB" | SPECIES_CD_1 == "SE" | 
                                            SPECIES_CD_1 == "SS" | SPECIES_CD_1 == "SW" | SPECIES_CD_1 == "SX" | 
                                            SPECIES_CD_1 == "SXB" | SPECIES_CD_1 == "SXE" | SPECIES_CD_1 == "SXL" | 
                                            SPECIES_CD_1 == "SXS" | SPECIES_CD_1 == "SXW" | SPECIES_CD_1 == "SXX" | 
                                           SPECIES_CD_1 == "YC") %>%
                                    filter (PROJ_HEIGHT_1 >= 15) %>%
                                    filter (AREA_HA >= 5) %>%
                                    group_by (WATERSHED_FEATURE_ID) %>%
                                    summarise (sum_area_ha = sum (AREA_HA)))
data.awa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_fawa", 
                         quiet = T)
data.awa.height <- merge (data.awa, data.area.height, all = TRUE)
data.awa.height$prop_moose_cover <- data.awa.height$sum_area_ha / data.awa.height$AREA_HA # proportion moose cover
st_write (data.awa.height, 
          dsn = here ("wildlife/data/okanagan/data_awa_cover.shp"))
```

```{r, map current forest cover, fig.cap = "Figure 6. Map of current percentage of conifer forest cover greater than 14 m tall and 5 ha patch size by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cover <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_cover.shp"),
                               quiet = T)
data.awa.cover$prp_ [is.na (data.awa.cover$prp_)] <- 0
data.awa.cover$brks_prop_moose <- cut (data.awa.cover$prp_, 
                                       breaks = c (-0.1, 0.2, 0.4, 0.6, 0.8, Inf), 
                                       labels = c ("0 - 0.19", "0.20 - 0.39", "0.40 - 0.59", "0.60 - 0.79",
                                                   "0.80 - 1.00"))
ggplot (data = data.awa.cover) +
        geom_sf (data = data.awa.cover, aes (fill = brks_prop_moose)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Proportion of Conifer Cover 
> 14m tall and 5 ha in Size") 
```

### Early Seral Forest for Moose
The amount of early seral forest (5 to 30 years old) from cutblocks was highest in the west central and southeast parts of the TSA (Fig. 7). Some AWA's in these areas exceeded a proportion of 0.40 early seral cutblocks. The majority of AWA's were less than 0.20 early seral cutblocks.

```{r, calc area of cutblocks 5 to 30 yo, eval = F, echo = F, message = F, warning = F, include = F}
data.cutblocks <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                               layer = "union_cutblocks", 
                               quiet = T)
data.cutblocks <- as.data.frame (data.cutblocks)
data.cut.5.30 <- as.data.table (data.cutblocks %>%
                                    filter (HARVEST_YEAR >= 1988 & HARVEST_YEAR <= 2013) %>% # cutblocks 5 to 30 yo
                                    group_by (WATERSHED_FEATURE_ID) %>%
                                    summarise (sum_area_ha = sum (AREA_HA)))
data.cut.5.30 <- data.cut.5.30 [-1, ]
data.awa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_fawa", 
                         quiet = T)
data.awa.cut <- merge (data.awa, data.cut.5.30, all = TRUE)
data.awa.cut$prop_moose_food <- data.awa.cut$sum_area_ha / data.awa.cut$AREA_HA # proportion moose forage (5 to 30 yo cut)
st_write (data.awa.cut, 
          dsn = here ("wildlife/data/okanagan/data_awa_cut5to30.shp"))
```

```{r, map cutblocks 5 to 30 yo, fig.cap = "Figure 7. Map of proportion of area of cutblocks 5 to 30 years old by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cut <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_cut5to30.shp"),
                             quiet = T)
data.awa.cut$prp__ [is.na (data.awa.cut$prp__)] <- 0
data.awa.cut$brks_prop_cut <- cut (data.awa.cut$prp__, 
                                     breaks = c (-0.1, 0.2, 0.4, 0.6, 0.8, 1), 
                                     labels = c ("0 - 0.19", "0.20 - 0.39", "0.40 - 0.59", "0.60 - 0.79", 
                                                 "0.80 - 1.00"))
ggplot (data = data.awa.cut, aes (fill = brks_prop_cut)) +
        geom_sf (data = data.awa.cut, aes (fill = brks_prop_cut)) +
       # scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Proportion of 5 to 30 year 
old cutblocks") 
```

## 5. Moose Habitat and Population Management Tools 
Previous research suggests that moose require a diversity of forest stand types to meet their forage and cover needs, for example, a mix of older and younger, open and more dense and conifer and deciduous stands. This can make habitat management for moose challenging, as it is difficult to develop and apply a simple forest management regime to meet all of these needs and maximize moose habitat. 

Roads are typically found to have a negative effect on moose, by increasing human accessiblity into moose habitat, and thus hunting pressure and disturbance that can negatively affect moose survival. Therefore, a management regime to minimize road density may be useful for minimizing forestry impacts on moose. However, moose populations across the Okanagan TSA seem to overall be stable or increasing, despite relatively high (greater than 2 km/km^2^) road densities.   

Currently, there are no alternative proposed large-scale forest management regimes designed to maximize moose habitat quality in the Okanagan TSA, and no simulated alternative regimes were applied within the TSR model. 

## 6. Simulated Future States of Moose Habitat and Populations Under Different Management Scenarios in the Timber Supply Area 
### Base Case Scenario
#### Simulated Road Density
Future forestry road density in the TSA was simulated as a function of simulated future cutblock density as outputted from the TSR model, based on a modeled statistical relationship between cutblock density and road density in AWA's ( [Muhly 2016](https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/stewardship/forest-analysis-inventory/tsr-annual-allowable-cut/wildlife-analysis/road_cutblock_model_20161102.pdf)).

Under the base case scenario, simulated road density increased steadily but at small increments from 2018 to 2078 across the TSA (Fig. 8). The median road density in AWA's in 2018 was 2.45 km/km^2^, which increased to 2.66 km/km^2^ by 2078. The increase in road density primarily occurred in a limited number of AWA's mostly along the edges of low road density areas in the northeast and southwest potions of the TSA (Fig. 9, Fig. 10). 

```{r, calc area new cut and new roads by decade, eval = F, echo = F, message = F, warning = F}
## current cutblocks
data.cutblocks <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                               layer = "union_cutblocks", 
                               quiet = T)
data.cutblocks <- as.data.frame (data.cutblocks)
cut.start <- as.data.table (data.cutblocks %>%
                              group_by (WATERSHED_FEATURE_ID) %>%
                                summarise (sum_area_ha = sum (AREA_HA)))
cut.start <- cut.start [-1, ]

## simulated future cutblocks, new area, not overlapping existing cb's
data.future.cutblocks <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                                      layer = "union_future_cutblocks_new_area", 
                                      quiet = T)
data.future.cutblocks <- data.frame (data.future.cutblocks)

# annual new area cut
cut.2028 <- data.future.cutblocks %>%
              filter (log_2027 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))  
cut.2028 <- cut.2028 [-1, ]

cut.2038 <- data.future.cutblocks %>%
              filter (log_2037 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2038 <- cut.2038 [-1, ]

cut.2048 <- data.future.cutblocks %>%
              filter (log_2047 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2048 <- cut.2048 [-1, ]

cut.2058 <- data.future.cutblocks %>%
              filter (log_2057 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2058 <- cut.2058 [-1, ]

cut.2068 <- data.future.cutblocks %>%
              filter (log_2067 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2068 <- cut.2068 [-1, ]

cut.2078 <- data.future.cutblocks %>%
              filter (log_2077 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2078 <- cut.2078 [-1, ]

cut.2088 <- data.future.cutblocks %>%
              filter (log_2087 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2088 <- cut.2088 [-1, ]

cut.2098 <- data.future.cutblocks %>%
              filter (log_2097 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2098 <- cut.2098 [-1, ]

cut.2108 <- data.future.cutblocks %>%
              filter (log_2107 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2108 <- cut.2108 [-1, ]

cut.2118 <- data.future.cutblocks %>%
              filter (log_2117 == 1) %>%
                group_by (WATERSHED_FEATURE_ID) %>%
                  summarise (area.ha = sum (area_ha_1, na.rm = T))                  
cut.2118 <- cut.2118 [-1, ]

# join new are cut to AWA data
data.awa <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                         layer = "clip_fawa")
data.awa <- data.awa %>%
              select (WATERSHED_FEATURE_ID, AREA_HA, GEOMETRY)

data.awa.cut.decade <- data.awa %>% 
                          full_join (cut.start, by = "WATERSHED_FEATURE_ID") %>%
                          rename (all_area_ha_cut_2018 = sum_area_ha) %>%
                            full_join (cut.2028, by = "WATERSHED_FEATURE_ID") %>%
                            rename (new_area_ha_cut_2028 = area.ha) %>%
                             full_join (cut.2038, by = "WATERSHED_FEATURE_ID") %>%
                             rename (new_area_ha_cut_2038 = area.ha) %>%
                              full_join (cut.2048, by = "WATERSHED_FEATURE_ID") %>%
                              rename (new_area_ha_cut_2048 = area.ha) %>%
                                full_join (cut.2058, by = "WATERSHED_FEATURE_ID") %>%
                                rename (new_area_ha_cut_2058 = area.ha) %>%
                                  full_join (cut.2068, by = "WATERSHED_FEATURE_ID") %>%
                                  rename (new_area_ha_cut_2068 = area.ha) %>%
                                      full_join (cut.2078, by = "WATERSHED_FEATURE_ID") %>%
                                      rename (new_area_ha_cut_2078 = area.ha) %>%
                                        full_join (cut.2088, by = "WATERSHED_FEATURE_ID") %>%
                                        rename (new_area_ha_cut_2088 = area.ha) %>%
                                          full_join (cut.2098, by = "WATERSHED_FEATURE_ID") %>%
                                          rename (new_area_ha_cut_2098 = area.ha) %>%
                                            full_join (cut.2108, by = "WATERSHED_FEATURE_ID") %>%
                                            rename (new_area_ha_cut_2108 = area.ha) %>%
                                              full_join (cut.2118, by = "WATERSHED_FEATURE_ID") %>%
                                              rename (new_area_ha_cut_2118 = area.ha)
data.awa.cut.decade [is.na (data.awa.cut.decade)] <- 0

# from new area calc new proportion/density cut   
data.awa.cut.decade$prop_all_cut_2018 <- data.awa.cut.decade$all_area_ha_cut_2018 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2028 <- data.awa.cut.decade$new_area_ha_cut_2028 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2038 <- data.awa.cut.decade$new_area_ha_cut_2038 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2048 <- data.awa.cut.decade$new_area_ha_cut_2048 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2058 <- data.awa.cut.decade$new_area_ha_cut_2058 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2068 <- data.awa.cut.decade$new_area_ha_cut_2068 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2078 <- data.awa.cut.decade$new_area_ha_cut_2078 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2088 <- data.awa.cut.decade$new_area_ha_cut_2088 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2098 <- data.awa.cut.decade$new_area_ha_cut_2098 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2108 <- data.awa.cut.decade$new_area_ha_cut_2108 / data.awa.cut.decade$AREA_HA 
data.awa.cut.decade$prop_new_cut_2118 <- data.awa.cut.decade$new_area_ha_cut_2118 / data.awa.cut.decade$AREA_HA 

# add current road density data to table
data.roads <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                               layer = "union_roads_final")
data.roads <- as.data.frame (data.roads)
roads.start <- as.data.table (data.roads %>%
                                group_by (WATERSHED_FEATURE_ID) %>%
                                  summarise (road_length_km_start = sum (Length_km)))
data.awa.cut.decade <- data.awa.cut.decade %>% 
                          full_join (roads.start, by = "WATERSHED_FEATURE_ID") 
data.awa.cut.decade$road_density_start <- data.awa.cut.decade$road_length_km_start / (data.awa.cut.decade$AREA_HA * 0.01)

# use road-cutblock model to esitmate new road density from new cut density
data.awa.cut.decade$road_density_2028 <- data.awa.cut.decade$road_density_start + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2028)
# coefficients from https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/stewardship/forest-analysis-inventory/tsr-annual-allowable-cut/wildlife-analysis/road_cutblock_model_20161102.pdf; used fixed effect + unit random effect
# don't use intercept because using current road density in awa as 'intercept'
data.awa.cut.decade$road_density_2038 <- data.awa.cut.decade$road_density_2028 + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2038)
data.awa.cut.decade$road_density_2048 <- data.awa.cut.decade$road_density_2038 + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2048)
data.awa.cut.decade$road_density_2058 <- data.awa.cut.decade$road_density_2048 + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2058)
data.awa.cut.decade$road_density_2068 <- data.awa.cut.decade$road_density_2058 + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2068)
data.awa.cut.decade$road_density_2078 <- data.awa.cut.decade$road_density_2068 + 
                                          ((1.16 - 0.327) * data.awa.cut.decade$prop_new_cut_2078)

# calculate estimate of area of 5 to 30 yo cutblocks, by decade
cut.10.1998 <- as.data.table (data.cutblocks %>% # 
                                  group_by (WATERSHED_FEATURE_ID) %>%
                                    filter (HARVEST_YEAR >= 1989 & HARVEST_YEAR <= 1998) %>% # cutblocks in decade
                                      summarise (cut_10_1998_area_ha = sum (AREA_HA)))
cut.10.1998 <- cut.10.1998 [-1, ]  

cut.10.2008 <- as.data.table (data.cutblocks %>% 
                                  group_by (WATERSHED_FEATURE_ID) %>%
                                    filter (HARVEST_YEAR >= 1999 & HARVEST_YEAR <= 2008) %>% 
                                      summarise (cut_10_2008_area_ha = sum (AREA_HA)))
cut.10.2008 <- cut.10.2008 [-1, ]  

cut.10.2018 <- as.data.table (data.cutblocks %>%
                                  group_by (WATERSHED_FEATURE_ID) %>%
                                    filter (HARVEST_YEAR >= 2009 & HARVEST_YEAR <= 2018) %>% 
                                      summarise (cut_10_2018_area_ha = sum (AREA_HA)))
cut.10.2018 <- cut.10.2018 [-1, ]  

total.0to30cut.2018 <- as.data.table (data.cutblocks %>%
                                        group_by (WATERSHED_FEATURE_ID) %>%
                                          filter (HARVEST_YEAR >= 1989 & HARVEST_YEAR <= 2018) %>% 
                                            summarise (total_0to30cut_2018 = sum (AREA_HA)))
total.0to30cut.2018 <- total.0to30cut.2018 [-1, ]  

data.awa.cut.decade <- data.awa.cut.decade %>% 
                          full_join (cut.10.1998, by = "WATERSHED_FEATURE_ID") %>%
                            full_join (cut.10.2008, by = "WATERSHED_FEATURE_ID") %>%
                             full_join (cut.10.2018, by = "WATERSHED_FEATURE_ID") %>%
                              full_join (total.0to30cut.2018, by = "WATERSHED_FEATURE_ID")
                              
data.awa.cut.decade [is.na (data.awa.cut.decade)] <- 0

## simulated future cutblocks, new area, not overlapping existing cb's
data.future.cutblocks.all <- sf::st_read (dsn = "C:\\Work\\tsr\\data\\okanagan_tsr_20200407.gdb",
                                          layer = "union_future_cutblocks")
data.future.cutblocks.all <- data.frame (data.future.cutblocks.all)

cut.2028.all <- data.future.cutblocks.all %>%
                  filter (log_2027 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2028_area_ha = sum (area_ha_1, na.rm = T))  
cut.2028.all <- cut.2028.all [-1, ]

cut.2038.all <- data.future.cutblocks.all %>%
                  filter (log_2037 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2038_area_ha = sum (area_ha_1, na.rm = T))  
cut.2038.all <- cut.2038.all [-1, ]

cut.2048.all <- data.future.cutblocks.all %>%
                  filter (log_2047 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2048_area_ha = sum (area_ha_1, na.rm = T))  
cut.2048.all <- cut.2048.all [-1, ]

cut.2058.all <- data.future.cutblocks.all %>%
                  filter (log_2057 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2058_area_ha = sum (area_ha_1, na.rm = T))  
cut.2058.all <- cut.2058.all [-1, ]

cut.2068.all <- data.future.cutblocks.all %>%
                  filter (log_2067 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2068_area_ha = sum (area_ha_1, na.rm = T))  
cut.2068.all <- cut.2068.all [-1, ]

cut.2078.all <- data.future.cutblocks.all %>%
                  filter (log_2077 == 1) %>%
                    group_by (WATERSHED_FEATURE_ID) %>%
                      summarise (cut_all_2078_area_ha = sum (area_ha_1, na.rm = T))  
cut.2078.all <- cut.2078.all [-1, ]

data.awa.cut.decade <- data.awa.cut.decade %>% 
                          full_join (cut.2028.all, by = "WATERSHED_FEATURE_ID") %>%
                            full_join (cut.2038.all, by = "WATERSHED_FEATURE_ID") %>%
                             full_join (cut.2048.all, by = "WATERSHED_FEATURE_ID") %>%
                              full_join (cut.2058.all, by = "WATERSHED_FEATURE_ID") %>%
                                full_join (cut.2068.all, by = "WATERSHED_FEATURE_ID") %>%
                                  full_join (cut.2078.all, by = "WATERSHED_FEATURE_ID")
data.awa.cut.decade [is.na (data.awa.cut.decade)] <- 0
  
data.awa.cut.decade$total_0to30cut_2028 <- data.awa.cut.decade$total_0to30cut_2018 + 
                                           data.awa.cut.decade$cut_all_2028_area_ha -
                                           data.awa.cut.decade$cut_10_1998_area_ha
data.awa.cut.decade$total_0to30cut_2028 [data.awa.cut.decade$total_0to30cut_2028 < 0] <- 0 # some values slightly <0
  
data.awa.cut.decade$total_0to30cut_2038 <- data.awa.cut.decade$total_0to30cut_2028 + 
                                           data.awa.cut.decade$cut_all_2038_area_ha -
                                           data.awa.cut.decade$cut_10_2008_area_ha
data.awa.cut.decade$total_0to30cut_2038 [data.awa.cut.decade$total_0to30cut_2038 < 0] <- 0 # some values slightly <0

data.awa.cut.decade$total_0to30cut_2048 <- data.awa.cut.decade$total_0to30cut_2038 + 
                                           data.awa.cut.decade$cut_all_2048_area_ha -
                                           data.awa.cut.decade$cut_10_2018_area_ha
data.awa.cut.decade$total_0to30cut_2048 [data.awa.cut.decade$total_0to30cut_2048 < 0] <- 0 # some values slightly <0

data.awa.cut.decade$total_0to30cut_2058 <- data.awa.cut.decade$total_0to30cut_2048 + 
                                           data.awa.cut.decade$cut_all_2058_area_ha -
                                           data.awa.cut.decade$cut_all_2028_area_ha  
data.awa.cut.decade$total_0to30cut_2058 [data.awa.cut.decade$total_0to30cut_2058 < 0] <- 0 # some values slightly <0

data.awa.cut.decade$total_0to30cut_2068 <- data.awa.cut.decade$total_0to30cut_2058 + 
                                           data.awa.cut.decade$cut_all_2068_area_ha -
                                           data.awa.cut.decade$cut_all_2038_area_ha  
data.awa.cut.decade$total_0to30cut_2068 [data.awa.cut.decade$total_0to30cut_2068 < 0] <- 0 # some values slightly <0

data.awa.cut.decade$total_0to30cut_2078 <- data.awa.cut.decade$total_0to30cut_2068 + 
                                           data.awa.cut.decade$cut_all_2078_area_ha -
                                           data.awa.cut.decade$cut_all_2048_area_ha  
data.awa.cut.decade$total_0to30cut_2078 [data.awa.cut.decade$total_0to30cut_2078 < 0] <- 0 # some values slightly <0


data.awa.cut.decade$prop_0t30_2018 <- data.awa.cut.decade$total_0to30cut_2018/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2028 <- data.awa.cut.decade$total_0to30cut_2028/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2038 <- data.awa.cut.decade$total_0to30cut_2038/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2048 <- data.awa.cut.decade$total_0to30cut_2048/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2058 <- data.awa.cut.decade$total_0to30cut_2058/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2068 <- data.awa.cut.decade$total_0to30cut_2068/data.awa.cut.decade$AREA_HA
data.awa.cut.decade$prop_0t30_2078 <- data.awa.cut.decade$total_0to30cut_2078/data.awa.cut.decade$AREA_HA

# some clean up
data.awa.cut.decade$new_area_ha_cut_2028 <- NULL
data.awa.cut.decade$new_area_ha_cut_2038 <- NULL
data.awa.cut.decade$new_area_ha_cut_2038 <- NULL
data.awa.cut.decade$new_area_ha_cut_2048 <- NULL
data.awa.cut.decade$new_area_ha_cut_2058 <- NULL
data.awa.cut.decade$new_area_ha_cut_2068 <- NULL
data.awa.cut.decade$new_area_ha_cut_2078 <- NULL
data.awa.cut.decade$new_area_ha_cut_2088 <- NULL
data.awa.cut.decade$new_area_ha_cut_2098 <- NULL
data.awa.cut.decade$new_area_ha_cut_2108 <- NULL
data.awa.cut.decade$new_area_ha_cut_2118 <- NULL
data.awa.cut.decade$prop_new_cut_2028 <- NULL
data.awa.cut.decade$prop_new_cut_2038 <- NULL
data.awa.cut.decade$prop_new_cut_2048 <- NULL
data.awa.cut.decade$prop_new_cut_2058 <- NULL
data.awa.cut.decade$prop_new_cut_2068 <- NULL
data.awa.cut.decade$prop_new_cut_2078 <- NULL
data.awa.cut.decade$prop_new_cut_2088 <- NULL
data.awa.cut.decade$prop_new_cut_2098 <- NULL
data.awa.cut.decade$prop_new_cut_2108 <- NULL
data.awa.cut.decade$prop_new_cut_2118 <- NULL
data.awa.cut.decade$road_length_km_start <- NULL
data.awa.cut.decade$cut_10_1998_area_ha <- NULL
data.awa.cut.decade$cut_10_2008_area_ha <- NULL
data.awa.cut.decade$cut_10_2018_area_ha <- NULL
data.awa.cut.decade$total_0to30cut_2018 <- NULL
data.awa.cut.decade$cut_all_2028_area_ha <- NULL
data.awa.cut.decade$cut_all_2038_area_ha <- NULL
data.awa.cut.decade$cut_all_2048_area_ha <- NULL
data.awa.cut.decade$cut_all_2058_area_ha <- NULL
data.awa.cut.decade$cut_all_2068_area_ha <- NULL
data.awa.cut.decade$cut_all_2078_area_ha <- NULL
data.awa.cut.decade$total_0to30cut_2028 <- NULL
data.awa.cut.decade$total_0to30cut_2038 <- NULL
data.awa.cut.decade$total_0to30cut_2048 <- NULL
data.awa.cut.decade$total_0to30cut_2058 <- NULL
data.awa.cut.decade$total_0to30cut_2068 <- NULL
data.awa.cut.decade$total_0to30cut_2078 <- NULL

st_write (data.awa.cut.decade, 
          dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"))
```

```{r, boxplot current and sim road density, fig.cap = "Figure 8. Current and simulated future road density by decade in freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.road <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.road <- as.data.table (data.awa.road)
data.awa.road.box <- melt (data.awa.road,  
                           id.vars = c("WATERSH"),
                           measure.vars = c("rd_dns_", "r__2028", "r__2038", "r__2048", "r__2058", 
                                            "r__2068", "r__2078"),
                           variable.name = "decade", 
                           value.name = "road_density")

data.awa.road.box [decade == "rd_dns_", decade := "2018"]
data.awa.road.box [decade == "r__2028", decade := "2028"]
data.awa.road.box [decade == "r__2038", decade := "2038"]
data.awa.road.box [decade == "r__2048", decade := "2048"]
data.awa.road.box [decade == "r__2058", decade := "2058"]
data.awa.road.box [decade == "r__2068", decade := "2068"]
data.awa.road.box [decade == "r__2078", decade := "2078"]

ggplot (data = data.awa.road.box,
        aes (x = decade, y = road_density)) +
        geom_boxplot () + 
        scale_y_continuous (name = "Road Density (km/km2)", 
                            breaks = seq (0, 14, 2)) +
        scale_x_discrete (name = "Decade") +
        coord_flip () +
        geom_hline (yintercept = 2.45, colour = "green")

```

```{r, map sim road density 2048, fig.cap = "Figure 9. Map of simulated road density in 2048 by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.road <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.road$brks_rd_dn_2048 <- cut (data.awa.road$r__2048, 
                                      breaks = c (-0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, Inf), 
                                      labels = c ("0 - 0.19", "0.20 - 0.39", "0.40 - 0.59", "0.60 - 0.79",
                                                  "0.80 - 0.99", "1.00 - 1.19", "1.20 - 1.39", "1.40 - 1.59",
                                                  "1.60 - 1.79", "1.80 - 1.99", ">2.00"))
ggplot (data = data.awa.road) +
        geom_sf (data = data.awa.road, aes (fill = brks_rd_dn_2048)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Road Density (km/km2)") 
```

```{r, map sim road density 2078, fig.cap = "Figure 10. Map of simulated road density in 2078 by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.road <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.road$brks_rd_dn_2078 <- cut (data.awa.road$r__2068, 
                                      breaks = c (-0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, Inf), 
                                      labels = c ("0 - 0.19", "0.20 - 0.39", "0.40 - 0.59", "0.60 - 0.79",
                                                  "0.80 - 0.99", "1.00 - 1.19", "1.20 - 1.39", "1.40 - 1.59",
                                                  "1.60 - 1.79", "1.80 - 1.99", ">2.00"))
ggplot (data = data.awa.road) +
        geom_sf (data = data.awa.road, aes (fill = brks_rd_dn_2078)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Road Density (km/km2)") 
```

#### Simulated Early Seral Forest
The median density of 0 to 30 year old cutblocks in AWA's was 0.09 km^2^/km^2^ in 2018. This increased slightly to 0.10 km^2^/km^2^ in 2028 and then declined to approximately 0.06 km^2^/km^2^ in subsequent decades (Fig. 11). Cutblock density became more disbursed from 2018 to 2078 (i.e., there were more lower density cutblock AWA's and fewer high density AWA's). The location of cutblocks shifted from the west central portion of the TSA (Fig. 11), to east central portions in 2014 (Fig. 12) and northern portions (Fig. 13) of the TSA in 2078. 

```{r, boxplot current and sim 0 to 30 cutblock density, fig.cap = "Figure 11. Distribtuion of current cutblock density by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cut <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.cut <- as.data.table (data.awa.cut)
data.awa.cut.box <- melt (data.awa.cut,  
                           id.vars = c("WATERSH"),
                           measure.vars = c("p_030_201", "p_030_202", "p_030_203", "p_030_204", "p_030_205", 
                                            "p_030_206", "p_030_207"),
                           variable.name = "decade", 
                           value.name = "cutblock_density")

data.awa.cut.box [decade == "p_030_201", decade := "2018"]
data.awa.cut.box [decade == "p_030_202", decade := "2028"]
data.awa.cut.box [decade == "p_030_203", decade := "2038"]
data.awa.cut.box [decade == "p_030_204", decade := "2048"]
data.awa.cut.box [decade == "p_030_205", decade := "2058"]
data.awa.cut.box [decade == "p_030_206", decade := "2068"]
data.awa.cut.box [decade == "p_030_207", decade := "2078"]

ggplot (data = data.awa.cut.box,
        aes (x = decade, y = cutblock_density)) +
        geom_boxplot () + 
        scale_y_continuous (name = "Cutblock Density (km2/km2)", 
                            breaks = seq (0, 1, 0.1)) +
        scale_x_discrete (name = "Decade") +
        coord_flip () +
        geom_hline (yintercept = 0.09, colour = "green")

```

```{r, map sim cutblock density 2018, fig.cap = "Figure 12. Map of simulated cutblock density in 2018 by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cut <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.cut$brks_cut_2018 <- cut (data.awa.cut$p_030_201, 
                                   breaks = c (-0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, Inf), 
                                   labels = c ("0 - 0.09", "0.10 - 0.19", "0.20 - 0.29", "0.30 - 0.39",
                                               "0.40 - 0.49", "0.50 - 0.59", "0.60 - 0.69", "0.70 - 0.79",
                                               "0.80 - 0.89", "0.90 - 1.00"))
ggplot (data = data.awa.cut) +
        geom_sf (data = data.awa.cut, aes (fill = brks_cut_2018)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Cutblock Density (km2/km2)") 
```

```{r, map sim cutblock density 2048, fig.cap = "Figure 13. Map of simulated cutblock density in 2048 by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cut <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.cut$brks_cut_2048 <- cut (data.awa.cut$p_030_204, 
                                   breaks = c (-0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, Inf), 
                                   labels = c ("0 - 0.09", "0.10 - 0.19", "0.20 - 0.29", "0.30 - 0.39",
                                               "0.40 - 0.49", "0.50 - 0.59", "0.60 - 0.69", "0.70 - 0.79",
                                               "0.80 - 0.89", "0.90 - 1.00"))
ggplot (data = data.awa.cut) +
        geom_sf (data = data.awa.cut, aes (fill = brks_cut_2048)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Cutblock Density (km2/km2)") 
```

```{r, map sim cutblock density 2078, fig.cap = "Figure 14. Map of simulated cutblock density in 2078 by freshwater assessment watershed area in the Okanagan Timber Supply Area.", eval = T, echo = F, message = F, warning = F}
data.awa.cut <- sf::st_read (dsn = here ("wildlife/data/okanagan/data_awa_sim_cut0to30_road.shp"),
                              quiet = T)
data.awa.cut$brks_cut_2078 <- cut (data.awa.cut$p_030_207, 
                                   breaks = c (-0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, Inf), 
                                   labels = c ("0 - 0.09", "0.10 - 0.19", "0.20 - 0.29", "0.30 - 0.39",
                                               "0.40 - 0.49", "0.50 - 0.59", "0.60 - 0.69", "0.70 - 0.79",
                                               "0.80 - 0.89", "0.90 - 1.00"))
ggplot (data = data.awa.cut) +
        geom_sf (data = data.awa.cut, aes (fill = brks_cut_2078)) +
        scale_fill_discrete () +
        scale_fill_brewer (type = "div",
                           palette = "RdBu",
                           direction = -1,
                           aesthetics = "fill",
                           na.value = "orange") +
        annotation_scale (location = "br", 
                          width_hint = 0.5, 
                          height = unit (0.1, "cm")) +
        annotation_north_arrow (location = "tl", 
                                which_north = "true", 
                                style = north_arrow_fancy_orienteering,
                                height = unit (1.25, "cm"),
                                width = unit (1, "cm")) +
        labs (fill = "Cutblock Density (km2/km2)") 
```

## 7. Conclusions
In general, the coarse-scale moose population and habitat indicators used here suggest that moose populations are at least stable across the majority of the TSA. There are no clear indications that previous forestry and other land use activites are negatively impacting the sustainability of moose populations in the region. Interestingly, road densities are high across much of the TSA, but they do not appear to be correlated with declining moose populations or high hunting pressure. In addition, there are high densities of cutblocks in the west-central portion of the TSA, but moose population indicators do not suggest a declining population there.       

Simulated future forestry in the base case suggests that road density may increase in portions of the TSA. However, increases were relatively small (10%), as road densities were already high, indicating that future roads may not be a particularly large concern for moose management. Simulated future cutblock densities were on average lower across the TSA, and distribution was less dispersed. it is possible that future lower cutblock densities may create a shortage of forage for moose in some areas.

Overall, the indicators suggest that at a coarse-scale, previous and future forestry activity has not, and potentially will not, have a clear negative impact on moose.