Leftover fragments from catchment reconciliation

[mikejohnson51]

Hey @dblodgett-usgs,

This is something that is being cleaned up in the ngen work but that I wanted to bring to your attention. Here I am using a GPKG containing the upstream network and catchments of gage 01013500 (CAMELS basin). Within reconcile_catchment_divides fragments of a split catchment are assigned to the disjoint split element. If this is an artifact you care about addressing in hyRefactor I can send what I have via PR (it is a function to run on the output).

It also addresses other geometry challenges we have discussed.

Happy to share the GPKG as well but this issue is intended to be more illustrative.

library(hyRefactor)
library(sf)
library(dplyr)

# Path to network subset
path   <- '/Volumes/Transcend/ngen/refactor-tests/CAMELS/base-runs/gage_01013500.gpkg'
#Path to FACFDR tif files 
facfdr <- '/Volumes/Transcend/ngen/refactor-tests/base-data/fdrfac'

flowpaths  <- read_sf(path, "flowpaths")
catchments <- read_sf(path, "catchments")

tf <- tempfile(pattern = "refactored", fileext = ".gpkg")
tr <- tempfile(pattern = "reconciled", fileext = ".gpkg")

refactor_nhdplus(nhdplus_flines              = flowpaths, 
                 split_flines_meters         = 10000, 
                 split_flines_cores          = 1, 
                 collapse_flines_meters      = 1000,
                 collapse_flines_main_meters = 1000,
                 out_refactored = tf, 
                 out_reconciled = tr, 
                 three_pass          = TRUE, 
                 purge_non_dendritic = FALSE, 
                 warn = FALSE)

rpus <-unique(flowpaths$RPUID)
rpus <-rpus[!is.na(rpus)]

fdrfac_files <- list.files(facfdr, pattern = rpus, full.names = TRUE)
fdr          <- raster::raster(grep("_fdr", fdrfac_files, value = TRUE))
fac          <- raster::raster(grep("_fac", fdrfac_files, value = TRUE))
catchments   <-  st_transform(catchments, st_crs(fdr)) 
st_precision(catchments) <- raster::res(fdr)[1]

reconciled <- st_transform(read_sf(tr),  st_crs(fdr)) 
refactored <- st_transform(read_sf(tf),  st_crs(fdr)) 

divides    <- reconcile_catchment_divides(catchment = catchments,
                                          fline_ref = refactored,
                                          fline_rec = reconciled,
                                          fdr       = fdr,
                                          fac       = fac,
                                          para      = 1, 
                                          cache     = NULL) 
#> Loading required namespace: rgeos

plot(filter(divides, ID == 489)$geom) 

plot(filter(catchments, FEATUREID == 724556)$geom)
plot(filter(divides, ID == 489)$geom, col = "red", border = "red", add = TRUE) 

^{Created on 2021-09-14 by the reprex package (v2.0.1)}

[dblodgett-usgs]

Oh -- yes please. Maybe just use your example here as an example / test? Feel free to package up the example / test artifacts in inst/extdata/ Thanks!!

[mikejohnson51]

Right on! Do you want it as an exported function, or one that is called at the end/within of reconcile_catchment_divides? Or both?

[dblodgett-usgs]

Up to you -- I would guess it won't hurt to just call it at the end?

[mikejohnson51]

I like that personally... hides the mischief :)