From 07ab0d4184311cdbbdd3b5da2d39ecfc03480dfd Mon Sep 17 00:00:00 2001
From: David Blodgett <dblodgett@usgs.gov>
Date: Tue, 16 Aug 2022 14:25:50 -0500
Subject: [PATCH] fixes #87 fixes #92
---
workspace/01_NHD_prep.Rmd | 4 +-
workspace/05_hyRefactor_flines.Rmd | 8 ++--
workspace/06-1_hyRefactor_cats.Rmd | 10 ++---
workspace/07-1_merge.Rmd | 27 ++++++++++--
workspace/R/utils.R | 68 ++++++++++++++++++++++++++++++
5 files changed, 103 insertions(+), 14 deletions(-)
diff --git a/workspace/01_NHD_prep.Rmd b/workspace/01_NHD_prep.Rmd
index d21e195..b48bcc6 100644
--- a/workspace/01_NHD_prep.Rmd
+++ b/workspace/01_NHD_prep.Rmd
@@ -44,7 +44,7 @@ fline <- sf::st_cast(fline, "LINESTRING")
catchment <- sf::read_sf(file.path(data_paths$nhdplus_dir, "reference_catchments.gpkg"))
-catchment <- sf::st_transform(catchment, 5070)
+catchment <- sf::st_transform(catchment, crs)
# we can remove truely degenerate COMIDs
# for 0 upstream area and no catchment area
@@ -97,7 +97,7 @@ for(VPU in vpu_codes) {
reachcode, frommeas, tomeas, pathlength, arbolatesu,
ftype, fcode, vpuid, rpuid, wbareacomi)
- write_sf(nhd, nav_gpkg, nhd_flowline)
+ write_sf(sf::st_transform(nhd, crs), nav_gpkg, nhd_flowline)
write_sf(cat_network, nav_gpkg, nhd_network)
cats_rpu <- full_cat_table %>%
diff --git a/workspace/05_hyRefactor_flines.Rmd b/workspace/05_hyRefactor_flines.Rmd
index 3e80585..d262403 100644
--- a/workspace/05_hyRefactor_flines.Rmd
+++ b/workspace/05_hyRefactor_flines.Rmd
@@ -118,8 +118,8 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) {
exclude_cats = c(outlets$COMID, avoid$COMID, POI_downstream$COMID),
warn = TRUE)
- write_sf(read_sf(tf), out_refac_gpkg, refactored_layer)
- write_sf(read_sf(tr), out_refac_gpkg, reconciled_layer)
+ write_sf(st_transform(read_sf(tf), crs), out_refac_gpkg, refactored_layer)
+ write_sf(st_transform(read_sf(tr), crs), out_refac_gpkg, reconciled_layer)
unlink(list(tf, tr))
rm(tf, tr)
@@ -128,7 +128,8 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) {
```{r lookup}
# create lookup for ref flowlines to use in the non-dendritic steps
-refactor_lookup <- dplyr::select(st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)), ID, member_COMID) %>%
+refactor_lookup <- st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)) %>%
+ dplyr::select(ID, member_COMID) %>%
dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
tidyr::unnest(cols = member_COMID) %>%
dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
@@ -139,7 +140,6 @@ if(is.character(refactor_lookup$reconciled_ID)) refactor_lookup$reconciled_ID <-
lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID")
-
readr::write_csv(lookup_table, lookup_table_file)
write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
```
diff --git a/workspace/06-1_hyRefactor_cats.Rmd b/workspace/06-1_hyRefactor_cats.Rmd
index d55bde2..4fb8fc5 100644
--- a/workspace/06-1_hyRefactor_cats.Rmd
+++ b/workspace/06-1_hyRefactor_cats.Rmd
@@ -35,16 +35,16 @@ if(needs_layer(out_refac_gpkg, divide_layer)) {
fac <- data_paths$fac[[rpu_string]]
fdr_temp <- terra::rast(fdr)
- crs <- terra::crs(fdr_temp)
+ raster_crs <- terra::crs(fdr_temp)
# reconciled have no ID applied in refactoring flowlines
- reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
+ reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), raster_crs)
# refactored has the original COMIDs and other attributes
- refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), crs) %>%
+ refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), raster_crs) %>%
distinct(.keep_all = T)
- cats <- st_transform(read_sf(out_refac_gpkg, nhd_catchment), crs)
+ cats <- st_transform(read_sf(out_refac_gpkg, nhd_catchment), raster_crs)
divides <- reconcile_catchment_divides(catchment = cats,
fline_ref = refactored,
@@ -59,7 +59,7 @@ if(needs_layer(out_refac_gpkg, divide_layer)) {
if(exists("divides")) unlink(cache_split)
- write_sf(divides, out_refac_gpkg, divide_layer)
+ write_sf(sf::st_transform(divides, crs), out_refac_gpkg, divide_layer)
}
sf::st_layers(out_refac_gpkg)
diff --git a/workspace/07-1_merge.Rmd b/workspace/07-1_merge.Rmd
index f9f4cb6..ac6c51b 100644
--- a/workspace/07-1_merge.Rmd
+++ b/workspace/07-1_merge.Rmd
@@ -62,12 +62,33 @@ if(needs_layer(gf_gpkg, agg_cats_layer)) {
sf::write_sf(merged_layers[[mapped_outlets_layer]], gf_gpkg, mapped_outlets_layer)
+ merged_layers[[agg_cats_layer]] <- merged_layers[[agg_cats_layer]] %>%
+ mutate(areasqkm = as.numeric(units::set_units(sf::st_area(.), "km^2")))
+
sf::write_sf(merged_layers[[agg_cats_layer]], gf_gpkg, agg_cats_layer)
-
- sf::write_sf(merged_layers[[agg_fline_layer]], gf_gpkg, agg_fline_layer)
+
+ merged_layers[[agg_fline_layer]] <- merged_layers[[agg_fline_layer]] %>%
+ mutate(lengthkm = as.numeric(units::set_units(sf::st_length(.), "km"))) %>%
+ left_join(select(sf::st_drop_geometry(merged_layers[[agg_cats_layer]]), ID, areasqkm), by = "ID")
+
+ sf::write_sf(merged_layers[[agg_fline_layer]], gf_gpkg, agg_fline_layer)
- sf::write_sf(merged_layers[[lookup_table_refactor]],
+ sf::write_sf(merged_layers[[lookup_table_refactor]],
gf_gpkg, lookup_table_refactor)
+ st_drop_geometry(merged_layers[[reconciled_layer]]) %>%
+ select(id = ID, toID = toID, lengthkm = LENGTHKM, levelpathid = LevelPathID) %>%
+ left_join(select(st_drop_geometry(merged_layers[[divide_layer]]),
+ id = ID, areasqkm), by = "id") %>%
+ relocate(id, toID, lengthkm, areasqkm, levelpathid) %>%
+ write_sf(rfc_gpkg, catchment_network_table)
+
+ st_drop_geometry(merged_layers[[agg_fline_layer]]) %>%
+ select(id = ID, toid = toID, lengthkm, areasqkm) %>%
+ left_join(distinct(select(merged_layers[[lookup_table_refactor]],
+ id = aggregated_flowpath_ID,
+ levelpathid = LevelPathID))) %>%
+ sf::write_sf(gf_gpkg, catchment_network_table)
+
}
```
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index 84561ad..1d2cbfc 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -466,6 +466,7 @@ merge_refactor <- function(rpus, rpu_vpu_out,
out
}
+# FIXME: THESE FUNCTIONS SHOULD COME FROM hyFab
#' Extract nexus locations for Reference POIs
#' @param gpkg a reference hydrofabric gpkg
#' @param type the type of desired POIs
@@ -505,3 +506,70 @@ nexus_from_poi = function(mapped_POIs,
nexus_locations
}
+
+#' Generate Lookup table for Aggregated Nextwork
+#' @param gpkg to the aggregated network. The look up table will be added here as a layer called "lookup_table"
+#' @param flowpath_name the name of the flowpath layer in gpkg
+#' @param refactored_fabric the gpkg for the corresponding VPU reference fabric
+#' @return data.frame
+#' @export
+#' @importFrom sf read_sf st_drop_geometry write_sf
+#' @importFrom dplyr mutate select full_join left_join
+#' @importFrom tidyr unnest
+
+generate_lookup_table = function(nl,
+ mapped_POIs) {
+
+ # THIS IS THE CODE ABOVE!!!!
+ nexus_locations = nexus_from_poi(mapped_POIs, verbose = FALSE)
+
+ lu = nl$lookup_table %>%
+ select(reconciled_ID, NHDPlusV2_COMID)
+
+ names(nl$flowpaths) <- tolower(names(nl$flowpaths))
+
+ nl$flowpaths <- rename(nl$flowpaths, member_COMID = set)
+
+ lu2 = nl$flowpaths %>%
+ st_drop_geometry() %>%
+ select(
+ aggregated_ID = id,
+ toID = toid,
+ member_COMID = member_comid,
+ divide_ID = id,
+ POI_ID = poi_id,
+ mainstem_id = levelpathid
+ ) %>%
+ mutate(member_COMID = strsplit(member_COMID, ","),
+ POI_ID = as.integer(POI_ID)) %>%
+ unnest(col = 'member_COMID') %>%
+ mutate(NHDPlusV2_COMID_part = sapply( strsplit(member_COMID, "[.]"),
+ FUN = function(x){ x[2] }),
+ NHDPlusV2_COMID_part = ifelse(is.na(NHDPlusV2_COMID_part), 1L,
+ as.integer(NHDPlusV2_COMID_part)),
+ NHDPlusV2_COMID = sapply( strsplit(member_COMID, "[.]"),
+ FUN = function(x){ as.numeric(x[1]) })
+ ) %>%
+ full_join(lu,
+ by = "NHDPlusV2_COMID")
+
+
+ group_by(aggregated_ID) %>%
+ arrange(hydroseq) %>%
+ mutate(POI_ID = c(POI_ID[1], rep(NA, n()-1))) %>%
+ ungroup() %>%
+ select(-hydroseq, - member_COMID) %>%
+ left_join(select(nexus_locations, POI_ID = poi_id, value, type),
+ by = "POI_ID") %>%
+ select(NHDPlusV2_COMID, NHDPlusV2_COMID_part,
+ reconciled_ID,
+ aggregated_ID, toID, mainstem = mainstem_id,
+ POI_ID, POI_TYPE = type, POI_VALUE = value)
+
+
+ write_sf(lu2, gpkg, "lookup_table")
+
+ return(gpkg)
+
+}
+
--
GitLab