diff --git a/NAMESPACE b/NAMESPACE index db30bcb..5ea6c28 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -205,3 +205,5 @@ importFrom(units,drop_units) importFrom(units,set_units) importFrom(utils,capture.output) importFrom(utils,str) +importFrom(yyjsonr,read_geojson_file) +importFrom(yyjsonr,write_geojson_file) diff --git a/R/aggregate_along_mainstems.R b/R/aggregate_along_mainstems.R index d784009..871b09b 100644 --- a/R/aggregate_along_mainstems.R +++ b/R/aggregate_along_mainstems.R @@ -68,8 +68,7 @@ aggregate_along_mainstems = function(network_list, hl_un, ideal_size_sqkm, min_area_sqkm, - min_length_km - ) + min_length_km) ) %>% ungroup() %>% group_by(levelpathid, ind) %>% @@ -79,6 +78,13 @@ aggregate_along_mainstems = function(network_list, hydroseq, member_comid, poi_id, n) + index_table |> + group_by(id) |> + mutate(n = n()) |> + filter(n > 1) + + + v = aggregate_sets(network_list, index_table) v = add_network_type(v, verbose = verbose) @@ -349,71 +355,94 @@ aggregate_sets = function(network_list, index_table) { #### - single_flowpaths = filter(index_table, n == 1) %>% + single_flowpaths = tryCatch({ + filter(index_table, n == 1) %>% + inner_join(network_list$flowpaths, by = "id") %>% + st_as_sf() %>% + select(set) %>% + rename_geometry("geometry") + }, error = function(e){ + NULL + }) + + flowpaths_out = tryCatch({ + filter(index_table, n > 1) %>% inner_join(network_list$flowpaths, by = "id") %>% st_as_sf() %>% select(set) %>% + union_linestrings('set') %>% rename_geometry("geometry") + }, error = function(e){ + NULL + }) - flowpaths_out = filter(index_table, n > 1) %>% - inner_join(network_list$flowpaths, by = "id") %>% - st_as_sf() %>% - select(set) %>% - union_linestrings('set') %>% - rename_geometry("geometry") %>% + flowpaths_out = flowpaths_out |> bind_rows(single_flowpaths) %>% left_join(set_topo_fin, by = "set") %>% rename(id = set, toid = toset) #### + single_catchments = tryCatch({ + filter(index_table, n == 1) %>% + inner_join(network_list$catchments, by = "id") %>% + st_as_sf() %>% + select(set) %>% + rename_geometry("geometry") + }, error = function(e){ + NULL + }) + + catchments_out = tryCatch({ + filter(index_table, n != 1) %>% + inner_join(network_list$catchments, by = "id") %>% + st_as_sf() %>% + select(set) %>% + union_polygons('set') %>% + mutate(areasqkm = add_areasqkm(.)) + }, error = function(e){ + NULL + }) - single_catchments = filter(index_table, n == 1) %>% - inner_join(network_list$catchments, by = "id") %>% - st_as_sf() %>% - select(set) %>% - rename_geometry("geometry") - - catchments_out = filter(index_table, n != 1) %>% - inner_join(network_list$catchments, by = "id") %>% - st_as_sf() %>% - select(set) %>% - union_polygons('set') %>% - mutate(areasqkm = add_areasqkm(.)) %>% + catchments_out = catchments_out %>% bind_rows(single_catchments) %>% left_join(set_topo_fin, by = "set") %>% select(id = set, toid = toset) - mps = suppressWarnings({ - catchments_out %>% - st_cast("MULTIPOLYGON") %>% - st_cast("POLYGON") %>% - fast_validity_check() %>% - add_count(id) - }) - - if(nrow(mps) > nrow(catchments_out)){ - fixers = filter(mps, n > 1) %>% - mutate(areasqkm = add_areasqkm(.), - tmpID = 1:n()) %>% - group_by(id) - - keep = slice_max(fixers, areasqkm) %>% - ungroup() - - to_merge = filter(fixers, !tmpID %in% keep$tmpID) %>% - ungroup() - - good_to_go = filter(mps, n == 1) %>% - bind_rows(select(keep, id, toid)) %>% - fast_validity_check() - - catchments_out = suppressWarnings({ - terra::combineGeoms(vect(good_to_go), vect(to_merge)) %>% - st_as_sf() %>% - st_cast("POLYGON") + #### + + if(!is.null(catchments_out)){ + mps = suppressWarnings({ + catchments_out %>% + st_cast("MULTIPOLYGON") %>% + st_cast("POLYGON") %>% + fast_validity_check() %>% + add_count(id) }) - } - + + if(nrow(mps) > nrow(catchments_out)){ + fixers = filter(mps, n > 1) %>% + mutate(areasqkm = add_areasqkm(.), + tmpID = 1:n()) %>% + group_by(id) + + keep = slice_max(fixers, areasqkm) %>% + ungroup() + + to_merge = filter(fixers, !tmpID %in% keep$tmpID) %>% + ungroup() + + good_to_go = filter(mps, n == 1) %>% + bind_rows(select(keep, id, toid)) %>% + fast_validity_check() + + catchments_out = suppressWarnings({ + terra::combineGeoms(vect(good_to_go), vect(to_merge)) %>% + st_as_sf() %>% + st_cast("POLYGON") + }) + } + } + prepare_network(network_list = list(flowpaths = flowpaths_out, catchments = catchments_out)) } diff --git a/R/aggregate_to_distribution.R b/R/aggregate_to_distribution.R index 0060648..9278e36 100644 --- a/R/aggregate_to_distribution.R +++ b/R/aggregate_to_distribution.R @@ -31,7 +31,8 @@ aggregate_to_distribution = function(gpkg = NULL, vpu = NULL, flowpath = NULL, divide = NULL, - hydrolocations = NULL, + crs = 5070, + pois = NULL, ideal_size_sqkm = 10, min_length_km = 1, min_area_sqkm = 3, @@ -73,27 +74,37 @@ aggregate_to_distribution = function(gpkg = NULL, network_list <- read_hydrofabric(gpkg, catchments = divide, flowpaths = flowpath, - crs = 5070) |> + crs = crs) |> prepare_network() |> add_network_type(verbose = FALSE) + if(!"member_comid" %in% names(network_list$flowpaths)){ + network_list$flowpaths$member_comid = network_list$flowpaths$id + } + # Add outlets - if (!is.null(hydrolocations)) { + if (!is.null(pois)) { + + names(pois) = tolower(names(pois)) - names(hydrolocations) = tolower(names(hydrolocations)) + bad = filter(pois, !id %in% network_list$flowpaths$id) - outflows = hydrolocations %>% + outflows = pois %>% st_drop_geometry() %>% - select(poi_id, id) %>% - filter(!is.na(poi_id)) |> - group_by(id) %>% - mutate(poi_id = paste(na.omit(poi_id), collapse = ",")) %>% - slice(1) %>% - ungroup() + select(poi_id, id) network_list$flowpaths = left_join(mutate(network_list$flowpaths, poi_id = NULL), st_drop_geometry(outflows), by = 'id') + + unique(pois$poi_id) |> length() + unique(network_list$flowpaths$poi_id) |> sort() + + if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))){ + warning("All POIs are not indexable: reference", call. = FALSE ) + } + + } else { network_list$flowpaths$poi_id = NA network_list$flowpaths$hl_uri = NA @@ -113,48 +124,70 @@ aggregate_to_distribution = function(gpkg = NULL, rm(tmp) } - if(!"member_comid" %in% names(network_list$flowpaths)){ - network_list$flowpaths$member_comid = NA + main_agg = + aggregate_along_mainstems( + network_list, + ideal_size_sqkm, + min_area_sqkm, + min_length_km, + verbose = verbose, + cache_file = cache_file + ) + + if(length(unique(pois$poi_id)) != length(unique(na.omit(main_agg$flowpaths$poi_id)))){ + warning("All POIs are not indexable: mainstem agg",call. = FALSE ) } - network_list = network_list |> - aggregate_along_mainstems( - ideal_size_sqkm, - min_area_sqkm, - min_length_km, - verbose = verbose, - cache_file = cache_file - ) |> collapse_headwaters2( + collapse_agg = collapse_headwaters2( + network_list = main_agg, min_area_sqkm, min_length_km, verbose = verbose, cache_file = cache_file) - - network_list$catchments = clean_geometry(network_list$catchments, ID = "id", keep = NULL) - -if(!is.null(hydrolocations)){ - network_list$hydrolocations = network_list$flowpaths %>% + collapse_agg$catchments = clean_geometry(collapse_agg$catchments, + ID = "id", + keep = NULL) + + if(length(unique(pois$poi_id)) != length(unique(na.omit(collapse_agg$flowpaths$poi_id)))){ + warning("All POIs are not indexable: collapse", call. = FALSE ) + } + + network_list = collapse_agg + rm(collapse_agg) + rm(main_agg) + +if(!is.null(pois)){ + network_list$pois = network_list$flowpaths %>% st_drop_geometry() %>% select(id, poi_id) %>% - filter(!is.na(poi_id)) %>% - tidyr::separate_longer_delim(poi_id, delim = ",") %>% - left_join(mutate(select(hydrolocations, -id), poi_id = as.character(poi_id)), by = "poi_id",relationship = "many-to-many") %>% - distinct() - + mutate(poi_id = as.integer(poi_id)) |> + filter(!is.na(poi_id)) %>% + left_join(select(pois, poi_id), by = "poi_id") } - network_list$flowpaths = hydroloom::add_streamorder(network_list$flowpaths) + network_list$flowpaths = suppressWarnings({ + hydroloom::add_streamorder(network_list$flowpaths) + }) network_list$flowpaths = select(network_list$flowpaths, id, toid, mainstem = levelpathid, order = stream_order, - member_comid, poi_id, hydroseq, lengthkm, - areasqkm, tot_drainage_areasqkm = tot_drainage_area, has_divide) %>% + member_comid, + poi_id, + hydroseq, + lengthkm, + areasqkm, + tot_drainage_areasqkm = tot_drainage_area, + has_divide) %>% mutate(divide_id = ifelse(id %in% network_list$catchments$id, id, NA)) + if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))){ + warning("All POIs are not indexable: final",call. = FALSE ) + } + topo = st_drop_geometry(network_list$flowpaths) %>% select(divide_id, toid) @@ -178,13 +211,18 @@ if(!is.null(hydrolocations)){ tot_drainage_areasqkm) %>% separate_longer_delim(col = 'member', delim = ",") %>% mutate(hf_id_part = sapply( strsplit(member, "[.]"), FUN = function(x){ x[2] }), - hf_id_part = ifelse(is.na(hf_id_part), 1L, as.integer(hf_id_part)), + hf_id_part = ifelse(is.na(hf_id_part), 1L, floor(as.numeric((hf_id_part)))), hf_id = sapply( strsplit(member, "[.]"), FUN = function(x){ as.numeric(x[1]) }), member = NULL, hf_source = "NHDPlusV2" ) %>% - left_join(st_drop_geometry(select(network_list$divides, divide_id, type, ds_id)), by = "divide_id") - + left_join(st_drop_geometry(select(network_list$divides, divide_id, type, ds_id)), by = "divide_id") |> + distinct() + + if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$network$poi_id)))){ + warning("All POIs are not indexable: network", call. = FALSE ) + } + if(!is.null(vpu)){ network_list$network$vpu = vpu } else { @@ -195,7 +233,6 @@ if(!is.null(hydrolocations)){ if(!all(st_geometry_type(network_list$divides) == "POLYGON")){ warning("MULTIPOLYGONS FOUND VPU: ", vpu) } - if (!is.null(outfile)) { diff --git a/R/catchment_geometry.R b/R/catchment_geometry.R index f41ab2e..b1b9e05 100644 --- a/R/catchment_geometry.R +++ b/R/catchment_geometry.R @@ -159,8 +159,10 @@ clean_geometry <- function(catchments, extra_parts = filter(polygons, n != 1) # dissolve, and explode if necessary - try( - extra_parts <- ms_explode(ms_dissolve(extra_parts, ID, copy_fields = names(extra_parts), sys = sys), sys = sys), + extra_parts <- try({ + extra_parts <- ms_dissolve(extra_parts, ID, copy_fields = names(extra_parts), sys = sys) + ms_explode(extra_parts, sys = sys) + }, silent = TRUE ) @@ -177,10 +179,9 @@ clean_geometry <- function(catchments, l = lengths(imap) df = data.frame( - tmpID = rep(extra_parts[['tmpID']], times = l), - uid = rep(extra_parts[[ID]], times = l), - touch_id = fl[[fl_ID]][unlist(imap)] - ) %>% + tmpID = rep(extra_parts[['tmpID']], times = l), + uid = rep(extra_parts[[ID]], times = l), + touch_id = fl[[fl_ID]][unlist(imap)] ) %>% group_by(tmpID) %>% summarize(prime = any(uid == touch_id)) @@ -192,8 +193,7 @@ clean_geometry <- function(catchments, } # recalculate area - extra_parts <- mutate(extra_parts, - areasqkm = add_areasqkm(extra_parts)) %>% + extra_parts <- mutate(extra_parts, areasqkm = add_areasqkm(extra_parts)) %>% arrange(desc(prime), desc(areasqkm)) %>% mutate(newID = row_number()) @@ -208,7 +208,9 @@ clean_geometry <- function(catchments, extra_parts %>% filter(!newID %in% main_parts$newID) - if(!sum(nrow(main_parts)) + nrow(filter(polygons, n == 1)) == MASTER_COUNT){ stop() } + if(!sum(nrow(main_parts)) + nrow(filter(polygons, n == 1)) == MASTER_COUNT){ + stop() + } main_parts = bind_rows(main_parts, filter(polygons, n == 1)) @@ -218,8 +220,14 @@ clean_geometry <- function(catchments, if(nrow(small_parts) > 0){ # dissolve, and explode if necessary - small_parts <- tryCatch( - ms_explode(ms_dissolve(small_parts, ID, copy_fields = names(small_parts), sys = sys), sys = sys), + small_parts <- tryCatch({ + tmp <- ms_dissolve(small_parts, + ID, + copy_fields = names(small_parts), + sys = sys) + + ms_explode(tmp, sys = sys) + }, error = function(e){ NULL}, warning = function(w){ NULL } ) @@ -279,7 +287,6 @@ clean_geometry <- function(catchments, in_cat = fast_validity_check(main_parts) } - if(!is.null(keep)){ in_cat = tryCatch({ @@ -295,12 +302,12 @@ clean_geometry <- function(catchments, x = select(catchments, -any_of(c('areasqkm'))) %>% st_drop_geometry() %>% - mutate(ID = as.numeric(ID)) + mutate("{ID}" := as.numeric(get(!!ID))) x2 = mutate(in_cat, areasqkm = add_areasqkm(in_cat)) |> st_transform(crs) |> - mutate(ID = as.numeric(ID)) |> - select("{ID}" := ID, areasqkm) |> + mutate("{ID}" := as.numeric(get(!!ID))) |> + select(ID, areasqkm) |> left_join(x, by = ID) return(x2) @@ -309,7 +316,11 @@ clean_geometry <- function(catchments, if(!is.null(keep)){ polygons = tryCatch({ - simplify_process(polygons, keep, sys, force = force, gb = gb) + simplify_process(polygons, + keep = keep, + sys = sys, + force = force, + gb = gb) }, error = function(e){ if(force){ message("Even when using the system mapshaper, an error has been found.") @@ -363,7 +374,7 @@ simplify_process = function(catchments, keep, sys, gb = 8, force = TRUE){ unlink(tmp2) } else { - catchments = ms_simplify(catchments, keep = keep, keep_shapes = TRUE, sys = sys) + cats = ms_simplify(catchments, keep = keep, keep_shapes = TRUE, sys = sys) } tt = fast_validity_check(cats) diff --git a/R/collapse_inward.R b/R/collapse_inward.R index 94a788c..519fa34 100644 --- a/R/collapse_inward.R +++ b/R/collapse_inward.R @@ -10,6 +10,11 @@ define_touch_id = function(flowpaths, term_cut = 1e9){ tmp = st_cast(flowpaths, "MULTILINESTRING") + st_is_empty(tmp) + + flowpaths[786,] |> + t() + ends = rename_geometry(tmp, 'geometry') %>% mutate(geometry = get_node(., "end")$geometry) @@ -30,7 +35,7 @@ define_touch_id = function(flowpaths, term_cut = 1e9){ toid = rep(ends2$toid, times = lengths(tmap)), type = rep(ends2$type, times = lengths(tmap)), hs = rep(ends2$hydroseq, times = lengths(tmap)), - poi_id = rep(ends2$poi_id, times = lengths(tmap)), + poi_id = rep(ends2$poi_id, times = lengths(tmap)), touches = tmp$id[unlist(tmap)], touches_toID = tmp$toid[unlist(tmap)], touches_hs = tmp$hydroseq[unlist(tmap)] diff --git a/R/hyaggregate_utils.R b/R/hyaggregate_utils.R index b367413..d108b86 100644 --- a/R/hyaggregate_utils.R +++ b/R/hyaggregate_utils.R @@ -39,13 +39,23 @@ prepare_network = function(network_list) { } # Add a hydrosequence to the flowpaths - network_list$flowpaths = add_hydroseq(flowpaths = network_list$flowpaths) + + network_list$flowpaths = network_list$flowpaths |> + select(-any_of(c('topo_sort', 'hydroseq'))) + + network_list$flowpaths = hydroloom::add_topo_sort(network_list$flowpaths) |> + rename(hydroseq = topo_sort) # Add area and length measures to the network list network_list = add_measures(network_list$flowpaths, network_list$catchments) network_list$flowpaths = mutate(network_list$flowpaths, areasqkm = ifelse(is.na(areasqkm), 0, areasqkm)) #network_list$flowpaths$order = get_streamorder(st_drop_geometry(mutate(select(network_list$flowpaths, ID = id, toID = toid), divergence = 0))) - network_list$flowpaths$tot_drainage_area = calculate_total_drainage_area(select(network_list$flowpaths, ID = id, toID = toid, area = areasqkm)) + + + network_list$flowpaths$tot_drainage_area = accumulate_downstream(network_list$flowpaths, "areasqkm") + + calculate_total_drainage_area(rename(network_list$flowpaths, + area = areasqkm)) check_network_validity(network_list) } diff --git a/R/hydrofabric_io.R b/R/hydrofabric_io.R index 6ba97dd..f919d86 100644 --- a/R/hydrofabric_io.R +++ b/R/hydrofabric_io.R @@ -200,24 +200,29 @@ write_hydrofabric = function(network_list, dm = names(hf_dm$flowlines), outfile, "flowpaths") + write_dm_model(data = network_list$divides, dm = names(hf_dm$divides), outfile, "divides") - write_dm_model(data = network_list$hydrolocations, - dm = names(hf_dm$hydrolocations), + write_dm_model(data = network_list$pois, + dm = names(hf_dm$pois), outfile, - "hydrolocations") - write_dm_model(data = network_list$hydrolocations_lookup, dm = names(hf_dm$hydrolocation_lookup), outfile, "hydrolocations_lookup") + "pois") + + # write_dm_model(data = network_list$pois, + # dm = names(hf_dm$hydrolocation_lookup), + # outfile, + # "hydrolocations_lookup") write_dm_model(data = network_list$network, dm = names(hf_dm$network), outfile, "network") - write_dm_model(data = network_list$network_lookup, - dm = names(hf_dm$network_lookup), - outfile, "network_lookup") + # write_dm_model(data = network_list$network_lookup, + # dm = names(hf_dm$network_lookup), + # outfile, "network_lookup") if("WB" %in% names(network_list)){ write_dm_model(data = network_list$WB, dm = names(hf_dm$WB), outfile, "WB") diff --git a/R/merge_hydrofabrics.R b/R/merge_hydrofabrics.R index d33b47b..f33a5b8 100644 --- a/R/merge_hydrofabrics.R +++ b/R/merge_hydrofabrics.R @@ -162,7 +162,7 @@ assign_global_identifiers <- function(gpkgs = NULL, network_layer = "network", overwrite = FALSE, term_add = 1e9, - modifications = NULL, + modifications = NULL, verbose = TRUE) { @@ -182,10 +182,13 @@ assign_global_identifiers <- function(gpkgs = NULL, meta$outfiles = outfiles } - modifications = modifications %>% - select(from, to) %>% - mutate(connection = 1:n()) %>% - pivot_longer(-connection, names_to = "type", values_to = 'hf_id') + if(!is.null(modifications)){ + modifications = modifications %>% + select(from, to) %>% + mutate(connection = 1:n()) %>% + pivot_longer(-connection, names_to = "type", values_to = 'hf_id') + } + ll = lapply( 1:nrow(meta), @@ -203,15 +206,22 @@ assign_global_identifiers <- function(gpkgs = NULL, hyaggregate_log("FATAL", glue("Some terminal not found."), verbose) } - if(nrow(filter(ll, !is.na(hf_id))) != nrow(modifications)){ - hyaggregate_log("FATAL", glue("Some modification connections not found."), verbose) + if(!is.null(modifications)){ + if(nrow(filter(ll, !is.na(hf_id))) != nrow(modifications)){ + hyaggregate_log("FATAL", glue("Some modification connections not found."), verbose) + } + + conn = filter(ll, !is.na(connection)) %>% + arrange(connection) %>% + select(newID, connection, type) %>% + pivot_wider(id_cols = connection, names_from = type, values_from = newID) %>% + select(connection, from, to) + + } else { + conn = NULL } - conn = filter(ll, !is.na(connection)) %>% - arrange(connection) %>% - select(newID, connection, type) %>% - pivot_wider(id_cols = connection, names_from = type, values_from = newID) %>% - select(connection, from, to) + for(i in 1:nrow(meta)){ lyrs = st_layers(meta$path[i])$name diff --git a/R/modify_collapse_headwaters.R b/R/modify_collapse_headwaters.R index a84bd1c..c98aae0 100644 --- a/R/modify_collapse_headwaters.R +++ b/R/modify_collapse_headwaters.R @@ -52,11 +52,10 @@ build_headwater_collapse = function(network_list, collapse_headwaters2 = function(network_list, - min_area_sqkm = 3, - min_length_km = 1, - verbose = TRUE, - cache_file = NULL) { - + min_area_sqkm = 3, + min_length_km = 1, + verbose = TRUE, + cache_file = NULL) { hyaggregate_log("INFO", "\n--- Collapse Network Inward ---\n", verbose) diff --git a/R/nexus_topology.R b/R/nexus_topology.R new file mode 100644 index 0000000..4c4bf12 --- /dev/null +++ b/R/nexus_topology.R @@ -0,0 +1,580 @@ +add_prefix = function(topo, + hf_prefix = "cat-", + nexus_prefix = "nex-", + terminal_nexus_prefix = "tnx-", + coastal_nexus_prefix = "cnx-", + internal_nexus_prefix = "inx-"){ + + + t1 = filter(topo, topo_type == "network") %>% + mutate(nex_pre = case_when( + type == "terminal" ~ terminal_nexus_prefix, + type == "coastal" ~ coastal_nexus_prefix, + type == "internal" ~ internal_nexus_prefix, + TRUE ~ nexus_prefix + )) %>% + mutate(id = paste0(hf_prefix, id), + toid = paste0(nex_pre, toid), + nex_pre = NULL) + + t2 = filter(topo, topo_type == "nexus") %>% + mutate(nex_pre = case_when( + type == "terminal" ~ terminal_nexus_prefix, + type == "coastal" ~ coastal_nexus_prefix, + type == "internal" ~ internal_nexus_prefix, + TRUE ~ nexus_prefix + )) %>% + mutate(id = paste0(nex_pre,id), + toid = paste0(hf_prefix, toid), + nex_pre = NULL) + + bind_rows(t1,t2) %>% + select(id, toid, type, contains("vpu")) + +} + + + +flush_prefix = function (input, col) { + for (i in col) { + input[[i]] = as.numeric(gsub(".*-", "", input[[i]])) + } + input +} + +add_nonnetwork_nexus_location = function(divides, + coastal_nexus_prefix = "cnx-", + internal_nexus_prefix = "inx-", + waterbody_prefix = "wb-"){ + + coastal_nex = suppressWarnings({ st_point_on_surface(filter(divides, type == "coastal") ) }) %>% + flush_prefix(c('divide_id', 'toid')) %>% + mutate(divide_id = paste0(coastal_nexus_prefix, divide_id), + type = "coastal", + toid = paste0(waterbody_prefix, 0)) %>% + select(divide_id, toid, type) %>% + rename_geometry("geometry") + + internal_nex = suppressWarnings({ st_point_on_surface(filter(divides, type == "internal") ) }) %>% + flush_prefix(c('divide_id', 'toid')) %>% + mutate(divide_id = paste0(internal_nexus_prefix, divide_id), + type = "internal", + toid = paste0(waterbody_prefix, 0)) %>% + select(divide_id, toid, type) %>% + rename_geometry("geometry") + + return(bind_rows(coastal_nex, internal_nex) %>% + rename(id = divide_id)) + +} + + +#' Realign Topology to a nexus network +#' @param network_list list containing flowpath and catchment `sf` objects +#' @param nexus_prefix character prefix for nexus IDs +#' @param terminal_nexus_prefix character prefix for terminal nexus IDs +#' @param coastal_nexus_prefix character prefix for coastal nexus IDs +#' @param internal_nexus_prefix character prefix for internal nexus IDs +#' @param catchment_prefix character prefix for catchment IDs +#' @param waterbody_prefix character prefix for catchment IDs +#' @return list +#' @export +#' @importFrom dplyr select mutate left_join bind_rows group_by mutate ungroup filter distinct bind_rows slice_max rename case_when +#' @importFrom nhdplusTools rename_geometry get_node +#' @importFrom sf st_drop_geometry st_intersects st_as_sf +# +# filter(network_list$flowpaths, id == 3548) %>% +# mutate(l = add_lengthkm(.)) |> +# mapview::mapview() +# +# filter(topo, id == 3548) + +realign_topology = function(network_list, + nexus_prefix = NULL, + terminal_nexus_prefix = NULL, + coastal_nexus_prefix = NULL, + internal_nexus_prefix = NULL, + catchment_prefix = NULL, + waterbody_prefix = NULL, + term_add = 1e9, + term_filter = NULL){ + + + net = select(network_list$flowpaths, id, toid) |> + st_drop_geometry() + + if(is.null(term_filter)){ + + if(is.null(term_add)){ + term_net = filter(net, toid == 0 ) + } else { + term_net = filter(net, toid >= term_add) + } + + ee2 = filter(network_list$flowpaths, id %in% term_net$id) + + net = bind_cols(id = ee2$id, toid = ee2$toid, st_coordinates(get_node(ee2, "end"))) |> + group_by(X, Y) |> + mutate(toid = cur_group_id() + term_add) |> + ungroup() |> + bind_rows(filter(net, toid < term_add)) + + network_list$flowpaths$toid = net$toid[match(network_list$flowpaths$id, net$id)] + network_list$catchments$toid = net$toid[match(network_list$catchments$id, net$id)] + } else { + term_add = term_filter + } + + # Isolate the flow network + iso = select(network_list$flowpaths, id, toid, hydroseq, poi_id) + + # Cast flow network to end nodes, these are the outlets of the flowpaths + ends = iso %>% + rename_geometry('geometry') %>% + mutate(geometry = get_node(., "end")$geometry) %>% + left_join(st_drop_geometry(select(iso, id)), by = "id", multiple = "all") + + # Get all start and end node geometries + starts_ends = bind_rows(get_node(iso, "start"), get_node(iso, "end")) + # write_sf(starts_ends, "test.gpkg") + # Find the locations where the end points intersect with starting/ending points + emap = st_intersects(ends, starts_ends) + + # If more then one intersection occurs its a nexus, + # otherwise it is a junction + ends$type = ifelse(lengths(emap) > 1, "nex", "jun") + + # Now, intersect the typed ends with the isolated flow network + tmap = st_intersects(ends, iso) + + # Build a data.frame that stores the following: + # 1. ID - the ID of the end node + # 2. type - the type of the end node + # 3. touches - the flowlines the end node touches + # 4. hs - the hydrosequence of the end node + # 5. touches_toID - the toID of the flowline touched by the endnode + + # The data.frame is grouped by the ID (1) and the total entries are tabulated. + # The data.frame is then joined to the isolated flownetwork to append the topologic toID + + df = data.frame( + id = rep(ends$id, times = lengths(tmap)), + type = rep(ends$type, times = lengths(tmap)), + hs = rep(ends$hydroseq, times = lengths(tmap)), + touches = iso$id[unlist(tmap)], + touches_toID = iso$toid[unlist(tmap)]) %>% + group_by(id) %>% + mutate(n = dplyr::n()) %>% + ungroup() %>% + left_join(st_drop_geometry(select(iso, id, toid)), by = "id", multiple = "all") + + ### --- UPDATE DF --- ### + + # tmp_topo = df %>% + # select(id, toid) %>% + # distinct() %>% + # rename(new_toid = toid) + # + # df2 = left_join(df, tmp_topo, by = "id") + # + # filter(df2, id == 1527) + # + ### --- TERMINALS --- ### + # Terminals are those where an end point only touches itself + # ID/toID topology persists + + terminals = filter(df, toid >= term_add) %>% + select(id, toid = toid) %>% + mutate(type = "terminal") + + straight_flow = filter(df, n == 1 & id == touches) |> + select(id, toid) |> + mutate(type = 'nexus') + + ### --- NEXUS 1 --- ### + # nexuses are those typed as nex where the id and touches ID are not equal + # ID/toID topology persists + + nexus = filter(df, type == "nex") %>% + filter(id != touches) %>% + select(id, toid) %>% + distinct() %>% + mutate(type = "nexus") + + # Terminals and easy nexuses make the first set of nexus locations + non_junction = bind_rows(terminals, straight_flow, nexus) + + ### --- Junctions --- ### + # When flowlines have junctions involved (e.g. more then one incoming flowline) + # we need to ensure the most upstream is selected as the nexus AND + # that the others are topologically pushed downstream + + # First, all junctions are found + juns = filter(df, type == "jun") %>% + filter(id != touches) |> + filter(!id %in% non_junction$id) + + multi_ids = filter(juns, duplicated(id))$id + + juns = filter(juns, id %in% multi_ids) %>% + filter(touches_toID == toid) %>% + bind_rows(filter(juns, !id %in% multi_ids)) + + # The "top" junctions are those that do not touch any IDs in tmp1 + # AND that have the largest hydrosequnece + # Here we shift the toID to the flowline it touches. + + top_juns = juns %>% + filter(!touches %in% non_junction$toid) %>% + group_by(touches) %>% + slice_max(hs) %>% + ungroup() %>% + select(id, toid = touches) %>% + distinct() %>% + mutate(type = "nexus") + + inner_term = juns %>% + filter(!id %in% top_juns$id) %>% + filter(touches_toID > term_add) %>% + left_join( + select(df, id, ds_term = touches_toID), + by = c('toid' = "id"), + relationship = "many-to-many") %>% + left_join( + select(df, toid, hs), + by = c('ds_term' = "toid"), + relationship = "many-to-many") %>% + distinct() %>% + filter(touches_toID != ds_term) %>% + group_by(id) %>% + slice_min(hs.y) %>% + select(id, toid = ds_term) %>% + mutate(type = "junction") %>% + ungroup() + + # The "inner" junctions are those that are not top junctions + # Here we shift the toID to the toID of the flowline it touches + + inner_juns = juns %>% + filter(!id %in% c(top_juns$id, inner_term$id)) %>% + select(id, toid = touches_toID) %>% + distinct() %>% + mutate(type = "junction") + + if(nrow(top_juns) + nrow(inner_juns) + nrow(inner_term) != nrow(juns)){ + stop("To many junctions produced!", call. = FALSE) + } + + # The complete nexus topo network is the combination of the first set, + # The top junctions and the inner junctions + # Collectively, these define the fl --> nex network topology + + topo = bind_rows(non_junction, top_juns, inner_term, inner_juns) %>% + mutate(topo_type = "network") |> + select(id, toid, topo_type, type) |> + distinct() + + + ## NEW!!! ## + xx = select(topo, toid, type) %>% + group_by(toid) %>% + mutate(type = ifelse(any(type == 'terminal'), "terminal", type)) %>% + slice(1) %>% + ungroup() + + topo = topo %>% + mutate(type = NULL) %>% + left_join(xx, by = "toid", multiple = "all") %>% + filter(id != toid) + + # We'll use the fl --> nex topo to modify the input flow network toIDs + # Additionally we will add the NextGen required prefixes. + fl = left_join(select(network_list$flowpaths, -toid), + topo, + by = "id", + multiple = "all") %>% + rename_geometry('geometry') %>% + mutate(id = paste0(waterbody_prefix, id), + toid = paste0(ifelse(type == "terminal", terminal_nexus_prefix, nexus_prefix), toid), + realized_catchment = gsub(waterbody_prefix, catchment_prefix, id)) |> + filter(!duplicated(id)) + + + + divide = left_join(select(network_list$catchments, -toid), + select(topo, id, toid, net_type = type), + by = "id") %>% + rename_geometry('geometry') %>% + #filter(id == 24873) |> + mutate(toid = ifelse(type %in% c("coastal", "internal"), divide_id, toid), + net_type = ifelse(is.na(net_type), type, net_type), + type = ifelse(net_type == "terminal", "terminal", type)) %>% + mutate(nex_pre = case_when( + type == "terminal" ~ terminal_nexus_prefix, + type == "coastal" ~ coastal_nexus_prefix, + type == "internal" ~ internal_nexus_prefix, + TRUE ~ nexus_prefix + )) %>% + mutate(divide_id = paste0(catchment_prefix, divide_id), + id = paste0(waterbody_prefix, id), + toid = paste0(nex_pre, toid), + topo_type = NULL, + nex_pre = NULL) %>% + select(divide_id, id, ds_id, toid, type, areasqkm, has_flowline) |> + filter(!duplicated(divide_id)) + + # The nexuses defined so far are those part of the dendritic network, + # We also want to add those that are coastal or inland. + # We need to make POINT locations for + # all nexus and terminal, coastal and inland divides + + nex = filter(fl, type %in% c("nexus", "terminal")) %>% + group_by(toid) |> + #ADDED!! + slice_max(order) %>% + slice_max(hydroseq) %>% + ungroup() %>% + mutate(geometry = get_node(., "end")$geometry, + id = toid) %>% + select(id, toid, poi_id, type = type) %>% + flush_prefix(col = c("id", "toid")) %>% + distinct() %>% + mutate(nex_pre = case_when( + type == "terminal" ~ terminal_nexus_prefix, + type == "coastal" ~ coastal_nexus_prefix, + type == "internal" ~ internal_nexus_prefix, + TRUE ~ nexus_prefix + )) %>% + mutate(id = paste0(nex_pre, id), + toid = paste0(waterbody_prefix, toid), + nex_pre = NULL) %>% + mutate(toid = ifelse(type == "terminal", paste0(waterbody_prefix, 0), toid)) %>% + bind_rows( + add_nonnetwork_nexus_location( + divide, + coastal_nexus_prefix = coastal_nexus_prefix, + internal_nexus_prefix = internal_nexus_prefix, + waterbody_prefix = waterbody_prefix + )) + + # We then add the nex --> fl topo to the existing fl --> nex topo + topo = suppressWarnings({ + nex %>% + flush_prefix(col = c('id', 'toid')) %>% + select(id, toid, type) %>% + st_drop_geometry() %>% + mutate(topo_type = "nexus") %>% + bind_rows(topo) %>% + distinct() + }) + + + + if(sum(!divide$toid %in% nex$id) != 0){ + stop('Divides flow to nexus locations that do not exist!\n', + paste(st_drop_geometry(divide[which(!divide$toid %in% nex$id),"divide_id"]),"-->", st_drop_geometry(divide[which(!divide$toid %in% nex$id),"toid"]), "\n"), call. = FALSE) + } + + filter(divide, divide_id == 'cat-212673') + filter(divide, divide_id == 'cat-212675') + + if(sum(!fl$toid %in% nex$id) != 0 ){ + stop('Flowpaths flow to nexus locations that do not exist!\n', + paste(st_drop_geometry(fl[which(!fl$toid %in% nex$id),"id"]),"-->", st_drop_geometry(fl[which(!fl$toid %in% nex$id),"toid"]), "\n"), call. = FALSE) + } + + if(sum(duplicated(divide$divide_id)) != 0 ){ + stop('Divides have duplicated IDs!!') + } + + if(sum(duplicated(fl$id)) != 0 ){ + stop('Flowlines have duplicated IDs!!') + } + + n = list(flowpaths = select(fl, -type, -topo_type), + divides = divide, + nexus = nex, + topo = topo + ) + + #write_hydrofabric(n, "data/test_16.gpkg", enforce_dm = FALSE) + + return(n) + +} + + +#' @title Apply Nexus Topology +#' @description This function enforces the nexus-->flowpath topology and adds nexus locations, +#' a catchment edge list, a flowpath edge list, and a lookup_table to the +#' network_list object. +#' @param network_list list containing flowpath and catchment `sf` objects +#' @param nexus_prefix character prefix for nexus IDs +#' @param terminal_nexus_prefix character prefix for terminal nexus IDs +#' @param coastal_nexus_prefix character prefix for coastal nexus IDs +#' @param internal_nexus_prefix character prefix for internal nexus IDs +#' @param catchment_prefix character prefix for catchment IDs +#' @param waterbody_prefix character prefix for catchment IDs +#' @param enforce_dm should the data model be validated prior to writing? +#' @param export_gpkg file path to write new data. If NULL list object is returned +#' @return list or file path +#' @importFrom sf read_sf st_point_on_surface +#' @importFrom dplyr select mutate left_join everything distinct contains slice n +#' @export + +apply_nexus_topology = function(gpkg, + vpu = NULL, + nexus_prefix = "nex-", + terminal_nexus_prefix = "tnx-", + coastal_nexus_prefix = "cnx-", + internal_nexus_prefix = "inx-", + catchment_prefix = "cat-", + waterbody_prefix = "wb-", + term_add = 1e9, + term_filter = NULL, + verbose = TRUE, + enforce_dm = FALSE, + export_gpkg = NULL){ + + hyaggregate_log("INFO", "\n--- Applying HY_feature topology ---\n", verbose) + + network_list = read_hydrofabric(gpkg, verbose = verbose) + + ngen_flows = realign_topology(network_list, + nexus_prefix = nexus_prefix, + terminal_nexus_prefix = terminal_nexus_prefix, + coastal_nexus_prefix = coastal_nexus_prefix, + internal_nexus_prefix = internal_nexus_prefix, + catchment_prefix = catchment_prefix, + waterbody_prefix = waterbody_prefix, + term_add = term_add, + term_filter = term_filter) + + if(!is.null(term_filter)){ + term_add = term_filter + } + + if(!is.null(vpu)) { ngen_flows$topo$vpu = as.character(vpu) } + + ngen_flows$flowpaths = ngen_flows$flowpaths %>% + select(id, toid, + mainstem, order, hydroseq, + lengthkm, areasqkm, tot_drainage_areasqkm, + has_divide, divide_id = realized_catchment, + poi_id, + contains("vpu")) + + ngen_flows$divides = ngen_flows$divides %>% + select(divide_id, toid, type, ds_id, areasqkm, contains("vpu")) %>% + left_join(st_drop_geometry(select(ngen_flows$flowpaths, id, divide_id, lengthkm, tot_drainage_areasqkm )), by = "divide_id") %>% + mutate(areasqkm = add_areasqkm(.), + has_flowline = !is.na(id)) + + + if(layer_exists(gpkg, "pois")){ + + x = read_sf(gpkg, "pois") |> + select(poi_id) |> + mutate(poi_id = as.integer(poi_id)) + + y = st_drop_geometry(ngen_flows$flowpaths) |> + select(poi_id, id, nex_id = toid) |> + mutate(poi_id = as.integer(poi_id)) + + ngen_flows$pois = left_join(x,y, by = "poi_id") + } + + + network = add_prefix(ngen_flows$topo, + hf_prefix = waterbody_prefix, + nexus_prefix = nexus_prefix, + terminal_nexus_prefix = terminal_nexus_prefix, + coastal_nexus_prefix = coastal_nexus_prefix, + internal_nexus_prefix = internal_nexus_prefix) %>% + full_join(st_drop_geometry(select(ngen_flows$flowpaths, + has_divide, + id, divide_id, + lengthkm, + tot_drainage_areasqkm, + mainstem, + has_divide, + contains("vpu"))), + by = "id", + multiple = "all") %>% + full_join(st_drop_geometry(select(ngen_flows$divides, divide_id, areasqkm, n2 = type)), by = "divide_id") %>% + mutate(type = ifelse(is.na(type), n2, type), n2 = NULL) + + if(layer_exists(gpkg, "pois")){ + t = st_drop_geometry(ngen_flows$pois) |> + select(id, poi_id) |> + distinct() + + network = left_join(network, t, by = "id", relationship = "many-to-many") + } + + if(layer_exists(gpkg, "network")){ + tmp = read_sf(gpkg, "network") %>% + mutate(id = paste0(waterbody_prefix, id), + divide_id = paste0(catchment_prefix, divide_id), + #hl_uri = as.character(hl_uri), + toid = NULL, + t = rowSums(is.na(.))) %>% + group_by(id) %>% + slice_min(t) %>% + ungroup() + + sinks = filter(network, type %in% c(c('coastal', "internal"))) %>% + select(divide_id) %>% + filter(complete.cases(.)) %>% + left_join(select(st_drop_geometry(ngen_flows$divides), divide_id, toid), by = 'divide_id') %>% + left_join(select(tmp, -id), by = 'divide_id') %>% + mutate(vpu = as.character(vpu), + poi_id = as.integer(poi_id)) + + + ngen_flows$network = filter(network, !type %in% c(c('coastal', "internal"))) %>% + left_join(select(tmp,id, hf_source, hf_id, hf_id_part, hydroseq), by = 'id', relationship = "many-to-many") %>% + bind_rows(sinks) %>% + select(id, toid, divide_id, ds_id, mainstem, poi_id, hydroseq, poi_id, #hl_uri, + hf_source, hf_id, lengthkm, areasqkm, tot_drainage_areasqkm, + type, vpu) + } + + # if(layer_exists(gpkg, "hydrolocations")){ + # ngen_flows$nexus = ngen_flows$nexus %>% + # left_join(select(st_drop_geometry(ngen_flows$hydrolocations), #hl_uri, + # poi_id), + # by = "poi_id", + # relationship = "many-to-many") %>% + # mutate(type = type, type = NULL) %>% + # group_by(id) %>% + # #mutate(hl_uri = paste(hl_uri, collapse = ",")) %>% + # distinct() %>% + # ungroup() %>% + # mutate(type = case_when( + # !is.na(poi_id) ~ "poi", + # TRUE ~ "nexus" + # )) + # } + + ngen_flows$topo = NULL + + ## Add Nels Check ## + ## >>> ndf[ ~ndf['id'].isin( cdf['toid'] ) ] + baddies = filter(ngen_flows$nexus, !id %in% ngen_flows$divides$toid) %>% + filter(toid != "wb-0") + + write_sf(baddies, export_gpkg, "error") + + if(nrow(baddies) > 0){ + message("Foiled again ... ID/toIDs: \n\t", + paste(paste0(baddies$id, "-->", baddies$toid), collapse = "\n\t")) + } + + if(!is.null(export_gpkg)){ + write_hydrofabric(ngen_flows, export_gpkg, enforce_dm = enforce_dm) + return(export_gpkg) + } else { + return(ngen_flows) + } +} diff --git a/R/prep_split_events.R b/R/prep_split_events.R index 4e2328d..c7b3840 100644 --- a/R/prep_split_events.R +++ b/R/prep_split_events.R @@ -6,25 +6,24 @@ #' @return sf POINT object #' @export -prep_split_events = function(pois, fline, divides, theshold = 25) { +prep_split_events = function(pois, fline, divides, threshold = 25) { if("id" %in% names(fline)){ - fline = rename(fline, COMID = id) + fline$COMID = fline$id + fline$id = NULL flip_id = TRUE } else { flip_id = FALSE } if("frommeas" %in% names(fline)){ - fline = rename(fline, - FromMeas = frommeas, - ToMeas = tomeas) + fline$FromMeas = fline$frommeas + fline$ToMeas = fline$tomeas flip_meas = TRUE } else { flip_meas = FALSE } - if(inherits(pois, "sf")){ pts = st_filter(pois,divides) } else { @@ -42,7 +41,8 @@ prep_split_events = function(pois, fline, divides, theshold = 25) { left_join(select(st_drop_geometry(fline), COMID, FromMeas, ToMeas), by = "COMID") split_sites <- reference_poi |> - filter((100 * (REACH_meas - FromMeas) / (ToMeas - FromMeas)) > theshold) + mutate(m = (100 * (REACH_meas - FromMeas) / (ToMeas - FromMeas))) |> + filter(m < threshold) if(flip_id){ split_sites = rename(split_sites, id = COMID) diff --git a/R/reconcile.R b/R/reconcile.R index 6a08e0d..3344981 100644 --- a/R/reconcile.R +++ b/R/reconcile.R @@ -177,7 +177,7 @@ reconcile_catchment_divides <- function(catchment, simplify_tolerance_m = 40, vector_crs = 5070, fix_catchments = TRUE, - keep = .9) { + keep = NULL) { in_crs <- st_crs(catchment) catchment <- rename_geometry(catchment, "geom") @@ -187,12 +187,12 @@ reconcile_catchment_divides <- function(catchment, if(!is.null(fdr) & !is.null(fac)){ fdr_temp <- fdr + if(!inherits(fdr_temp, "SpatRaster")){ fdr_temp <- terra::rast(fdr_temp) } catchment <- st_transform(catchment, terra::crs(fdr_temp)) - #st_precision(catchment) <- terra::res(fdr_temp)[1] fline_ref <- st_transform(fline_ref, terra::crs(fdr_temp)) fline_rec <- st_transform(fline_rec, terra::crs(fdr_temp)) } @@ -200,30 +200,30 @@ reconcile_catchment_divides <- function(catchment, reconciled <- st_drop_geometry(fline_rec) %>% dplyr::select(ID, member_COMID) - rm(fline_rec) + #rm(fline_rec) # Not all catchments have flowlines. Remove the flowlines without. comid <- fline_ref$COMID # Just being explicit here. featureid <- catchment$FEATUREID # type conversion below is annoying. # as.integer removes the .1, .2, semantic part but the response retains # the semantic component. If you don't know what this means, stop hacking. - comid_with_catchment <- comid[as.integer(comid) %in% featureid] + comid_with_catchment <- comid[as.numeric(comid) %in% featureid] reconciled <- distinct(reconciled) %>% # had dups from prior steps. tidyr::separate_rows(member_COMID, sep = ",") %>% # Make long form - dplyr::filter(member_COMID %in% comid_with_catchment) %>% # + #dplyr::filter(member_COMID %in% comid_with_catchment) %>% # dplyr::group_by(ID) %>% dplyr::summarise(member_COMID = paste(member_COMID, collapse = ",")) %>% dplyr::ungroup() - fline_ref <- fline_ref[as.integer(fline_ref$COMID) %in% catchment$FEATUREID, ] + fline_ref <- fline_ref[floor(as.numeric((fline_ref$COMID))) %in% catchment$FEATUREID, ] to_split_bool <- as.numeric(fline_ref$COMID) != - as.integer(fline_ref$COMID) + floor(as.numeric((fline_ref$COMID))) to_split_ids <- fline_ref$COMID[which(to_split_bool)] - to_split_featureids <- unique(as.integer(to_split_ids)) + to_split_featureids <- unique(floor(as.numeric((to_split_ids)))) cl <- NULL @@ -245,7 +245,7 @@ reconcile_catchment_divides <- function(catchment, catchment = catchment, fdr = fdr, fac = fac, - min_area_m = min_area_m, + min_area_m = min_area_m , snap_distance_m = snap_distance_m, simplify_tolerance_m = simplify_tolerance_m, vector_crs = vector_crs, @@ -296,30 +296,36 @@ reconcile_catchment_divides <- function(catchment, ))) } - out <- st_sf(right_join(dplyr::select(split_cats, member_COMID = FEATUREID), - reconciled, - by = "member_COMID")) - - missing <- is.na(st_dimension(out$geom)) - - if (any(missing)) { + if(nrow(split_cats) > 0){ + + out <- st_sf(right_join(dplyr::select(split_cats, member_COMID = FEATUREID), + reconciled, + by = "member_COMID")) - out_mp <- filter(out, !missing) %>% - st_cast("MULTIPOLYGON") + missing <- is.na(st_dimension(out$geom)) - out <- select(catchment, member_COMID = FEATUREID) %>% - filter(member_COMID %in% unique(as.integer(out$member_COMID[missing]))) %>% - mutate(member_COMID = paste0(member_COMID, ".1")) %>% - mutate(ID = out$ID[match(member_COMID, out$member_COMID)]) %>% - select(ID, member_COMID) %>% - nhdplusTools::rename_geometry(attr(out_mp, "sf_column")) %>% - bind_rows(out_mp) - } + if (any(missing)) { + + out_mp <- filter(out, !missing) %>% + st_cast("MULTIPOLYGON") + + out <- select(catchment, member_COMID = FEATUREID) %>% + filter(member_COMID %in% unique(floor(as.numeric(out$member_COMID[missing])))) %>% + mutate(member_COMID = paste0(member_COMID, ".1")) %>% + mutate(ID = out$ID[match(member_COMID, out$member_COMID)]) %>% + select(ID, member_COMID) %>% + nhdplusTools::rename_geometry(attr(out_mp, "sf_column")) %>% + bind_rows(out_mp) + } + } + if(fix_catchments){ - # cat("Fixing Catchment Geometries...\n") - clean_geometry(catchments = out, "ID", keep) %>% - sf::st_transform(in_crs) + cat("Fixing Catchment Geometries...\n") + clean_geometry(catchments = out, + ID = "ID", + crs = in_crs, + keep = keep) } else { sf::st_transform(out, in_crs) } diff --git a/R/refactor_nhdplus.R b/R/refactor_nhdplus.R index 25406aa..66e9ffc 100644 --- a/R/refactor_nhdplus.R +++ b/R/refactor_nhdplus.R @@ -1,3 +1,11 @@ +drop_geometry <- function(x) { + if("sf" %in% class(x)) { + sf::st_drop_geometry(x) + } else { + x + } +} + #' @title Refactor NHDPlus #' @description A complete network refactor workflow has been packaged #' into this function. Builds a set of normalized catchment-flowpaths from @@ -84,14 +92,12 @@ refactor_nhdplus <- function(nhdplus_flines, flines <- nhdplus_flines %>% sf::st_cast("LINESTRING", warn = warn) %>% sf::st_transform(5070) %>% - split_flowlines(split_flines_meters, para = split_flines_cores, + split_flowlines(max_length = split_flines_meters, para = split_flines_cores, avoid = exclude_cats, events = events) rm(nhdplus_flines) - if (warn) { - message("flowlines split complete, collapsing") - } + if (warn) { message("flowlines split complete, collapsing") } exclude_cats <- c(exclude_cats, dplyr::filter(flines, !is.na(event_REACH_meas))$COMID, dplyr::filter(flines, !is.na(event_REACH_meas))$toCOMID) diff --git a/R/refactor_wrapper.R b/R/refactor_wrapper.R index 138866a..d2bbfc5 100644 --- a/R/refactor_wrapper.R +++ b/R/refactor_wrapper.R @@ -25,6 +25,10 @@ refactor = function (gpkg = NULL, split_flines_meters = 10000, collapse_flines_meters = 1000, collapse_flines_main_meters = 1000, + threshold = 25, + min_area_m = 800, + snap_distance_m = 100, + simplify_tolerance_m = 40, cores = 1, fac = NULL, fdr = NULL, @@ -33,11 +37,17 @@ refactor = function (gpkg = NULL, outfile = NULL) { network_list = read_hydrofabric(gpkg, catchments, flowpaths) + network_list$catchments = rename_geometry(network_list$catchments, 'geometry') + network_list$flow = rename_geometry(network_list$flowpaths, 'geometry') + names(network_list$catchments) = tolower(names(network_list$catchments)) + names(network_list$flowpaths) = tolower(names(network_list$flowpaths)) fl_lookup <- c(id = "comid", toid = "tocomid", levelpathi = "mainstemlp") div_lookup <- c(featureid = "divide_id", toid = "tocomid", levelpathi = "mainstemlp") - network_list$flowpaths = rename(network_list$flowpaths, any_of(fl_lookup)) + geom = st_geometry(network_list$flowpaths) + network_list$flowpaths = dplyr::rename(as.data.frame(network_list$flowpaths), any_of(fl_lookup)) + st_geometry( network_list$flowpaths ) = geom tf <- tempfile(pattern = "refactored", fileext = ".gpkg") tr <- tempfile(pattern = "reconciled", fileext = ".gpkg") @@ -47,7 +57,11 @@ refactor = function (gpkg = NULL, avoid = avoid[avoid %in% network_list$flowpaths$id] if(!is.null(pois)){ - events = prep_split_events(pois, network_list$flowpaths, network_list$catchments, 25) %>% + + events = prep_split_events(pois, + fline = network_list$flowpaths, + divides = network_list$catchments, + threshold = threshold) %>% mutate(event_identifier = as.character(row_number())) outlets <- pois %>% @@ -161,6 +175,7 @@ refactor = function (gpkg = NULL, outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets), fill = TRUE) %>% st_as_sf() + } else { if (!is.null(outlets)) { outlets_ref_COMID <- outlets %>% @@ -206,10 +221,27 @@ refactor = function (gpkg = NULL, div_lookup <- c(FEATUREID = "divide_id", FEATUREID = "featureid") - network_list$catchments = rename(network_list$catchments, any_of(div_lookup)) + geom = st_geometry(network_list$catchments) + network_list$catchments = rename(as.data.frame(network_list$catchments), any_of(div_lookup)) + st_geometry( network_list$catchments ) = geom + + r = rast(fac) + + fac_open = crop(r, + project(vect(network_list$catchments), crs(r)), + snap = "out") - fac_open = climateR::dap(URL = fac, AOI = network_list$catchments) - fdr_open = climateR::dap(URL = fdr, AOI = network_list$catchments) + r = rast(fdr) + + fdr_open = crop(r, + project(vect(network_list$catchments), crs(r)), + snap = "out") + + if(!is.null(keep)){ + fix_catchments = TRUE + } else { + fix_catchments = FALSE + } divides <- reconcile_catchment_divides( @@ -219,9 +251,12 @@ refactor = function (gpkg = NULL, fdr = fdr_open, fac = fac_open, para = cores, + min_area_m = min_area_m, + snap_distance_m = snap_distance_m, + simplify_tolerance_m = simplify_tolerance_m, cache = NULL, keep = keep, - fix_catchments = FALSE + fix_catchments = fix_catchments ) %>% rename_geometry("geometry") } else { diff --git a/R/split_flowline.R b/R/split_flowline.R index 7bbe165..c6ede86 100644 --- a/R/split_flowline.R +++ b/R/split_flowline.R @@ -47,8 +47,8 @@ split_flowlines <- function(flines, max_length = NULL, select(sf::st_set_geometry(flines, NULL), COMID, toCOMID, LevelPathI, Hydroseq, TotDASqKM), by = "COMID") - split$part <- unlist(lapply(strsplit(split$split_fID, "\\."), - function(x) x[2])) + split$part <- as.numeric(unlist(lapply(strsplit(split$split_fID, "\\."), + function(x) x[2]))) split <- group_by(split, COMID) @@ -62,14 +62,11 @@ split_flowlines <- function(flines, max_length = NULL, lead(part), sep = ".")))) split <- mutate(split, COMID = paste(COMID, part, sep = "."), - LENGTHKM = sf::st_length(sf::st_geometry(split)) / 1000) + LENGTHKM = add_lengthkm(split)) split <- sf::st_as_sf(select(split, -part, -split_fID)) - attr(split$LENGTHKM, "units") <- NULL - split[["LENGTHKM"]] <- as.numeric(split[["LENGTHKM"]]) - - remove_comid <- unique(as.integer(split[["COMID"]])) + remove_comid <- unique(round(as.numeric(split[["COMID"]]))) not_split <- dplyr::filter(select(flines, COMID, toCOMID, LENGTHKM, LevelPathI, diff --git a/man/aggregate_to_distribution.Rd b/man/aggregate_to_distribution.Rd index 51cbc28..ecea327 100644 --- a/man/aggregate_to_distribution.Rd +++ b/man/aggregate_to_distribution.Rd @@ -9,7 +9,8 @@ aggregate_to_distribution( vpu = NULL, flowpath = NULL, divide = NULL, - hydrolocations = NULL, + crs = 5070, + pois = NULL, ideal_size_sqkm = 10, min_length_km = 1, min_area_sqkm = 3, diff --git a/man/clean_geometry.Rd b/man/clean_geometry.Rd index e323681..160d4e1 100644 --- a/man/clean_geometry.Rd +++ b/man/clean_geometry.Rd @@ -12,6 +12,8 @@ clean_geometry( keep = NULL, crs = 5070, grid = 9e-04, + gb = 8, + force = FALSE, sys = NULL ) } @@ -31,6 +33,10 @@ If NULL, then no simplification will be executed.} \item{crs}{integer or object compatible with sf::st_crs coordinate reference. Should be a projection that supports area-calculations.} +\item{gb}{The amount of heap memory to be allocated when force = TRUE} + +\item{force}{should the mapshaper/mapshaper-xl binaries be used directly for simplification?} + \item{sys}{logical should the mapshaper system library be used. If NULL the system library will be used if available.} } @@ -38,11 +44,10 @@ the system library will be used if available.} sf object } \description{ -Fixes geometry issues present in catchments that originate in the -CatchmentSP layers, or from the reconcile_catchments hydrofab preocess. +Fixes geometry issues present in catchments derived from DEMs. These include, but are not limited to disjoint polygon fragments, artifacts -from the 30m DEM used to generate the catchments, and non-valid geometry topolgies. -A goal of this functions is also to provide means to reduce the data column +from the DEM used to generate the catchments, and non-valid geometry topologies. +A secondary goal of this functions is to provide a way to reduce the data column of the catchments by offering a topology preserving simplification through \code{\link[rmapshaper]{ms_simplify}}. Generally a "keep" parameter of .9 seems appropriate for the resolution of diff --git a/man/prep_split_events.Rd b/man/prep_split_events.Rd index cbf1917..0d1c109 100644 --- a/man/prep_split_events.Rd +++ b/man/prep_split_events.Rd @@ -4,16 +4,16 @@ \alias{prep_split_events} \title{Prep Split Events} \usage{ -prep_split_events(pois, fline, divides, theshold = 25) +prep_split_events(pois, fline, divides, threshold = 25) } \arguments{ \item{pois}{a set of POIs with a poi_id, X and (in 5070)} \item{divides}{a set of divides geometries (EPSG:5070)} -\item{theshold}{a percentage (0-100) a POI must be upstream before splitting} - \item{flines}{a set of flowlines geometries (EPSG:5070)} + +\item{theshold}{a percentage (0-100) a POI must be upstream before splitting} } \value{ sf POINT object diff --git a/man/reconcile_catchment_divides.Rd b/man/reconcile_catchment_divides.Rd index a7d00eb..a37ad06 100644 --- a/man/reconcile_catchment_divides.Rd +++ b/man/reconcile_catchment_divides.Rd @@ -17,7 +17,7 @@ reconcile_catchment_divides( simplify_tolerance_m = 40, vector_crs = 5070, fix_catchments = TRUE, - keep = 0.9 + keep = NULL ) } \arguments{ diff --git a/man/refactor.Rd b/man/refactor.Rd index cf241a6..058602b 100644 --- a/man/refactor.Rd +++ b/man/refactor.Rd @@ -13,6 +13,10 @@ refactor( split_flines_meters = 10000, collapse_flines_meters = 1000, collapse_flines_main_meters = 1000, + threshold = 25, + min_area_m = 800, + snap_distance_m = 100, + simplify_tolerance_m = 40, cores = 1, fac = NULL, fdr = NULL,