wsep.t2: WSEP Tier 2

test_that("readme works", {


  #---------------------------------
  # CREATE THE STREAM LAYER

  # Load the WSEP Tier 2 R-package
  library(wsep.t2)
  library(dplyr)

  # Load the sample data set for the
  # Tsolum River near Courtenay BC
  data(TsolumStreams)

  # Will must project our data so that
  # x,y are cartesian coordinates with units of meters
  strm <- utm_projection(data = TsolumStreams)


  # Then constrain the sampling frame to remove 1st
  # order tribs less than 600m in length
  # and clip the upper 300m off of other tributaries
  c_ctrm <- constrain_streams(strm = strm,
                          length_remove = 600,
                          length_trim = 300)



  # Remove lakes and other lotic reaches
  c_strm <- remove_lentic_bcfwa(strm = c_ctrm,
                                 EDGE_TYPE = "EDGE_TYPE")

  # Consider removing any other lakes manually or through a simple filter

  # Remove alpine areas - all segments over 800m
  ca_strm <- remove_alpine_bcfwa(strm = c_strm,
                              elevation_threshold = 800)


  # (Optional) Visualize original (raw) and constrained streams
  strm1_plot <- sf::st_zm(strm)
  strm2_plot <- sf::st_zm(ca_strm)
  plot(sf::st_geometry(strm1_plot), col = "lightgrey",
       main = "Adjustments")
  plot(sf::st_geometry(strm2_plot), col = "blue", add = TRUE)
  legend("topright", c("original", "adjusted"),
         col = c("lightgrey", "blue"), lwd = 1)

  # Finalize and adjust with any additional filters


  #---------------------------------
  # DEFINE STRATA
  # Create a new field called `strata` that divides the remaining
  # stream network based on stream order: `stratum_1` < 3rd order
  # streams and `stratum_2` ≥ 3rd order streams.

  ca_strm$strata <- NA
  ca_strm$strata <- ifelse(ca_strm$STREAM_ORDER < 3, "stratum_1", ca_strm$strata)
  ca_strm$strata <- ifelse(ca_strm$STREAM_ORDER >= 3, "stratum_2", ca_strm$strata)
  plot(ca_strm["strata"], main = "Sampling Stratum")


  #---------------------------------
  # STREAM CROSSINGS
  # A list of stream crossings in the watershed can
  # be generated by taking the intersection of the
  # stream layer and the road layer, providing a list
  # of all mapped stream crossings

  data("TsolumRoads")
  # Ensure roads match projection
  roads <- utm_projection(data = TsolumRoads)
  # Define crossings as the intersection
  crossings <- suppressWarnings({ sf::st_intersection(ca_strm, roads) })


  #--------------------------------------------
  # Site Type A (stream crossing):
  # 1.	Generate a random sample  from the list
  # of stream crossings for each of the two strata
  # (< 3rd order vs. ≥ 3rd order).
  # 2.	Create a field checklist with at least the
  # following fields: unique identifier (e.g.,
  # WatershedName_SD_A_001), coordinates, and Strata.

  site_type_a <- grouped_random_sample(
                    data = crossings,
                    group_name = "strata",
                    n = 20)




  # (Optional) visualize
  strm_plot <- sf::st_zm(ca_strm)
  road_plot <- sf::st_zm(roads)
  plot(sf::st_geometry(strm_plot), col = "darkblue", main = "Site Type A (stream crossing)")
  plot(sf::st_geometry(road_plot), add = TRUE, col = "burlywood")
  plot(sf::st_geometry(site_type_a), add = TRUE, col = ifelse(site_type_a$strata == "stratum_1", "black", "red"), pch = 19)
  legend("topright",
         c("roads", "streams", "stratum 1", "stratum 2"),
         col = c("burlywood", "darkblue", "black", "red"),
         lwd = c(1, 1, NA, NA),
         pch = c(NA, NA, 19, 19))

  #--------------------------------------------
  # Site Type B (road proximity):
  # 1.	Apply a buffer (20m for <3rd order and 40 m for ≥3rd order) streams.
  # 2.	Taking the intersection of roads and this buffer.
  # 3.	Removing any stream crossings (site Type A),
  # by excluding any cases within 100 m of a crossing to avoid double counting.
  # 4.	Removing any segments < 50 m in length unless they are near a
  # switch-back (determined by manual review of the map).
  # [manual step not automated in R code]
  # 5.	Providing the start point of the segment as well as the segment
  # length and strata id (<3rd order vs. ≥ 3rd order) and mapping the entire
  # segment on the field maps to facilitate sampling.
  # 6.	Generate a random sample from the Site Type B list the strata
  # with <3rd order streams.
  # 7.	Append the complete list of Site Type B from the ≥ 3rd order strata.
  # 8.	Create a field checklist with at least the following fields:
  # unique identifier (e.g., WatershedName_SD_B_001), coordinates of
  # start point and end point, segment length, and Strata.

  type_b <- road_proximity_sample(
    n = 20,
    strm = ca_strm,
    roads = roads,
    buffer_s1_m = 50,
    buffer_s2_m = 90,
    buffer_crossings_m = 100,
    small_strm_segment_m = 50,
    stream_order = "STREAM_ORDER"
  )

  site_type_b <- type_b$points


  # (Optional) visualize
  strm_plot <- sf::st_zm(ca_strm)
  road_plot <- sf::st_zm(roads)
  plot(sf::st_geometry(strm_plot), col = "darkblue", main = "Site Type B (road proximity)")
  plot(sf::st_geometry(road_plot), add = TRUE, col = "burlywood")
  plot(sf::st_geometry(site_type_b), add = TRUE, col = ifelse(site_type_a$strata == "stratum_1", "black", "red"), pch = 19)
  legend("topright",
         c("roads", "streams", "stratum 1", "stratum 2"),
         col = c("burlywood", "darkblue", "black", "red"),
         lwd = c(1, 1, NA, NA),
         pch = c(NA, NA, 19, 19))


  site_type_c <- strm_grts(n = 20, strm = ca_strm, stream_order = 'STREAM_ORDER')

   plot(sf::st_geometry(strm_plot), col = "darkblue", main = "Site Type C (riparian-crossings)")
  plot(sf::st_geometry(road_plot), add = TRUE, col = "burlywood")
  plot(sf::st_geometry(site_type_c), add = TRUE, col = ifelse(site_type_b$strata == "stratum_1", "black", "red"), pch = 19)
  legend("topright", c("roads", "streams", "stratum 1", "stratum 2"), col = c("burlywood", "darkblue", "black", "red"), lwd = c(1, 1, NA, NA), pch = c(NA, NA, 19, 19))


  # Run tests
  testthat::expect_true(nrow(site_type_a) > 0)
  testthat::expect_true(all(names(type_b) == c("points", "line_segments")))
  testthat::expect_true(nrow(site_type_c) > 0)



  if(FALSE) {

    # Output directory
    output_dir <- "C:/Users/mbayly/Desktop/delete/my_sites"

    export_sites(output_dir = output_dir,
                 site_type_a = site_type_a,
                 type_b = type_b,
                 site_type_c = site_type_c,
                 export_csv = TRUE,
                 export_shp = TRUE,
                 export_kml = TRUE)
  }



})