03b-SiteRegina.Rmd

```{r setupsiteregina, include=FALSE}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
library(raster)
library(leaflet)
library(sf)
library(lidR)
library(lwgeom)
library(osmdata)
theme_set(bayesplot::theme_default())
opts_chunk$set(
  echo = F, message = F, warning = F, fig.height = 6, fig.width = 8,
  cache = T, cache.lazy = F)
```

```{r datasiteregina, include = F}
crs <- '+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0' # global crs definition
dem <- raster("data/Paracou/lidar/DEM_1m_2015.tif") # for CRS
first_candidates <- st_read("data/Regina/trees/Regina TREEMUTATION candidates.gpx", layer = "waypoints")
angela <- st_read("data/Regina/trees/Angela.gpx")
las <- lapply(list.files("data/Regina/lidar/ONF/LAZ_zoi_trees/", pattern = ".laz$", full.names = T), readLAS)
dcm <- lapply(las, grid_canopy, res = 1, algorithm = p2r())
dcm <- do.call(merge, dcm)
projection(dcm) <- projection(dem)
dcm <- projectRaster(dcm, crs = crs)
# rm(dem) 
inventory <- lapply(c("AG_ONF.gpx", "EBR_ONF.gpx", "EBV_ONF.gpx", "MN_ONF.gpx"),
                function(file)
                  st_read(file.path("data", "Regina", "trees", file), layer = "waypoints")) %>% 
  bind_rows() %>% 
  separate(name, c("species", "id", "dbh"), remove = F, convert = T)
road <- opq(bbox = st_bbox(inventory)) %>%
  add_osm_feature(key = 'highway', value = "primary") %>%
    add_osm_feature(key = 'name', value = "Route de l'Est") %>%
    osmdata_sf()
road <- road$osm_lines %>% 
  st_union()
inventory$distroad <- as.numeric(st_distance(inventory, road))
new_candidates <- filter(inventory, (distroad < 100 & !grepl("E", species)) | 
                           (distroad < 200 & grepl("E", species))) %>% 
  filter((species == "AG" & dbh > 90) | 
           (species != "AG" & dbh > 60))
# new_candidates %>% 
#   dplyr::select(-species, -id, -dbh, -distroad) %>% 
#   st_write(dsn = "data/Regina/trees/Regina_TREEMUTATION_selected_candidates.gpx", 
#            delete_dsn = T, driver = "GPX", layer = "waypoints")
candidates <- bind_rows(separate(first_candidates, name, c("species", "id", "dbh"), remove = F, convert = T),
                        filter(new_candidates, id %in% c("5508", "5505")),
                        mutate(angela, name = "angela", id = NA, species = "AG"))
# candidates %>% 
#   dplyr::select(name, geometry) %>% 
#   st_write(dsn = "data/Regina/trees/Regina_TREEMUTATION_candidates2.gpx",
#            delete_dsn = T, driver = "GPX", layer = "waypoints")
candidates$species[3] <- "MN"
rm(first_candidates, new_candidates, angela, inventory)
drone <- lapply(list.files("data/Regina/drone/20210309/", full.names = T), read_csv) %>% 
  bind_rows() %>% 
  st_as_sf(coords = c("longitude", "latitude"),
                 crs = '+proj=latlong +zone=22 +datum=WGS84')
# trajectosONF <- vroom::vroom("data/Regina/lidar/Trajectos/Trajectos_241116.txt", 
#              col_names = c("TempsGPS", "ROLL", "Pict", "Yaw", "X", "Y", "H")) %>% 
#   st_as_sf(coords = c("X", "Y", "H"), crs = projection(dem)) %>% 
#   st_combine() %>% 
#   st_cast("POLYGON") %>% 
#   st_transform(crs = crs)
```

## Régina

Plus:

* 1 LiDAR in 2016 (ONF, see Fig. \@ref(fig:reginaMap)
* Possibility of clearing for drone
* Destructive sampling for tree rings
* Taller trees (older? and more cell divisions)

Cons:

* Distance for phenology

```{r reginaMap, fig.cap="Regina individuals and MNC from ONF's LiDAR."}
leaflet() %>%
  addTiles(group = "OSM") %>%
  # addPolygons(st_zm(trajectosONF), group = "trajectos", col = "blue") %>% 
  addRasterImage(dcm, colors = "Greys", group = "DCM") %>%
  addCircles(data = candidates, label = ~name, col = "red", radius = 20, group = "Trees") %>% 
  addCircles(data = drone, col = "blue", radius = 2, label = ~filename, group = "Drone") %>% 
  addLayersControl(overlayGroups = c("OSM", "trajectos", "DCM", "Trees", "Drone"),
                   options = layersControlOptions(collapsed = T))  
```