Skip to content
Snippets Groups Projects
Commit eac9c02b authored by Bock, Andy's avatar Bock, Andy
Browse files

Cleanup, removed comment blocks, deprecated funs.

parent 3dc18f46
No related branches found
No related tags found
1 merge request!107Non-dendritic work and outlet mapping, SWIM in get data
Hide whitespace changes
Inline Side-by-side
Showing
with 46 additions and 148 deletions
workspace/R/non_dend.R
+ 46
− 148
View file @ eac9c02b
......@@ -77,12 +77,12 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
# HUC12s that match Catchment footprint
huc12_in_nhd <- nhd_wbd_int_sub %>%
# HUC_12.1 is the intersection HUC12, not the NAWQA Xwalk HUC12
group_by(HUC_12.1) %>%
# int_HUC_12 is the intersection HUC12, not the NAWQA Xwalk HUC12
group_by(HUC_12_int) %>%
# Which HUC12 has the largest intersecting area with a given catchment
slice(which.max(intArea)) %>%
ungroup() %>%
group_by(FEATUREID, HUC_12.1) %>%
group_by(FEATUREID, HUC_12_int) %>%
summarize(n_NHD = n(),
NHD_area = mean(AreaSqKM),
HUC_12_area = sum(AreaHUC12),
......@@ -95,9 +95,9 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
ungroup()
# HUC12 IDs of the HUC12s above
HUC12_in_cats <- filter(st_drop_geometry(miss_cats_sub), HUC_12 %in% huc12_in_nhd$HUC_12.1) %>% # was member_COMID/FEATUREID
select(member_COMID, HUC_12) %>%
rename(new_POI_ID = HUC_12)
HUC12_in_cats <- filter(st_drop_geometry(miss_cats_sub), HUC_12_int %in% huc12_in_nhd$HUC_12_int) %>% # was member_COMID/FEATUREID
select(member_COMID, HUC_12_int) %>%
rename(new_POI_ID = HUC_12_int)
# attribute cats POI_ID with HUC12 POI_IDs if match HUC12s
divides_poi <- divides_poi %>%
......@@ -118,18 +118,18 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
group_by(FEATUREID) %>%
filter(intArea == max(intArea)) %>%
ungroup() %>%
group_by(HUC_12.1) %>%
group_by(HUC_12_int) %>%
summarize(n_HUC12 = n(),
NHD_area = sum(AreaSqKM),
HUC_12_area = mean(AreaHUC12),
intarea = sum(intArea),
CAC = intarea / ((NHD_area - intarea) + (HUC_12_area - intarea) + intarea)) %>%
filter(CAC > CAC_num, !HUC_12.1 %in% HUC12_in_cats$new_POI_ID)
filter(CAC > CAC_num, !HUC_12_int %in% HUC12_in_cats$new_POI_ID)
# HUC12 IDs of the HUC12s above
cat_in_HUC12s <- filter(st_drop_geometry(miss_cats_sub), HUC_12 %in% nhd_in_huc12$HUC_12.1) %>%
select(member_COMID, HUC_12) %>%
rename(new_POI_ID = HUC_12)
cat_in_HUC12s <- filter(st_drop_geometry(miss_cats_sub), HUC_12_int %in% nhd_in_huc12$HUC_12_int) %>%
select(member_COMID, HUC_12_int) %>%
rename(new_POI_ID = HUC_12_int)
# attribute cats POI_ID with HUC12 POI_IDs if match HUC12s
divides_poi <- divides_poi %>%
......@@ -180,7 +180,6 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
divides_poi <- divides_poi %>%
left_join(cats_us_POI_ID, by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(new_POI_ID), new_POI_ID, POI_ID),
#POI_ID = ifelse(HUC_12.1 %in% HUC12_in_cats$HUC_12.1, HUC_12.1, POI_ID),
aggStep = ifelse(member_COMID %in% cats_us_POI_ID$member_COMID, "HUC_12", aggStep)) %>%
select(-new_POI_ID)
......@@ -204,12 +203,12 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
#' meet CAC criteria
assign_HUC10 <- function(divides, HUC12_xwalk, nhd, CAC_num){
miss_cats_sub <- filter(divides, is.na(POI_ID)) %>%
mutate(HUC_10 = substr(HUC_12, 1, 10))
mutate(HUC_10 = substr(HUC_12_int, 1, 10))
nhd <- mutate(nhd, COMID = as.character(COMID))
nhd_wbd_int_sub <- HUC12_xwalk %>%
mutate(HUC_10 = substr(HUC_12.1, 1, 10))
mutate(HUC_10 = substr(HUC_12_int, 1, 10))
HUC10_area <- nhd_wbd_int_sub %>%
distinct(HUC_10, AreaHUC12, .keep_all = T) %>%
......@@ -229,8 +228,8 @@ assign_HUC10 <- function(divides, HUC12_xwalk, nhd, CAC_num){
# HUC12s that match Catchment footprint
huc10_in_nhd <- nhd_wbd_int_sub %>%
# HUC_12.1 is the intersection HUC12, not the NAWQA Xwalk HUC12
group_by(HUC_12.1) %>%
# HUC_12_int is the intersection HUC12
group_by(HUC_12_int) %>%
# Which HUC12 has the largest intersecting area with a given catchment
slice(which.max(intArea)) %>%
ungroup() %>%
......@@ -331,7 +330,7 @@ coastal_cats <- function(divides_poi, full_nhd, vpu_HUC12){
# Extract unaggregated cats within frontal HUC12s that may not be assigned later
frontal_HUC12 <- filter(vpu_HUC12, HU_12_TYPE == "F" | HU_12_DS == "OCEAN")
frontal_HUC12_cats <- filter(divides_poi, HUC_12 %in% frontal_HUC12$HUC_12, is.na(POI_ID)) %>%
frontal_HUC12_cats <- filter(divides_poi, HUC_12_int %in% frontal_HUC12$HUC_12, is.na(POI_ID)) %>%
mutate(member_COMID = as.integer(member_COMID))
divides_miss <- rbind(coastal_cats, filter(frontal_HUC12_cats, !member_COMID %in% coastal_cats$member_COMID)) %>%
......@@ -346,10 +345,10 @@ coastal_cats <- function(divides_poi, full_nhd, vpu_HUC12){
#
# attribute cats POI_ID with HUC12 POI_IDs if match HUC12s
divides_poi <- divides_poi %>%
left_join(st_drop_geometry(divides_miss) %>% select(member_COMID, HUC_12_POI = HUC_12), by = "member_COMID") %>%
left_join(st_drop_geometry(divides_miss) %>% select(member_COMID, HUC_12_POI = HUC_12_int), by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(HUC_12_POI), HUC_12_POI, POI_ID)) %>%
mutate(aggStep = ifelse(!is.na(HUC_12_POI), "coastal", aggStep)) %>%
select(-HUC_12_POI, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12, AreaHUC12, intArea) %>%
select(-HUC_12_POI, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12_int, AreaHUC12, intArea) %>%
mutate(vpu_ID = substr(rpu, 1, 2), comid_char = as.character(member_COMID))
# FL of cats that "contribtue" to the coast "flowlines" but are not coastal themselves
......@@ -374,11 +373,11 @@ coastal_cats <- function(divides_poi, full_nhd, vpu_HUC12){
# Assign POIs for these guys
cat2coast_POI <- cat2coast %>%
inner_join(st_drop_geometry(cat_at_coast) %>% select(member_COMID, Hydroseq, HUC_12, POI_ID), by = c("DnHydroseq" = "Hydroseq")) %>%
inner_join(st_drop_geometry(cat_at_coast) %>% select(member_COMID, Hydroseq, HUC_12_int, POI_ID), by = c("DnHydroseq" = "Hydroseq")) %>%
# Assign POI_ID of attached coastal HUC if it exists
mutate(POI_ID.x = ifelse(!is.na(POI_ID.y), POI_ID.y, POI_ID.x)) %>%
mutate(POI_ID.x = ifelse(HUC_12.x == HUC_12.y & is.na(POI_ID.x), HUC_12.x, POI_ID.x)) %>%
mutate(POI_ID.x = ifelse(is.na(POI_ID.x), HUC_12.x, POI_ID.x))
mutate(POI_ID.x = ifelse(HUC_12_int.x == HUC_12_int.y & is.na(POI_ID.x), HUC_12_int.x, POI_ID.x)) %>%
mutate(POI_ID.x = ifelse(is.na(POI_ID.x), HUC_12_int.x, POI_ID.x))
# Get upstream flowlines contributing to the coast that are not assigned a POI
nhd2coast_POI <- filter(nhd2coast, COMID %in% cat2coast_POI$member_COMID.x)# & StartFlag != 1)
......@@ -406,29 +405,17 @@ coastal_cats <- function(divides_poi, full_nhd, vpu_HUC12){
left_join(st_drop_geometry(cats_us_POI_ID) %>% select(member_COMID, new_POI_ID), by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(new_POI_ID), new_POI_ID, POI_ID)) %>%
mutate(aggStep = ifelse(!is.na(new_POI_ID), "coastal", aggStep)) %>%
select(-new_POI_ID, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12, AreaHUC12, intArea) %>%
select(-new_POI_ID, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12_int, AreaHUC12, intArea) %>%
mutate(vpu_ID = substr(rpu, 1, 2), comid_char = as.character(member_COMID))
} else {
divides_poi <- divides_poi %>%
left_join(st_drop_geometry(divides_miss) %>% select(member_COMID, HUC_12_POI = HUC_12), by = "member_COMID") %>%
left_join(st_drop_geometry(divides_miss) %>% select(member_COMID, HUC_12_POI = HUC_12_int), by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(HUC_12_POI), HUC_12_POI, POI_ID)) %>%
mutate(aggStep = ifelse(!is.na(HUC_12_POI), "coastal", aggStep)) %>%
select(-HUC_12_POI, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12, AreaHUC12, intArea) %>%
select(-HUC_12_POI, member_COMID, ID, rpu, POI_ID, aggStep, HUC_12_int, AreaHUC12, intArea) %>%
mutate(vpu_ID = substr(rpu, 1, 2), comid_char = as.character(member_COMID))
}
# # Extract unaggregated cats within frontal HUC12s that may not be assigned later
# frontal_HUC12 <- filter(vpu_HUC12, HU_12_TYPE == "F" | HU_12_DS == "OCEAN")
# frontal_HUC12_cats <- filter(st_drop_geometry(divides), HUC_12 %in% frontal_HUC12$HUC_12, is.na(POI_ID)) %>%
# select(member_COMID, HUC_12_POI = HUC_12)
#
# # bring over IDs to full data frame
# divides <- divides %>%
# left_join(frontal_HUC12_cats, by = "member_COMID") %>%
# mutate(POI_ID = ifelse(!is.na(HUC_12_POI), HUC_12_POI, POI_ID),
# aggStep = ifelse(!is.na(HUC_12_POI), "coastal", aggStep)) %>%
# select(-HUC_12_POI)
divides_poi <- dissolve_holes(divides_poi)
......@@ -474,12 +461,10 @@ outlets_close <- function(nhd_full, nhd_sub, divides_poi){
st_transform(crs)
# All cats except those with same COMID as little_outlets
nonout_cats <- filter(divides_poi, !member_COMID %in% little_outlets_fin$COMID) #%>%
#mutate(ID = row_number())
nonout_cats <- filter(divides_poi, !member_COMID %in% little_outlets_fin$COMID)
# All divide catchments with aggregated or POI IDs
poi_cats <- filter(divides_poi, !is.na(POI_ID)) #%>%
#mutate(row_ID = row_number())
poi_cats <- filter(divides_poi, !is.na(POI_ID))
# cast to multlinestring sf feature for proximity analysis
cats_lines <- nonout_cats %>%
......@@ -500,9 +485,9 @@ outlets_close <- function(nhd_full, nhd_sub, divides_poi){
by_element = TRUE)), 2)) %>%
mutate(dist2hru = round(as.numeric(st_distance(., poi_cats[st_nearest_feature(., poi_cats),],
by_element = TRUE)), 2)) %>%
left_join(select(st_drop_geometry(poi_cats), member_COMID, HUC_12), by = c("COMID" = "member_COMID")) %>%
inner_join(select(st_drop_geometry(nonout_cats), member_COMID, HUC_12, POI_ID, row_ID), by = c("nearest_cat" = "member_COMID")) %>%
rename(HUC_12_outlet = HUC_12.x, HUC_12_nrHRU = HUC_12.y)
left_join(select(st_drop_geometry(poi_cats), member_COMID, HUC_12_int), by = c("COMID" = "member_COMID")) %>%
inner_join(select(st_drop_geometry(nonout_cats), member_COMID, HUC_12_int, POI_ID, row_ID), by = c("nearest_cat" = "member_COMID")) %>%
rename(HUC_12_outlet = HUC_12_int.x, HUC_12_nrHRU = HUC_12_int.y)
}
return(missing_outlets)
......@@ -516,20 +501,19 @@ perims <- function(divides_poi){
# Subset to data frame of existing divides with POI assignments
divides_wpoi <- filter(divides_poi, !is.na(POI_ID)) %>%
select(POI_ID, HUC_12, member_COMID, row_ID)
select(POI_ID, HUC_12_int, member_COMID, row_ID)
# subset divides with no POI_Id assignment
divides_nopoi <- filter(divides_poi, is.na(POI_ID)) %>%
select(POI_ID, HUC_12, member_COMID, row_ID) #%>%
#mutate(ind = row_number())
select(POI_ID, HUC_12_int, member_COMID, row_ID) #
# Buffer with 1 meter
divides_nopoi_buff <- st_buffer(divides_nopoi, 1)
divides_nopoi_int <- st_intersection(divides_nopoi_buff, divides_wpoi) %>%
mutate(perim = lwgeom::st_perimeter(.) / 2) %>%
select(POI_ID, HUC_12, member_COMID, row_ID, perim, bound_cat_POI_ID = POI_ID.1, bound_cat_HUC12 = HUC_12.1,
bound_cat = member_COMID.1, bound_cat_row = row_ID.1) %>%
select(POI_ID, HUC_12_int = HUC_12_int, member_COMID, row_ID, perim, bound_cat_POI_ID = POI_ID.1,
bound_cat_HUC12 = HUC_12_int.1, bound_cat = member_COMID.1, bound_cat_row = row_ID.1) %>%
st_drop_geometry()
return(divides_nopoi_int)
......@@ -609,11 +593,11 @@ assignable_cats <- function(inData){
# Summarize by each catchment and its shared perimeters
summary_HUC12 <- nhdsub %>%
filter(member_COMID == inCOMID, bound_cat_HUC12 == HUC_12) %>%
filter(member_COMID == inCOMID, bound_cat_HUC12 == HUC_12_int) %>%
distinct(member_COMID, bound_cat_POI_ID)
summary_HUC12_multi <- nhdsub %>%
group_by(member_COMID, HUC_12, bound_cat_HUC12, outlet_toHUC12, outlet_toPOI_ID) %>%
group_by(member_COMID, HUC_12_int, bound_cat_HUC12, outlet_toHUC12, outlet_toPOI_ID) %>%
filter(member_COMID == inCOMID, bound_cat_HUC12 == outlet_toHUC12) %>%
ungroup() %>%
distinct(member_COMID, outlet_toPOI_ID)
......@@ -665,7 +649,7 @@ assignable_cats <- function(inData){
NHD_sinks <- function(divides_poi, area_thresh, HUC12_table){
# populate HUC10 field
divides_poi$HUC_10 <- substr(divides_poi$HUC_12, 1, 10)
divides_poi$HUC_10 <- substr(divides_poi$HUC_12_int, 1, 10)
# Filter HUC12 crosswalk to HUC10s id'ed in divides and closed HUC types
HUC12_sinks <- readRDS(HUC12_table) %>%
......@@ -674,9 +658,9 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table){
# Determine which HUC12s have only one or less aggregated catcment
HUC12_poi <- filter(divides_poi, !is.na(POI_ID)) %>%
distinct(POI_ID, HUC_12) %>%
group_by(HUC_12) %>%
filter(n() < 2, HUC_12 %in% HUC12_sinks$HUC_12) %>%
distinct(POI_ID, HUC_12_int) %>%
group_by(HUC_12_int) %>%
filter(n() < 2, HUC_12_int %in% HUC12_sinks$HUC_12_int) %>%
ungroup() %>%
group_by(POI_ID) %>%
filter(n() < 2) %>%
......@@ -688,21 +672,22 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table){
# Get the WBD sinks
sink_cats <- filter(divides_poi, is.na(POI_ID), !member_COMID %in% sink_cats_nhd$member_COMID) %>%
filter(HUC_10 %in% HUC12_sinks$HUC_10 & HUC_12 %in% HUC12_sinks$HUC_12) %>%
filter(HUC_10 %in% HUC12_sinks$HUC_10 & HUC_12_int %in% HUC12_sinks$HUC_12) %>%
rbind(sink_cats_nhd) %>%
distinct()
# If no sinks in VPU return to main script
if(nrow(sink_cats) == 0){
print ("no large sinks in vpu")
return(divides_poi)
}
# Check to see if HUC12 alreay identified as POI_ID in HUC12 step
sink_cats <- sink_cats %>%
mutate(POI_ID = ifelse(HUC_12 %in% divides_poi$POI_ID, HUC_12, POI_ID))
mutate(POI_ID = ifelse(HUC_12_int %in% divides_poi$POI_ID, HUC_12_int, POI_ID))
# Aggretate sink cats to HUC12 and add area
sink_cats_agg <- union_polygons_geos(filter(sink_cats, is.na(POI_ID)), "HUC_12") %>%
sink_cats_agg <- union_polygons_geos(filter(sink_cats, is.na(POI_ID)), "HUC_12_int") %>%
mutate(agg_id = row_number())
# Cast to polygon
......@@ -766,13 +751,13 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table){
sink_cats_small$agg_id <- agg_cats_small[sink_cats_index,]$agg_id
sink_cats_nhd_table <- st_drop_geometry(sink_cats_small) %>%
inner_join(select(cat_bound, -HUC_12), by = "member_COMID") %>%
inner_join(select(cat_bound, -HUC_12_int), by = "member_COMID") %>%
distinct()
small_sinks_POI <- sink_cats_nhd_table %>%
distinct(agg_id, HUC_12, bound_cat_HUC12, bound_cat_POI_ID) %>%
distinct(agg_id, HUC_12_int, bound_cat_HUC12, bound_cat_POI_ID) %>%
#summarize(perim = sum(perim)) %>%
filter(HUC_12 == bound_cat_HUC12) %>%
filter(HUC_12_int == bound_cat_HUC12) %>%
group_by(agg_id) %>%
filter(n_distinct(bound_cat_POI_ID) == 1) %>%
ungroup()
......@@ -797,90 +782,3 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table){
}
#' Find and lump remaining areas where possible
#' @param divides_lu (sf data frame) full set of divides and catchments with no assigned POI_ID
#
#' @return (sf data frame) data frame of likely POI_ID assignments of reamining sinks
NHD_sinks_orig <- function(divides_poi){
# Subset to data frame of existing divides with POI assignments
cats2bind <- filter(divides_poi, !is.na(POI_ID)) %>%
select(POI_ID, HUC_12, member_COMID)
# subset divides with no POI_Id assignment
missing_cats <- filter(divides_poi, is.na(POI_ID))
# Summarize by HUC12
miss_cats_assign <- union_polygons_geos(missing_cats, "HUC_12")
# Cast to polygon
miss_cats_poly <- do.call(rbind, lapply(1:nrow(miss_cats_assign),
function(i){st_cast(miss_cats_assign[i,], "POLYGON")})) %>%
mutate(member_COMID = row_number()) %>%
select(-areasqkm)
miss_cats_poly <- st_compatibalize(miss_cats_poly, cats2bind)
# centroids of divides with no assigned POI_ID
cents <- st_point_on_surface(missing_cats)
# intersect centroids with the aggregated HUC12 polygons
cats_int <- st_intersects(cents, miss_cats_poly)
# Index of HRUs that missing cats may lie within with holes filled
member_COMID_feat <- sapply(cats_int, function(s) if (length(s) == 0) NA else s)
# Populate ID of HRUs where applicable
missing_cats$POI_ID_init <- miss_cats_poly[member_COMID_feat,]$member_COMID
# Build dataframe for aggregated missing areas and divides with POIs
miss_cats_neighbors <- miss_cats_poly %>%
mutate(POI_ID = NA) %>%
rbind(cats2bind) %>%
mutate(row_ID = row_number()) %>%
st_make_valid(.)
# Buffer with 1 meter
miss_cats_neighbors_buff <- st_buffer(filter(miss_cats_neighbors, is.na(POI_ID)), 1)
# Determine neighboring dvidies with assigned POIs
overlap_missing <- st_overlaps(miss_cats_neighbors_buff, miss_cats_neighbors) %>%
purrr::set_names(miss_cats_neighbors_buff$row_ID) %>%
stack() %>%
mutate(ind = as.integer(ind))
# limit the result to cats that have not been assigned POI_ID
adjacent_polys <- filter(miss_cats_neighbors, is.na(POI_ID)) %>%
left_join(overlap_missing, by = c("row_ID" = "ind")) %>%
inner_join(st_drop_geometry(miss_cats_neighbors), by = c("values" = "row_ID")) %>%
select(HUC_12_noPOI = HUC_12.x, member_COMID = member_COMID.x, member_COMID_bdHRU = member_COMID.y,
HUC_12_bdHRU = HUC_12.y, POI_ID_bdHRU = POI_ID.y, POI_ID.x)
# catchments that are partially surrounded
part_surrounded <- adjacent_polys %>%
# origin is the CAT from which the neighboring relationships are determined
group_by(POI_ID_bdHRU, member_COMID) %>%
mutate(nCats = length(unique(member_COMID)), nHRU = length(unique(POI_ID_bdHRU, na.rm = T))) %>%
distinct(member_COMID, POI_ID_bdHRU, .keep_all = T)
# Final determination of aggregated non-assigned divides
final_sinks <- part_surrounded %>%
group_by(HUC_12_noPOI, member_COMID, HUC_12_bdHRU, POI_ID_bdHRU) %>%
summarize() %>%
filter(HUC_12_noPOI == HUC_12_bdHRU) %>%
ungroup() %>%
group_by(member_COMID) %>%
filter(n() == 1) %>%
select(member_COMID, POI_ID_bdHRU) %>%
mutate(member_COMID = as.integer(member_COMID))
miss_cats_assign <- missing_cats %>%
inner_join(st_drop_geometry(final_sinks), by = c("POI_ID_init" = "member_COMID")) %>%
select(member_COMID, POI_ID_bdHRU)
# Bring over POI_ID to divides
divides_fin <- divides_lu %>%
left_join(st_drop_geometry(miss_cats_assign), by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(POI_ID_bdHRU), POI_ID_bdHRU, POI_ID),
aggStep = ifelse(!is.na(POI_ID_bdHRU), "hoho", aggStep)) %>%
return(divides_fin)
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment