RAS tifs RAS2FIM Illinois

The FIM4 Algorithm (via Codebase Decomposition)

```{{r, setup, warning = FALSE, echo = FALSE, results=‘hide’}} library(magrittr) library(rgl, quietly=T) knitr::knit_hooks$set(webgl = hook_webgl)

unit_to_map <- ‘01050002’

path_to_output <- file.path(“C:/Users/stick/root/Dropbox/root/projects/floodmapping/tools/FrankenFIM/FrankenFIM/vignettes/vis”,fsep=.Platform\(file.sep) path_to_tmp_output <- file.path(path_to_output,"tmp",fsep = .Platform\)file.sep) path_to_maine_stock <- file.path(“C:/Users/stick/root//Dropbox/root/projects/floodmapping/methods/OWP_FIM/outputs/FIM4_maine/”,fsep = .Platform$file.sep)

path_to_fim4_unit_outputs <- path_to_maine_stock

all_hucs <- sf::st_read(file.path(“C:/Users/stick/root//Dropbox/root/database/hosted/water/HUC8.fgb”,fsep = .Platform$file.sep))

wbd8_clp_inputs <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“wbd8_clp.gpkg”,fsep=.Platform\(file.sep)) levelpaths <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_subset_streams_levelPaths.gpkg",fsep = .Platform\)file.sep)) nearby_hucs <- sf::st_transform(all_hucs,sf::st_crs(wbd8_clp_inputs))[wbd8_clp_inputs %>% sf::st_buffer(1000),]

test_pkg <- file.path(“~/data/temp/test.gpkg”,fsep = .Platform$file.sep) hf_flowpaths <- sf::read_sf(test_pkg, “flowpaths”) hf_divides <- sf::read_sf(test_pkg, “divides”) hf_flowlines <- sf::read_sf(test_pkg, “flowlines”) hf_network <- sf::read_sf(test_pkg, “network”)

huc_flowpaths <- hf_flowpaths[wbd8_clp_inputs[4,],] huc_divides <- hf_divides[hf_divides\(id %in% huc_flowpaths\)divide_id,] huc_network <- hf_network[hf_network\(divide_id %in% huc_flowpaths\)divide_id,]

huc_bbox <- sf::st_bbox(sf::st_transform(huc_divides,sf::st_crs(‘EPSG:4326’))) %>% unlist() %>% unname()

hf_bounds <- sf::st_bbox(sf::st_transform(huc_divides,sf::st_crs(‘EPSG:4326’))) %>% unlist() %>% unname()

branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branch_polygons.gpkg”,fsep = .Platform$file.sep)) branch_polys_selction <- branch_polys[25,] #2,11, 19 aoi_bbox <- sf::st_bbox(sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’))) %>% unlist() %>% unname()


## Reminder: Last time {data-background-image="Water_bw.png"}

* Flowpath: The HY_Features (OGC-16-032r2) are defined such that we use the term to describe "whatever exists as part of a network of waterbody and surface depressions and surfance channels" and we represent as a linear centerline
* Divide: Divides are the polygon partitoning of the landscape
* Flowline: A flowline is a one-dimensional (linear) feature that represents a flowing body of water and is functional similar to a flowpath but does not realize the catchment concept.

![](./vis/flowline_flowpath.png){fig-align="center" height="80%"}

> From: Hydrologic Modeling and River Corridor Applications of HY_Features Concepts. https://docs.ogc.org/per/22-040.html. and Figure 9 reproduced.

## Hydrofabric for VPU 1

```{r, VPU, warning = FALSE, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("VPU1_HF.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  hf_bounds <- sf::st_bbox(sf::st_transform(hf_divides,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()
  hf_divide_values <- unique(hf_divides$divide_id)
  # hf_divide_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(hf_divide_values)),hf_divide_values,na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    leafem::addFeatures(data = sf::st_transform(hf_divides,sf::st_crs('EPSG:4326')),stroke = TRUE, opacity = 0.8,fillOpacity = 0,fill = FALSE,weight = 1) %>%
    leafem::addFeatures(data = sf::st_transform(hf_flowpaths,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="blue") %>%
    # leaflet::addLabelOnlyMarkers(data=sf::st_transform(forecast_basins %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
    #                              label = ~`ID`,
    #                              labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
    #                              group="Index Labels") %>%
    leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("VPU1_HF.png"),fsep = .Platform$file.sep))
}

Focus in on HUC:01050002

```{{r, tVPU, warning = FALSE, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“VPU1_HF_HUC.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% leafem::addFeatures(data = sf::st_transform(hf_flowpaths,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“blue”) %>% leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs[4,],sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 0,weight = 0.2,color=“black”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leaflet::addLabelOnlyMarkers(data=sf::st_transform(nearby_hucs %>% sf::st_centroid(quiet = TRUE),sf::st_crs(‘EPSG:4326’)), label = ~name, labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T), group=“Index Labels”) %>% leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4]) mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“VPU1_HF_HUC.png”),fsep = .Platform$file.sep)) }


![](./vis/VPU1_HF_HUC.png)

## Focus in on HUC:01050002

```{r, HUC, warning = FALSE, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  huc_flowpaths <- hf_flowpaths[wbd8_clp_inputs[4,],]
  huc_divides <- hf_divides[hf_divides$id %in% huc_flowpaths$divide_id,]
  huc_network <- hf_network[hf_network$divide_id %in% huc_flowpaths$divide_id,]
  
  hf_bounds <- sf::st_bbox(sf::st_transform(huc_divides,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()
  hf_divide_values <- unique(huc_divides$divide_id)
  hf_divide_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(hf_divide_values)),hf_divide_values,na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    leafem::addFeatures(data = sf::st_transform(huc_divides,sf::st_crs('EPSG:4326')),stroke = TRUE, opacity = 0.8,fillOpacity = 0,fill = FALSE,weight = 1) %>%
    leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs[4,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 0,weight = 0.2,color="black") %>%
    leafem::addFeatures(data = sf::st_transform(huc_flowpaths,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="blue") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leaflet::addLabelOnlyMarkers(data=sf::st_transform(nearby_hucs %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
                                 label = ~`name`,
                                 labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
                                 group="Index Labels") %>%
    leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002.png"),fsep = .Platform$file.sep))
}

And By Divide

```{{r,divide, warning = FALSE, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_divide.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

huc_flowpaths <- hf_flowpaths[wbd8_clp_inputs[4,],] huc_divides <- hf_divides[hf_divides\(id %in% huc_flowpaths\)divide_id,] huc_network <- hf_network[hf_network\(divide_id %in% huc_flowpaths\)divide_id,]

hf_bounds <- sf::st_bbox(sf::st_transform(huc_divides,sf::st_crs(‘EPSG:4326’))) %>% unlist() %>% unname() hf_divide_values <- unique(huc_divides\(divide_id) hf_divide_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(hf_divide_values)),hf_divide_values,na.color = "transparent") map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>% leafem::addFeatures(data = sf::st_transform(huc_divides,sf::st_crs('EPSG:4326')),stroke = TRUE, opacity = 0.8,fillOpacity = 1,fill = TRUE,weight = 1,fillColor=hf_divide_pal(huc_divides\)divide_id)) %>% leafem::addFeatures(data = sf::st_transform(huc_flowpaths,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“blue”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leaflet::addLabelOnlyMarkers(data=sf::st_transform(nearby_hucs %>% sf::st_centroid(quiet = TRUE),sf::st_crs(‘EPSG:4326’)), label = ~name, labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T), group=“Index Labels”) %>% leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4]) mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_divide.png”),fsep = .Platform$file.sep)) }


![](./vis/HUC01050002_divide.png)

## What does this network look like?

```{r, visNetwork1, echo = FALSE}
huc_flowpaths <- hf_flowpaths[wbd8_clp_inputs[4,],]
huc_divides <- hf_divides[hf_divides$id %in% huc_flowpaths$divide_id,]
huc_network <- hf_network[hf_network$divide_id %in% huc_flowpaths$divide_id,]
g = dplyr::select(huc_network, id, toid) %>%
  dplyr::distinct() %>%
  igraph::graph_from_data_frame(directed = TRUE) 
visNetwork::visIgraph(g) %>%
  visNetwork::visOptions(highlightNearest = list(enabled = T, hover = T), nodesIdSelection = T)

FIM Pipeline

FIM Pipeline

Stepping through the pipeline

The FIM Pipeline: Derive Levelpaths

```{{r, lelvepaths, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_levels.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branch_polygons.gpkg”,fsep = .Platform$file.sep)) map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs(‘EPSG:4326’)),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color=“black”) %>% leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“#000000”) %>% leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys %>% sf::st_centroid(quiet = TRUE),sf::st_crs(‘EPSG:4326’)), label = ~levpa_id, labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T), group=“Index Labels”) %>% leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_levels.png”),fsep = .Platform$file.sep)) }


![](./vis/HUC01050002_levels.png)

## Aside: What does SCHISM see?

```{r, bs, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_schism.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  basedir <- "~/data/raw/nohrsc/owp_files/nwm/nwm_parameters/NWM_coastal_parameters.tar/coastal_parameters/coastal/atlgulf/"
  hgrid.gr3_file <- file.path(basedir,"hgrid.gr3")
  hgrid.gr3 <- data.table::fread(hgrid.gr3_file,skip = 2,nrows = 10537611-2)
  hgrid.gr3_sf <- sfheaders::sf_point(
    obj = hgrid.gr3
    , x = "V2"
    , y = "V3"
    , z = "V4"
    , keep = TRUE
  ) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% 
    sf::st_transform(sf::st_crs("EPSG:4326"))
  schism_pts <- hgrid.gr3_sf[sf::st_transform(wbd8_clp_inputs,sf::st_crs(hgrid.gr3_sf)) %>% sf::st_make_valid(),]
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.1,color="black") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leaflet::addLabelOnlyMarkers(data=sf::st_transform(nearby_hucs %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
                                 label = ~`name`,
                                 labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
                                 group="Index Labels") %>%
    leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_schism.png"),fsep = .Platform$file.sep))
}

The FIM Pipeline: Target branch

```{{r, tlevelpath, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_levels_select.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branch_polygons.gpkg”,fsep = .Platform$file.sep)) branch_polys_selction <- branch_polys[25,] #2,11, 19 map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs(‘EPSG:4326’)),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color=“black”) %>% leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“yellow”) %>% leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys %>% sf::st_centroid(quiet = TRUE),sf::st_crs(‘EPSG:4326’)), label = ~levpa_id, labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T), group=“Index Labels”) %>% leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_levels_select.png”),fsep = .Platform$file.sep)) }


![](./vis/HUC01050002_levels_select.png)

## The FIM Pipeline: Target branch callout

```{r, slevelpath, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_levels_select_lone.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"branch_polygons.gpkg",fsep = .Platform$file.sep))
  branch_polys_selction <- branch_polys[25,] #2,11, 19
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color="black") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    # leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color="yellow") %>%
    # leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
    #                              label = ~`levpa_id`,
    #                              labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
    #                              group="Index Labels") %>%
    leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_levels_select_lone.png"),fsep = .Platform$file.sep))
}

The FIM Pipeline: Target branch

```{{r, tbranch, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branch_polygons.gpkg”,fsep = .Platform\(file.sep)) branch_polys_selction <- branch_polys[25,] #2,11, 19 aoi_bbox <- sf::st_bbox(sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname() map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>% # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color="black") %>% leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths\)levpa_id==branch_polys_selction$levpa_id,],sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% # leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“yellow”) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 0.2,weight = 1,color=“#000000”) %>% leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys_selction %>% sf::st_centroid(quiet = TRUE),sf::st_crs(‘EPSG:4326’)), label = ~levpa_id, labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T), group=“Index Labels”) %>% # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4]) leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])

mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_selectlevel.png”),fsep = .Platform$file.sep)) }


![](./vis/HUC01050002_selectlevel.png)

## The FIM Pipeline: DEM

```{r,dem, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_selectlevel_DEM.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("dem_meters_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))

  gg_df <- as.data.frame(dem_meters,xy = T)
  gglegend <- ggplot2::ggplot() +
    ggplot2::geom_raster(data = gg_df, ggplot2::aes(x = x, y = y, fill = dem_meters_2678000065)) +
    ggplot2::scale_fill_stepsn(name = "3DEP 10m (m)", 
                               colors =c("forestgreen","yellow","tan","brown"),
                               breaks = round(seq(plyr::round_any(min(na.omit(gg_df$dem_meters_2678000065)), 10, f = floor),
                                            plyr::round_any(max(na.omit(gg_df$dem_meters_2678000065)), 10, f = ceiling), length = 7), 1),
                               values = scales::rescale(
                                 seq(plyr::round_any(min(na.omit(gg_df$dem_meters_2678000065)), 10, f = floor),
                                     plyr::round_any(max(na.omit(gg_df$dem_meters_2678000065)), 10, f = ceiling), length = 6)))
  llegend <- cowplot::get_legend(gglegend)
  cowplot::save_plot(plot=cowplot::plot_grid(llegend),filename=file.path(path_to_tmp_output,glue::glue("test_legend.png"),fsep = .Platform$file.sep))
  pal <- leaflet::colorNumeric(c("forestgreen","yellow","tan","brown"), terra::values(dem_meters),na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths$levpa_id==branch_polys_selction$levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    # leafem::addGeoRaster(dem_meters,colorOptions = leafem:::colorOptions(palette = pal,
    #                                                                      # breaks = as.numeric(c(0:xx$max_ffreq)),
    #                                                                      domain = c(dem_meters_values_min,dem_meters_values_max))) %>%
    # leafem::addGeoRaster(dem_meters,color= pal) %>%
    leaflet::addRasterImage(dem_meters, colors = pal, opacity = 0.9,maxBytes=42*1024*1024) %>%
    # leaflet::addLegend("bottomleft",title = "3DEP 10m (m)", pal = pal, values = terra::values(dem_meters)) %>% 
    leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
  map <- magick::image_read(file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
  legend <- magick::image_read(file.path(path_to_tmp_output,glue::glue("test_legend.png"),fsep = .Platform$file.sep)) %>% 
    magick::image_trim() %>% 
    magick::image_fill(
      color = "transparent", 
      refcolor = "white", 
      fuzz = 4,
      point = "+1+1"
    ) %>% 
    magick::image_scale("200")
  magick::image_write(image = magick::image_mosaic(c(map, legend)),path=file.path(path_to_output,glue::glue("HUC01050002_selectlevel_DEM.png"),fsep = .Platform$file.sep))
}

The FIM Pipeline: Flow Girds (boolean)

```{{r, fgrids,echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_flows_grid_boolean.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

flows_grid_boolean <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,“branches”,branch_polys_selction\(levpa_id,glue::glue("flows_grid_boolean_{branch_polys_selction\)levpa_id}.tif”),fsep = .Platform\(file.sep)) %>% terra::classify(cbind(-Inf, 0, NA), right=TRUE) pal <- leaflet::colorFactor(c("black","blue"), c(0,1),na.color = "transparent") map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>% # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>% # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>% leaflet::addRasterImage(flows_grid_boolean, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>% leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4]) # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)") map1 mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep), border_size = 4,just = "center") ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_flows_grid_boolean.png"),fsep = .Platform\)file.sep)) }


![](./vis/HUC01050002_selectlevel_flows_grid_boolean.png)

## The FIM Pipeline: Headwaters

```{r, heads,echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_selectlevel_headwaters.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  headwaters <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_headwater_points_subset.gpkg",fsep = .Platform$file.sep))
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
    leafem::addFeatures(data = sf::st_transform(headwaters,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,radius = 2,color="#FF0000") %>%
    # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
    leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
  # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
  map1
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_headwaters.png"),fsep = .Platform$file.sep))
}

The FIM Pipeline: Flow Direction

```{{r,fd, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_flowdir.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

# flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,“flowdir_d8_burned_filled.tif”,fsep = .Platform\(file.sep)) flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction\)levpa_id,glue::glue(“flowdir_d8_burned_filled_{branch_polys_selction\(levpa_id}.tif"),fsep = .Platform\)file.sep)) pal <- leaflet::colorFactor(c(”#862627”,“#26207a”,“#6563b7”,“#6563b7”,“#cbc6c0”,“#fee08a”,“#f4a429”,“#c2631f”),c(1,2,3,4,5,6,7,8),na.color = “transparent”) map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths\(levpa_id==branch_polys_selction\)levpa_id,],sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leaflet::addRasterImage(flowdir_d8_burned_filled, colors = pal, opacity = 1,maxBytes=4210241024) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color=“black”) %>% # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4]) leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4]) # addLegend(pal = pal, values = terra::values(dem_meters), title = “Elevation data for Santander (mts)”) map1 mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_flowdir.png”),fsep = .Platform$file.sep))

}


![](./vis/HUC01050002_selectlevel_flowdir.png)

## The FIM Pipeline: Burned & Filled Flow Direction

```{r, bf,echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_selectlevel_BFflowdir.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  # flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"flowdir_d8_burned_filled.tif",fsep = .Platform$file.sep))
  flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("flowdir_d8_burned_filled_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
  pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths$levpa_id==branch_polys_selction$levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leaflet::addRasterImage(flowdir_d8_burned_filled, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
    leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
    # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
    leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
  # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
  map1
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_BFflowdir.png"),fsep = .Platform$file.sep))

}

The FIM Pipeline: Flow Accumulation

```{{r, fa, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_faBFflowdir.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

flowaccum_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,“branches”,“0”,“flowaccum_d8_burned_filled_0.tif”,fsep = .Platform\(file.sep)) pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent") map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>% # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>% leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths\)levpa_id==branch_polys_selction\(levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>% leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>% leaflet::addRasterImage(flowaccum_d8_burned_filled, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>% leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4]) # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)") map1 mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep), border_size = 4,just = "center") ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_faBFflowdir.png"),fsep = .Platform\)file.sep))

}


![](./vis/HUC01050002_selectlevel_faBFflowdir.png)

## The FIM Pipeline: Slope

```{r,slope, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_selectlevel_Slope.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

  # slopes_d8_dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","slopes_d8_dem_meters_masked_0.tif",fsep = .Platform$file.sep))
  slopes_d8_dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("slopes_d8_dem_meters_masked_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
  pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths$levpa_id==branch_polys_selction$levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
    leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leaflet::addRasterImage(slopes_d8_dem_meters, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
    leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
  # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
  map1
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_Slope.png"),fsep = .Platform$file.sep))

}

test

The FIM Pipeline: DEM Derived Reaches

```{{r, demreach, echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_snreach.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

# sn_demDerived_reaches <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branches”,“0”,“demDerived_reaches_0.shp”,fsep = .Platform\(file.sep)) # sn_catchemnt_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","sn_catchments_reaches_0.tif",fsep = .Platform\)file.sep)) sn_demDerived_reaches <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,“branches”,branch_polys_selction\(levpa_id,glue::glue("demDerived_reaches_{branch_polys_selction\)levpa_id}.shp”),fsep = .Platform\(file.sep)) sn_catchemnt_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction\)levpa_id,glue::glue(“sn_catchments_reaches_{branch_polys_selction\(levpa_id}.tif"),fsep = .Platform\)file.sep))

catch_values <- unique(terra::values(sn_catchemnt_reaches)) pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(catch_values)),catch_values,na.color = “transparent”) map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% leafem::addFeatures(data = sf::st_transform(sn_demDerived_reaches,sf::st_crs(‘EPSG:4326’)),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color=“#0327cd”) %>% # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leaflet::addRasterImage(sn_catchemnt_reaches, colors = pal, opacity = 1,maxBytes=4210241024) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color=“black”) %>% # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4]) leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4]) # addLegend(pal = pal, values = terra::values(dem_meters), title = “Elevation data for Santander (mts)”) mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_snreach.png”),fsep = .Platform$file.sep))

}


![](./vis/HUC01050002_selectlevel_snreach.png)

## The FIM Pipeline: Gage Weighted (gw) Reaches

```{r, gwreach, echo = FALSE}
if(!file.exists(file.path(path_to_output,glue::glue("HUC01050002_selectlevel_gwreach.png"),fsep = .Platform$file.sep))) {
  unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
  
  # gw_catchments_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","gw_catchments_reaches_0.tif",fsep = .Platform$file.sep))
  gw_catchments_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("gw_catchments_reaches_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
  
  catch_values <- unique(terra::values(gw_catchments_reaches))
  pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(catch_values)),catch_values,na.color = "transparent")
  map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
    leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
    # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
    leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
    leaflet::addRasterImage(gw_catchments_reaches, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
    leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
    # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
    leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
  # addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
  mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = 1080,vheight = 1080)
  magick::image_write(
    magick::image_read(
      cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
    ),file.path(path_to_output,glue::glue("HUC01050002_selectlevel_gwreach.png"),fsep = .Platform$file.sep))

}

The FIM Pipeline: HAND

```{{r, HAND,echo = FALSE}} if(!file.exists(file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_HAND.png”),fsep = .Platform\(file.sep))) { unlink(file.path(path_to_tmp_output,"*",fsep = .Platform\)file.sep))

# HAND <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,“branches”,“0”,“rem_zeroed_masked_0.tif”,fsep = .Platform\(file.sep)) HAND <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction\)levpa_id,glue::glue(“rem_zeroed_masked_{branch_polys_selction\(levpa_id}.tif"),fsep = .Platform\)file.sep))

HAND_values <- terra::values(HAND) pal <- leaflet::colorNumeric(RColorBrewer::brewer.pal(5, “Blues”), c(min(HAND_values,na.rm = TRUE),max(HAND_values,na.rm = TRUE)),na.color = “transparent”) map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>% leaflet::addProviderTiles(“CartoDB.Positron”,group = “CartoDB.Positron”) %>% # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color=“#000000”) %>% leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color=“#000000”) %>% leaflet::addRasterImage(HAND, colors = pal, opacity = 1,maxBytes=5010241024) %>% leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs(‘EPSG:4326’)),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color=“black”) %>% # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4]) leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4]) # addLegend(pal = pal, values = terra::values(dem_meters), title = “Elevation data for Santander (mts)”) mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue(“test.png”),fsep = .Platform\(file.sep),vwidth = 1080,vheight = 1080) magick::image_write( magick::image_read( cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform\)file.sep), border_size = 4,just = “center”) ),file.path(path_to_tmp_output,glue::glue(“test_cropped.png”),fsep = .Platform$file.sep))

gg_df <- as.data.frame(HAND,xy = T) gglegend <- ggplot2::ggplot() + ggplot2::geom_raster(data = gg_df, ggplot2::aes(x = x, y = y, fill = rem_zeroed_masked_2678000065)) + ggplot2::scale_fill_stepsn(name = “HAND (m)”, colors =RColorBrewer::brewer.pal(5, “Blues”), breaks = round(seq(plyr::round_any(min(na.omit(gg_df\(rem_zeroed_masked_2678000065)), 10, f = floor), plyr::round_any(max(na.omit(gg_df\)rem_zeroed_masked_2678000065)), 10, f = ceiling), length = 7), 1), values = scales::rescale( seq(plyr::round_any(min(na.omit(gg_df\(rem_zeroed_masked_2678000065)), 10, f = floor), plyr::round_any(max(na.omit(gg_df\)rem_zeroed_masked_2678000065)), 10, f = ceiling), length = 6))) llegend <- cowplot::get_legend(gglegend) cowplot::save_plot(plot=cowplot::plot_grid(llegend),filename=file.path(path_to_tmp_output,glue::glue(“test_legend.png”),fsep = .Platform$file.sep))

map <- magick::image_read(file.path(path_to_tmp_output,glue::glue(“test_cropped.png”),fsep = .Platform\(file.sep)) legend <- magick::image_read(file.path(path_to_tmp_output,glue::glue("test_legend.png"),fsep = .Platform\)file.sep)) %>% magick::image_trim() %>% magick::image_fill( color = “transparent”, refcolor = “white”, fuzz = 4, point = “+1+1” ) %>% magick::image_scale(“200”) magick::image_write(image = magick::image_mosaic(c(map, legend)),path=file.path(path_to_output,glue::glue(“HUC01050002_selectlevel_HAND.png”),fsep = .Platform$file.sep))

}

![](./vis/HUC01050002_selectlevel_HAND.png)

## These are all interactive too!

```{r, HANDleaf,echo = FALSE}
HAND <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("rem_zeroed_masked_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
HAND_values <- terra::values(HAND)
pal <- leaflet::colorNumeric(RColorBrewer::brewer.pal(5, "Blues"), c(min(HAND_values,na.rm = TRUE),max(HAND_values,na.rm = TRUE)),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths$levpa_id==branch_polys_selction$levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(HAND, colors = pal, opacity = 1,maxBytes=50*1024*1024) %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
map1

Images for squirrels

{r, gif,echo = FALSE} if(!file.exists(file.path(path_to_output,glue::glue("FIM4_gif"),fsep = .Platform$file.sep))) { m <- magick::image_animate(c(magick::image_read(file.path(path_to_output,"VPU1_HF_HUC.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_divide.png",fsep = .Platform$file.sep)), # magick::image_read(file.path(path_to_output,"HUC01050002_schism.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_levels.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_levels_select.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_levels_select_lone.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_DEM.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_flows_grid_boolean.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_headwaters.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_flowdir.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_BFflowdir.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_faBFflowdir.png",fsep = .Platform$file.sep)), # magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_Slope.png",fsep = .Platform$file.sep)), # magick::image_read(file.path(path_to_output,HUC01050002_selectlevel_snreach.png,fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,"HUC01050002_selectlevel_gwreach.png",fsep = .Platform$file.sep)), magick::image_read(file.path(path_to_output,HUC01050002_selectlevel_HAND.png,fsep = .Platform$file.sep)) ), fps = 1, loop = 1, dispose = "previous") magick::image_write(m,file.path(path_to_output,glue::glue("FIM4_gif"),fsep = .Platform$file.sep)) }

Early improvements

  • Sane(er) base data handling
  • 10 fewer steps
  • 68 seconds off pre-processing and 55 seconds off a branch

Too long

```{r}
# ===========================
# FIM4
# ===========================
library(magrittr)
# library(leaflet)
# library(leafgl)
# remotes::install_github("doehm/cropcircles")

unit_to_map <- '01050002'

path_to_current <- file.path("~/data/output/FIM4_stock/",fsep = .Platform$file.sep)
path_to_new <- file.path("~/data/output/FIM4_refrepl/",fsep = .Platform$file.sep)
path_to_output <- file.path("~/data/output/vis",fsep = .Platform$file.sep)
path_to_tmp_output <- file.path(path_to_output,"tmp",fsep = .Platform$file.sep)
path_to_maine_stock <- file.path("~/Dropbox/root/projects/floodmapping/methods/OWP_FIM/outputs/FIM4_maine/",fsep = .Platform$file.sep)

path_to_fim4_unit_outputs <- path_to_maine_stock

all_hucs <- sf::st_read(file.path("~/Dropbox/root/database/hosted/water/HUC8.fgb",fsep = .Platform$file.sep))

wbd8_clp_inputs <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"wbd8_clp.gpkg",fsep=.Platform$file.sep))
levelpaths <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_subset_streams_levelPaths.gpkg",fsep = .Platform$file.sep))
nearby_hucs <- sf::st_transform(all_hucs,sf::st_crs(wbd8_clp_inputs))[wbd8_clp_inputs %>% sf::st_buffer(1000),]

hf_flowlines <- arrow::open_dataset("~/data/raw/lynker-spatial/hydrofabric/v2.2/final/conus/conus_flowpaths/") %>% dplyr::collect() %>% sf::st_as_sf()
sf::st_crs(hf_flowlines) <- sf::st_crs("EPSG:5070")
hf_divides <- arrow::open_dataset("~/data/raw/lynker-spatial/hydrofabric/v2.2/final/conus/conus_divides/") %>% dplyr::collect() %>% sf::st_as_sf()
sf::st_crs(hf_divides) <- sf::st_crs("EPSG:5070")

select_divides <- hf_divides[sf::st_transform(wbd8_clp_inputs,sf::st_crs(hf_divides)),]
select_flowlines <- hf_flowlines[hf_flowlines$divide_id %in% select_divides$divide_id,]

# Hydrofrabric maps (boundary conditions)
# ===========================
forecast_basins <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_catchments_proj_subset.gpkg",fsep = .Platform$file.sep))
forecast_basins_values <- unique(forecast_basins$ID)
forecast_basins_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(forecast_basins_values)),forecast_basins_values,na.color = "transparent")
forecast_network <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_subset_streams.gpkg",fsep = .Platform$file.sep))
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(forecast_basins,sf::st_crs('EPSG:4326')),stroke = FALSE, opacity = 0.8,fillOpacity = 1,fill = TRUE,weight = 1,fillColor=forecast_basins_pal(forecast_basins$ID)) %>%
  leafem::addFeatures(data = sf::st_transform(forecast_network,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="blue") %>%
  # leaflet::addLabelOnlyMarkers(data=sf::st_transform(forecast_basins %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
  #                              label = ~`ID`,
  #                              labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
  #                              group="Index Labels") %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# Hydrofabric
hydrofabric_divides_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(select_divides$divide_id)),select_divides$divide_id,na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(select_divides,sf::st_crs('EPSG:4326')),stroke = FALSE, opacity = 0.8,fillOpacity = 1,fill = TRUE,weight = 1,fillColor=hydrofabric_divides_pal(select_divides$divide_id)) %>%
  leafem::addFeatures(data = sf::st_transform(select_flowlines,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="blue") %>%
  # leaflet::addLabelOnlyMarkers(data=sf::st_transform(forecast_basins %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
  #                              label = ~`ID`,
  #                              labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
  #                              group="Index Labels") %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# Test Flowlines
hydrofabric_flowlpaths <- sf::st_read("~/data/temp/test.gpkg",layer = 'flowpaths')
hydrofabric_divides <- sf::st_read("~/data/temp/test.gpkg",layer = 'divides')
hydrofabric_flowlines <- sf::st_read("~/data/temp/test.gpkg",layer = 'flowlines')
hydrofabric_divides_pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(select_divides$divide_id)),select_divides$divide_id,na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(select_divides,sf::st_crs('EPSG:4326')),stroke = FALSE, opacity = 0.8,fillOpacity = 1,fill = TRUE,weight = 1,fillColor=hydrofabric_divides_pal(select_divides$divide_id)) %>%
  leafem::addFeatures(data = sf::st_transform(select_flowlines,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="blue") %>%
  # leaflet::addLabelOnlyMarkers(data=sf::st_transform(forecast_basins %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
  #                              label = ~`ID`,
  #                              labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
  #                              group="Index Labels") %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===============


# Coastal
basedir <- "~/data/raw/nohrsc/owp_files/nwm/nwm_parameters/NWM_coastal_parameters.tar/coastal_parameters/coastal/atlgulf/"
hgrid.gr3_file <- file.path(basedir,"hgrid.gr3")
hgrid.gr3 <- data.table::fread(hgrid.gr3_file,skip = 2,nrows = 10537611-2)
hgrid.gr3_sf <- sfheaders::sf_point(
  obj = hgrid.gr3
  , x = "V2"
  , y = "V3"
  , z = "V4"
  , keep = TRUE
) |> sf::st_set_crs(sf::st_crs("EPSG:6349")) %>% 
  sf::st_transform(sf::st_crs("EPSG:4326"))
schism_pts <- hgrid.gr3_sf[sf::st_transform(wbd8_clp_inputs,sf::st_crs(hgrid.gr3_sf)) %>% sf::st_make_valid(),]

# Branch 0
# ===========================
wbd_buffered_inputs <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"wbd_buffered.gpkg",fsep=.Platform$file.sep))
branch_0_bbox <- sf::st_bbox(sf::st_transform(wbd_buffered_inputs,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()

map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color="black") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addLabelOnlyMarkers(data=sf::st_transform(nearby_hucs %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
                               label = ~`name`,
                               labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
                               group="Index Labels") %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
branch_polys <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"branch_polygons.gpkg",fsep = .Platform$file.sep))
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color="black") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color="#000000") %>%
  leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
                               label = ~`levpa_id`,
                               labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
                               group="Index Labels") %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# BranchPolySelector===========================
branch_polys_selction <- branch_polys[25,] #2,11, 19
aoi_bbox <- sf::st_bbox(sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(schism_pts,sf::st_crs('EPSG:4326')),stroke = FALSE,opacity = 1,fillOpacity = 1,weight = 0.02,color="black") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  # leafem::addFeatures(data = sf::st_transform(branch_polys,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=TRUE,opacity = 1,fillOpacity = 0.2,weight = 1,color="yellow") %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 0.2,weight = 1,color="#000000") %>%
  leaflet::addLabelOnlyMarkers(data=sf::st_transform(branch_polys_selction %>% sf::st_centroid(quiet = TRUE),sf::st_crs('EPSG:4326')),
                               label = ~`levpa_id`,
                               labelOptions = leaflet::labelOptions(noHide = T, sticky = T,textOnly = T),
                               group="Index Labels") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# ===========================
dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"dem_meters.tif",fsep = .Platform$file.sep))
dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("dem_meters_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))

gg_df <- as.data.frame(dem_meters,xy = T)
gglegend <- ggplot2::ggplot() +
  ggplot2::geom_raster(data = gg_df, ggplot2::aes(x = x, y = y, fill = dem_meters_2678000065)) +
  ggplot2::scale_fill_stepsn(name = "3DEP 10m (m)", 
                             colors =c("forestgreen","yellow","tan","brown"),
                             breaks = seq(plyr::round_any(min(na.omit(gg_df$dem_meters_2678000065)), 10, f = floor),
                                          plyr::round_any(max(na.omit(gg_df$dem_meters_2678000065)), 10, f = ceiling), length = 7),
                             values = scales::rescale(
                               seq(plyr::round_any(min(na.omit(gg_df$dem_meters_2678000065)), 10, f = floor),
                                   plyr::round_any(max(na.omit(gg_df$dem_meters_2678000065)), 10, f = ceiling), length = 6))) 
llegend <- cowplot::get_legend(gglegend)

pal <- leaflet::colorNumeric(c("forestgreen","yellow","tan","brown"), terra::values(dem_meters),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths[levelpaths$levpa_id==branch_polys_selction$levpa_id,],sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  # leafem::addGeoRaster(dem_meters,colorOptions = leafem:::colorOptions(palette = pal,
  #                                                                      # breaks = as.numeric(c(0:xx$max_ffreq)),
  #                                                                      domain = c(dem_meters_values_min,dem_meters_values_max))) %>%
  # leafem::addGeoRaster(dem_meters,color= pal) %>%
  leaflet::addRasterImage(dem_meters, colors = pal, opacity = 0.9,maxBytes=42*1024*1024) %>%
  # leaflet::addLegend("bottomleft",title = "3DEP 10m (m)", pal = pal, values = terra::values(dem_meters)) %>% 
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep),vwidth = input$dimension[1])
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))

ggplot2::ggsave(cowplot::plot_grid(llegend, rel_widths = c(1)),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep),)
holepuch_map <- cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
magick::image_display(holepuch_map)

magick::image_read(holepuch_map)
magick::image_read(test)

magick::image_mosaic(c(bigdata, logo, frink))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# ===========================
flows_grid_boolean <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"flows_grid_boolean.tif",fsep = .Platform$file.sep))
flows_grid_boolean <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("flows_grid_boolean_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep)) %>%
  terra::classify(cbind(-Inf, 0, NA), right=TRUE)
pal <- leaflet::colorFactor(c("black","blue"), c(0,1),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(flows_grid_boolean, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
headwaters <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"nwm_headwater_points_subset.gpkg",fsep = .Platform$file.sep))
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  leafem::addFeatures(data = sf::st_transform(headwaters,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,radius = 2,color="#FF0000") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"flowdir_d8_burned_filled.tif",fsep = .Platform$file.sep))
flowdir_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("flowdir_d8_burned_filled_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(flowdir_d8_burned_filled, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
flowaccum_d8_burned_filled <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","flowaccum_d8_burned_filled_0.tif",fsep = .Platform$file.sep))
pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(flowaccum_d8_burned_filled, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
slopes_d8_dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","slopes_d8_dem_meters_masked_0.tif",fsep = .Platform$file.sep))
slopes_d8_dem_meters <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("slopes_d8_dem_meters_masked_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))
pal <- leaflet::colorFactor(c("#862627","#26207a","#6563b7","#6563b7","#cbc6c0","#fee08a","#f4a429","#c2631f"),c(1,2,3,4,5,6,7,8),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(flowaccum_d8_burned_filled, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# ===========================
sn_demDerived_reaches <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","demDerived_reaches_0.shp",fsep = .Platform$file.sep))
sn_catchemnt_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","sn_catchments_reaches_0.tif",fsep = .Platform$file.sep))
sn_demDerived_reaches <- sf::st_read(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("demDerived_reaches_{branch_polys_selction$levpa_id}.shp"),fsep = .Platform$file.sep))
sn_catchemnt_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("sn_catchments_reaches_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))

catch_values <- unique(terra::values(sn_catchemnt_reaches))
pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(catch_values)),catch_values,na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(data = sf::st_transform(sn_demDerived_reaches,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(sn_catchemnt_reaches, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

# ===========================
gw_catchments_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","gw_catchments_reaches_0.tif",fsep = .Platform$file.sep))
gw_catchments_reaches <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("gw_catchments_reaches_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))

catch_values <- unique(terra::values(gw_catchments_reaches))
pal <- leaflet::colorFactor(randomcoloR::distinctColorPalette(length(catch_values)),catch_values,na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(gw_catchments_reaches, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))
# ===========================
HAND <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches","0","rem_zeroed_masked_0.tif",fsep = .Platform$file.sep))
HAND <- terra::rast(file.path(path_to_fim4_unit_outputs,unit_to_map,"branches",branch_polys_selction$levpa_id,glue::glue("rem_zeroed_masked_{branch_polys_selction$levpa_id}.tif"),fsep = .Platform$file.sep))

HAND_values <- terra::values(HAND)
pal <- leaflet::colorNumeric(RColorBrewer::brewer.pal(5, "Blues"), c(min(HAND_values,na.rm = TRUE),max(HAND_values,na.rm = TRUE)),na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(wbd8_clp_inputs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1.2,color="#000000") %>%
  leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>%
  leaflet::addRasterImage(HAND, colors = pal, opacity = 1,maxBytes=50*1024*1024) %>%
  leafem::addFeatures(data = sf::st_transform(branch_polys_selction,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 2,color="black") %>%
  # leaflet::fitBounds(branch_0_bbox[3],branch_0_bbox[2],branch_0_bbox[1],branch_0_bbox[4])
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
map1
mapview::mapshot2(map1,file=file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(path_to_tmp_output,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(path_to_tmp_output,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(path_to_tmp_output,"*",fsep = .Platform$file.sep))

```