Skip to content

Commit

Permalink
v2.2 cleanup and modifications
Browse files Browse the repository at this point in the history
  • Loading branch information
@mikejohnson51
mikejohnson51 committed Sep 24, 2024
1 parent f4a0848 commit 69bde0d
Show file tree
Hide file tree
Showing 20 changed files with 937 additions and 195 deletions.
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
131 changes: 80 additions & 51 deletions R/aggregate_along_mainstems.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) %>%
Expand All @@ -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)
Expand Down Expand Up @@ -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))
}

113 changes: 75 additions & 38 deletions R/aggregate_to_distribution.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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 {
Expand All @@ -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)) {

Expand Down
Loading

0 comments on commit 69bde0d

Please sign in to comment.