diff --git a/workspace/01_Gage_Selection.Rmd b/workspace/01_Gage_Selection.Rmd index c39ee2d293d9c984177529ab2416bb343301012b..f7c247f65542cdb8e751279cac302a5a039a6572 100644 --- a/workspace/01_Gage_Selection.Rmd +++ b/workspace/01_Gage_Selection.Rmd @@ -44,7 +44,6 @@ nhd <- read_sf(data_paths$ref_fl) if(file.exists("cache/dropped_gages.csv")){ file.remove("cache/dropped_gages.csv") } - # gages II gages_ii <- read_sf(data_paths$gagesii_lyr) ref_gages <- read_sf("data/nldi/usgs_nldi_gages.geojson") @@ -124,7 +123,7 @@ if(!file.exists("temp/GagestoCheck_GAGESIII.csv")){ # Subset further by streamgages below the minimum drainage area critera gagesIII_db_fin <- gagesIII_db %>% - filter(!site_no %in% potCanalsIII$site_no, drain_area > min_da_km_gages) + filter(!site_no %in% potCanalsIII$site_no, drain_area > (min_da_km_gages * sqkm_sqmi)) # Write out gages that dropped to document. gage_document(select(potCanalsIII, site_no), "GAGES-III", "ST; Pot. Canal") @@ -176,7 +175,7 @@ gageloc_sf_full <- read_sf(data_paths$nhdplus_gdb, "Gage") %>% filter(!site_no %in% gagesIII_sf_full$site_no) %>% # Update index information with NLDI gages left_join(select(st_drop_geometry(ref_gages), site_no = id, - comid = nhdpv2_COMID, reachcode = nhdpv2_COMID, + comid = nhdpv2_COMID, reachcode = nhdpv2_REACHCODE, reach_meas = nhdpv2_REACH_measure ), by = "site_no") # Subset to locations with streamflow @@ -224,7 +223,8 @@ if(!file.exists("temp/GagestoCheck_GageLoc.csv")){ } gageloc_db_fin <- gageloc_db %>% - filter(!site_no %in% potCanals_gageloc$site_no, drain_area > min_da_km_gages) + filter(!site_no %in% potCanals_gageloc$site_no, drain_area > + (min_da_km_gages * sqkm_sqmi)) # Gages that are ST site types, but reside on canals, artificial waterbodies. gage_document(select(potCanals_gageloc, site_no), "GageLoc", "ST; Pot. Canal") diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd index b085e005a809fe6725e28b1ac62ff006647ecfbb..6c1edc83dd32cafe33f79f842a9596e18470b592 100644 --- a/workspace/02_NHD_navigate.Rmd +++ b/workspace/02_NHD_navigate.Rmd @@ -91,7 +91,8 @@ if(!"Type_Gages" %in% names(tmp_POIs)) { rename(COMID = comid) %>% filter(COMID %in% nhd$COMID) %>% inner_join(st_drop_geometry(st_drop_geometry(nhd), COMID, Divergence), by = "COMID") %>% - switchDiv(., nhd) + switchDiv(., nhd) %>% + filter(drain_area > min_da_km_gages) streamgages <- streamgages_VPU %>% group_by(COMID) %>% @@ -138,17 +139,32 @@ if(!"Type_TE" %in% names(tmp_POIs)) { nhd$VPUID <- substr(nhd$RPUID, 1, 2) } - # Read in Thermoelectric shapefile - TE_COMIDs <- readxl::read_xlsx(file.path(data_paths$TE_points_path, - "nhdv1_2_comids_all_plants_all_years.xlsx")) %>% - rename(COMID = comid) %>% - #read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>% - #mutate(COMID = as.integer(COMID)) %>% - inner_join(., select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>% + # Read in Thermoelectric shapefile for 2015 estimates + TE_v1 <- read_sf(data_paths$TE_points_path, + "2015_TE_Model_Estimates_lat.long_COMIDs") + # Read in Thermoelectric Plants csv for updated estimates + TE_v2 <- read.csv(file.path(data_paths$TE_points_path, + "te_plants.csv"), + colClasses = "character") %>% + mutate(Plant.Code = as.integer(Plant.Code)) + + # Bind hydrographic address information together and get final COMIDs + TE_COMIDs_full <- TE_v2 %>% + filter(water_source %in% c("sw_fresh", "gw_fresh", "gw_sw_mix")) %>% + left_join(select(st_drop_geometry(TE_v1), + Plant.Code = EIA_PLANT_, COMID), by = "Plant.Code") %>% + left_join(rename(st_drop_geometry(HUC12_COMIDs), huc12_comid = COMID), + by = c("huc_12" = "HUC12")) %>% + mutate(COMID = ifelse(is.na(COMID), huc12_comid, COMID)) %>% + select(-huc12_comid) + + # Prepare TE for POI Creation + TE_COMIDs <- TE_COMIDs_full %>% + inner_join(select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>% filter(grepl(paste0("^", substr(vpu_codes, 1, 2), ".*"), .data$VPUID), COMID > 0) %>% switchDiv(., nhd) %>% group_by(COMID) %>% - summarize(eia_id = paste0(unique(eia_id), collapse = " "), count = n()) %>% + summarize(eia_id = paste0(unique(Plant.Code), collapse = " "), count = n()) %>% ungroup() if(nrow(TE_COMIDs) > 0){ @@ -188,14 +204,14 @@ WBs_VPU_all <- filter(readRDS(data_paths$waterbodies_path), st_transform(geo_crs) %>% st_make_valid() -saveRDS(WBs_VPU_all, data_paths$waterbodies_path) +write_sf(WBs_VPU_all, temp_gpkg, "wbs_draft") if(needs_layer(temp_gpkg, "WBs_layer_orig")){ sf_use_s2(FALSE) ref_WB <- WBs_VPU_all %>% group_by(wb_id, GNIS_NAME) %>% sfheaders::sf_remove_holes() %>% - summarize(do_union = T, member_comid = paste(COMID, collapse = ",")) %>% + summarize(do_union = TRUE, member_comid = paste(COMID, collapse = ",")) %>% st_make_valid() %>% ungroup() %>% mutate(area_sqkm = as.numeric(st_area(.)) / 1000000) %>% @@ -220,8 +236,9 @@ if(!"Type_resops" %in% names(tmp_POIs)){ mutate(member_comid = as.character(member_comid)) %>% distinct() + # TODO: get this file from data.json # ResOpsUS locations with attributes from the VPU NHDv2 wb set - resops_wb_df <- read.csv("data/reservoir_data/resops_crosswalk_AB.csv") %>% + resops_wb_df <- read.csv("data/reservoir_data/istarf_xwalk_final_48_gfv11.csv") %>% # Subset to VPU, only one DAMID per waterbody filter(flowlcomid %in% nhd$COMID | wbareacomi %in% wb_table$member_comid) %>% @@ -600,8 +617,9 @@ if(!"Type_WBIn" %in% names(tmp_POIs)) { tmp_POIs <- wb_in_col$POIs if(!all(is.na(wb_layers$events))) { - wb_inlet_events <- select(wb_layers$events, - -c(LevelPathI, Hydroseq, ToMeas)) + wb_inlet_events <- wb_layers$events %>% + select(COMID, REACHCODE, REACH_meas, POI_identifier, + event_type, nexus) # Write out events and outlets if(!needs_layer(temp_gpkg, split_layer)){ @@ -705,7 +723,7 @@ if(!"Type_DA" %in% names(tmp_POIs)){ This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m ```{r Elevation Break POIs} -if(!"Type_Elev" %in% names(tmp_POIs)){ +if(!"Type_elev" %in% names(tmp_POIs)){ inc_segs <- segment_increment(filter(nhd, minNet == 1), filter(st_drop_geometry(nhd), @@ -801,8 +819,8 @@ if(all(is.na(tmp_POIs$Type_Travel))) { if(nrow(tt_splits) > 0) { tmp_POIs <- POI_creation(select(tt_splits, COMID, tt_class), - filter(nhd, poi == 1), "tt") %>% - addType(., tmp_POIs, "tt", nexus = TRUE) + filter(nhd, poi == 1), "Travel") %>% + addType(., tmp_POIs, "Travel", nexus = TRUE) write_sf(tmp_POIs, temp_gpkg, nav_poi_layer) } @@ -908,86 +926,8 @@ final_POIs_prec <- mutate(final_POIs, id = row_number(), moved = NA) %>% write_sf(nav_gpkg, nav_poi_layer) collapse <- TRUE -# Load data -if(collapse){ - - # Move HUC12 to other POIs - moved_pois <- poi_move(final_POIs_prec, "Type_HUC12", poi_dar_move, - poi_distance_move, - c("Type_Gages", "Type_TE", "Type_NID")) - - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois$moved_points %>% - mutate(move_type = "huc12 to other") - } else { - final_POIs <- moved_POIs - } - - # Gages to confluences - moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, - poi_distance_move, "Type_Conf") - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois_table %>% - rbind(moved_pois$moved_points %>% - mutate(move_type = "gages to conf")) - } else { - final_POIs <- moved_POIs - } - - - # Gages to waterbody inlets - moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, - poi_distance_move, "Type_WBIn") - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois_table %>% - rbind(moved_pois$moved_points %>% - mutate(move_type = "gages to wbin")) - } else { - final_POIs <- moved_pois - } - - # Gages to waterbody outlets - moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, - poi_distance_move, "Type_WBOut") # 6 - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois_table %>% - rbind(moved_pois$moved_points %>% - mutate(move_type = "gages to wbout")) - } else { - final_POIs <- moved_pois - } - - # Streamgage to term outlet - moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, - poi_distance_move, "Type_Term") - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois_table %>% - rbind(moved_pois$moved_points %>% - mutate(move_type = "gages to term")) - } else { - final_POIs <- moved_pois - } - - # NID to waterbody outlet - moved_pois <- poi_move(final_POIs, "Type_hilarri", poi_dar_move/10, - poi_distance_move * 0.4, c("Type_WBOut", "Type_TE")) - if(!is.data.frame(moved_pois)){ - final_POIs <- moved_pois$final_pois - moved_pois_table <- moved_pois_table %>% - rbind(moved_pois$moved_points %>% - mutate(move_type = "nid to wb_out")) - } else { - final_POIs <- moved_pois - } - - write_sf(final_POIs, nav_gpkg, pois_all_layer) - write_sf(moved_pois_table, temp_gpkg, "pois_collapsed") -} +# TODO: document what this is doing! +source("R/poi_move.R") check_dups <- final_POIs %>% group_by(COMID) %>% @@ -1009,7 +949,6 @@ mapview(final_POIs, layer.name = "All POIs", col.regions = "red") + ```{r lookup} # Final POI layer POIs <- final_POIs %>% - mutate(identifier = row_number()) # Unique POI geometry diff --git a/workspace/05_hyRefactor_AK.Rmd b/workspace/03_hyRefactor_AK.Rmd similarity index 100% rename from workspace/05_hyRefactor_AK.Rmd rename to workspace/03_hyRefactor_AK.Rmd diff --git a/workspace/05_hyRefactor_HI.Rmd b/workspace/03_hyRefactor_HI.Rmd similarity index 100% rename from workspace/05_hyRefactor_HI.Rmd rename to workspace/03_hyRefactor_HI.Rmd diff --git a/workspace/05_hyRefactor_flines.Rmd b/workspace/03_hyRefactor_flines.Rmd similarity index 75% rename from workspace/05_hyRefactor_flines.Rmd rename to workspace/03_hyRefactor_flines.Rmd index 49910400ce16bab7cb0cebbceae5b06249c541f2..cabb6c3a42fccc4b00cc77f84ca78b55f1371437 100644 --- a/workspace/05_hyRefactor_flines.Rmd +++ b/workspace/03_hyRefactor_flines.Rmd @@ -25,19 +25,22 @@ See `hyRefactory_config.R` for all constants. ```{r setup} source("R/utils.R") source("R/config.R") +source("R/user_vars.R") ``` Load and filter source NHD Flowline network. Bring in POIs ```{r flowlines} nhd <- read_sf(out_refac_gpkg, nhd_flowline) -cats <- read_sf(out_refac_gpkg, nhd_catchment) +#cats <- read_sf(out_refac_gpkg, nhd_catchment) POIs_data <- read_sf(nav_gpkg, poi_data_table) POIs <- read_sf(nav_gpkg, poi_geometry_table) POI_id <- st_drop_geometry(POIs) + +events <- read_sf(nav_gpkg, event_geometry_table) ``` Read POIs and add nhd outlets. Save to a non-spatial table. @@ -71,15 +74,17 @@ nhdplus_flines <- sf::st_zm(filter(nhd, refactor == 1)) %>% st_as_sf() if(!needs_layer(nav_gpkg, event_geometry_table)){ - events <- read_sf(nav_gpkg, event_geometry_table) %>% - mutate(event_identifier = as.character(row_number())) + events <- read_sf(nav_gpkg, event_geometry_table) events_ref <- filter(events, hy_id %in% nhdplus_flines$COMID) %>% rename(COMID = hy_id) %>% - distinct(COMID, REACHCODE, REACH_meas, event_identifier) + distinct(COMID, REACHCODE, REACH_meas, + poi_id, geom) %>% + arrange(poi_id) %>% + mutate(event_identifier = row_number()) - outlets <- filter(outlets, !poi_id %in% events$poi_id) %>% - rename(COMID = hy_id) + outlets <- rename(outlets, COMID = hy_id) %>% + filter(!poi_id %in% events_ref$poi_id) } else { events_ref <- NULL } @@ -100,7 +105,8 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) { three_pass = TRUE, purge_non_dendritic = FALSE, events = events_ref, - exclude_cats = c(outlets$COMID, avoid$COMID, POI_downstream$COMID), + exclude_cats = c(outlets$COMID, + avoid$COMID, POI_downstream$COMID), warn = TRUE) write_sf(st_transform(read_sf(tf), proj_crs), out_refac_gpkg, refactored_layer) @@ -123,7 +129,7 @@ refactor_lookup <- st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)) % if(is.character(refactor_lookup$reconciled_ID)) refactor_lookup$reconciled_ID <- as.integer(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") + 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) @@ -135,49 +141,49 @@ if(needs_layer(out_refac_gpkg, outlets_layer)) { select(member_COMID = COMID, Hydroseq, event_identifier, event_REACHCODE) %>% inner_join(select(st_drop_geometry(nhd), orig_COMID = COMID, Hydroseq), by = "Hydroseq") - if(!is.null(events_ref)){ + if(nrow(events_ref) > 0){ # Subset for events refactored_events <- refactored %>% - filter(!is.na(event_REACHCODE), !is.na(event_identifier)) + filter(!is.na(event_REACHCODE), !is.na(event_identifier)) %>% + mutate(event_identifier = as.numeric(event_identifier)) - event_outlets <- events %>% + event_outlets <- events_ref %>% inner_join(st_drop_geometry(refactored_events), by = "event_identifier") %>% - select(hy_id, event_id, poi_id, member_COMID) + select(hy_id = COMID, event_id = event_identifier, poi_id, member_COMID)%>% + inner_join(select(lookup_table, -NHDPlusV2_COMID), by = "member_COMID") # subset for refactored outlets (non-events) - refactored_outlets <- filter(refactored, !member_COMID %in% event_outlets$member_COMID) + refactored_outlets <- filter(refactored, + !member_COMID %in% event_outlets$member_COMID) # get ref_COMId for other outlets outlets_ref <- outlets %>% - left_join(select(st_drop_geometry(refactored_outlets), member_COMID, orig_COMID), - by = c("COMID" = "orig_COMID")) %>% + filter(!poi_id %in% event_outlets$poi_id) %>% + inner_join(lookup_table, + by = c("COMID" = "NHDPlusV2_COMID")) %>% group_by(COMID) %>% filter(member_COMID == max(member_COMID)) %>% - select(hy_id = COMID, poi_id, member_COMID, type) + select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, type) - outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets), fill = TRUE) %>% + outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets), + fill = TRUE) %>% st_as_sf() } else { # get ref_COMId for other outlets outlets_ref_COMID <- outlets %>% - left_join(select(st_drop_geometry(refactored), member_COMID, orig_COMID), - by = c("COMID" = "orig_COMID")) %>% + inner_join(lookup_table, + by = c("COMID" = "NHDPlusV2_COMID")) %>% group_by(COMID) %>% filter(member_COMID == max(member_COMID)) %>% - select(hy_id = COMID, poi_id, member_COMID, type) + select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, type) } - - final_outlets <- outlets_ref_COMID %>% - st_as_sf() %>% - inner_join(select(lookup_table, member_COMID, reconciled_ID), - by = "member_COMID") - write_sf(final_outlets, out_refac_gpkg, outlets_layer) + write_sf(outlets_ref_COMID, out_refac_gpkg, outlets_layer) } else { - final_outlets <- read_sf(out_refac_gpkg, outlets_layer) + outlets_ref_COMID <- read_sf(out_refac_gpkg, outlets_layer) } -check_dups_poi <- final_outlets %>% +check_dups_poi <- outlets_ref_COMID %>% group_by(reconciled_ID) %>% filter(n() > 1) %>% ungroup() diff --git a/workspace/06-1_hyRefactor_cats.Rmd b/workspace/04-1_hyRefactor_cats.Rmd similarity index 98% rename from workspace/06-1_hyRefactor_cats.Rmd rename to workspace/04-1_hyRefactor_cats.Rmd index b5a0aa35de5dd9bc67202c7cfdd52a3dba145d77..e920d7a8b7af7953e907cf6cf0fcb49bf3d3fcca 100644 --- a/workspace/06-1_hyRefactor_cats.Rmd +++ b/workspace/04-1_hyRefactor_cats.Rmd @@ -22,6 +22,7 @@ See `config.R` for all constants. ```{r libs} source("R/utils.R") source("R/config.R") +source("R/user_vars_ww.R") source("R/NHD_navigate.R") source("R/user_vars.R") ``` diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/04-2_aggregate_cats.Rmd similarity index 96% rename from workspace/06-2_aggregate_cats.Rmd rename to workspace/04-2_aggregate_cats.Rmd index 32f90287b57346f7ade635cbb3476497d37e0130..bb49aa1f076d253620ad695f2bfbb06b7441571f 100644 --- a/workspace/06-2_aggregate_cats.Rmd +++ b/workspace/04-2_aggregate_cats.Rmd @@ -22,6 +22,7 @@ See `hyRefactory_config.R` for all constants. ```{r setup} source("R/utils.R") source("R/config.R") +source("R/user_vars.R") source("R/NHD_navigate.R") source("R/user_vars.R") ``` @@ -104,7 +105,9 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){ mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>% select(-toID) - outlets <- distinct(outlets) + outlets <- filter(outlets, ID %in% reconciled$ID) %>% + distinct() + agg_cats <- aggregate_catchments(flowpath = reconciled, divide = divides, outlets = outlets, @@ -140,9 +143,14 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){ Type_Con = ifelse(!is.na(toID), 1, NA)) %>% st_as_sf() - final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing, + + if(nrow(POIs_missing) > 0){ + final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing, -c(Type_Term, Type_Con))), fill = TRUE) %>% - st_as_sf() + st_as_sf() + } else { + final_POIs <- POIs_att + } unaggregated_outlets <- filter(outlets_POI, !reconciled_ID %in% final_POIs$ID) diff --git a/workspace/07-1_merge.Rmd b/workspace/05-1_merge.Rmd similarity index 100% rename from workspace/07-1_merge.Rmd rename to workspace/05-1_merge.Rmd diff --git a/workspace/07-2_NonDend.Rmd b/workspace/05-2_NonDend.Rmd similarity index 69% rename from workspace/07-2_NonDend.Rmd rename to workspace/05-2_NonDend.Rmd index 3b75f3103bc1f95bac8288a8463088fe46fde0b3..e70088705e9dc629d40968ab231ec677ea33ad74 100644 --- a/workspace/07-2_NonDend.Rmd +++ b/workspace/05-2_NonDend.Rmd @@ -28,38 +28,41 @@ debug <- FALSE # Load custom functions source("R/utils.R") source("R/config.R") +source("R/user_vars.R") source("R/non_dend.R") source("R/user_vars.R") # Read in full VPU NHD set -if (vpu == "20"){ - full_nhd <- readRDS(file.path(data_paths$islands_dir, "NHDPlusNationalData/nhdplus_flowline.rds")) %>% - st_transform(proj_crs) +if (vpu_codes == "20"){ + full_nhd <- readRDS(file.path(data_paths$islands_dir, + "NHDPlusNationalData/nhdplus_flowline.rds")) %>% + st_transform(crs) } else { - # FIXME: Using the old nhdplus attributes here because we need coastlines below. - full_nhd <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline.rds")) + full_nhd <- readRDS(data_paths$VAA_rds) } -elev <- data_paths$elev_cm[grepl(paste0("Ned", substr(vpu, 1, 2)), data_paths$elev_cm, ignore.case = TRUE)] +# Elevation Data +elev <- data_paths$elev_cm[grepl(paste0("Ned", substr(vpu_codes, 1, 2)), + data_paths$elev_cm, ignore.case = TRUE)] -if (vpu == "08"){ +if (vpu_codes == "08"){ elev$rpu_03g <- data_paths$elev_cm$rpu_03g } # HUC extraction for specific NHD+ vpus -if(vpu == "02"){ +if(vpu_codes == "02"){ grep_exp <-"^02|^04" -} else if (vpu == "08") { +} else if (vpu_codes == "08") { grep_exp <- "^03|^08" } else { - grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2)) + grep_exp <- paste0("^", substr(vpu_codes, start = 1, stop = 2)) elev <- append(elev, list(rpu_03g = data_paths$elev_cm$rpu_03g)) } cat_rpu_table <- readRDS(data_paths$fullcats_table) # get full flowline set including coastlines -full_nhd <- full_nhd %>% +vpu_nhd <- full_nhd %>% filter(grepl(paste0("^", grep_exp, ".*"), .data$VPUID)) %>% align_nhdplus_names(.) @@ -81,12 +84,11 @@ if(needs_layer(ND_gpkg, xwalk_layer)){ # Read in full NHD cats cats <- filter(cats, full_cats == 1) - cats <- sf::st_make_valid(cats) divides <- sf::st_make_valid(divides) # Intersect NHD catchments and divides with HUC12 - nhd_wbd_int <- get_HUC12_Xwalk(vpu, cats, divides, + nhd_wbd_int <- get_HUC12_Xwalk(vpu_codes, cats, divides, file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv"), file.path(data_paths$nhdplus_dir, "HUC12.rds")) @@ -96,7 +98,7 @@ if(needs_layer(ND_gpkg, xwalk_layer)){ select(-c(ACRES, HU_12_MOD)) divides <- divides %>% - left_join(xwalk_divides_wbd %>% + left_join(distinct(xwalk_divides_wbd) %>% group_by(ID) %>% filter(intArea == max(intArea)) %>% ungroup() %>% @@ -110,8 +112,7 @@ if(needs_layer(ND_gpkg, xwalk_layer)){ rm(nhd_wbd_int) cats <- cats %>% - left_join( - xwalk_nhd_wbd %>% + left_join(distinct(xwalk_nhd_wbd) %>% group_by(FEATUREID) %>% slice(which.max(intArea)) %>% ungroup() %>% @@ -146,31 +147,35 @@ if(needs_layer(ND_gpkg, divides_nd)){ divides_lu <- divides %>% left_join(distinct(select(lookup, reconciled_ID, aggregated_divide_ID)), by = c("ID" = "reconciled_ID")) %>% + filter(!is.na(ID) & !is.na(aggregated_divide_ID)) %>% rename(POI_ID = aggregated_divide_ID) %>% # attribute that hyrefactor was the step used to aggregate these catchments - mutate(aggStep = ifelse(!is.na(POI_ID), "hyRef", NA)) + mutate(aggStep = "hyRef") # Identify cats witin the full catchment set not refactored or aggregated, add # to divides data frame cats_miss <- cats %>% - left_join(select(lookup, NHDPlusV2_COMID, POI_ID = aggregated_divide_ID), + left_join(select(lookup, NHDPlusV2_COMID, member_COMID, reconciled_ID , + POI_ID = aggregated_divide_ID), by = c("FEATUREID" = "NHDPlusV2_COMID")) %>% - filter(is.na(POI_ID)) %>% - select(FEATUREID, HUC_12_int, intArea, AreaHUC12, POI_ID) %>% - mutate(ID = NA, aggStep = NA) %>% + filter(is.na(POI_ID), !member_COMID %in% divides_lu$member_COMID) %>% + mutate(aggStep = NA) %>% distinct() %>% left_join(select(cat_rpu_table, FEATUREID, rpu = RPUID)) %>% - mutate(member_COMID = as.character(FEATUREID)) %>% - select(-FEATUREID) + #filter(!reconciled_ID %in% divides_lu$ID) %>% + select(-member_COMID) %>% + mutate(areasqkm = as.numeric(st_area(.)/1000000)) %>% + select(ID = reconciled_ID, member_COMID = FEATUREID, rpu, HUC_12_int, intArea, AreaHUC12, POI_ID, aggStep) divides_lu <- divides_lu %>% select(-areasqkm) %>% rbind(cats_miss) %>% # resolve terminal non-POI outlets where we can - nhdplusTools:::check_valid() + nhdplusTools:::check_valid() %>% + st_cast("POLYGON") + # clean_geometry(divides_lu, ID = "POI_ID", keep = 1.0, crs = st_crs(divides_lu)) divides_lu <- dissolve_holes(divides_lu) - rm(cats_miss) # DEBUG: @@ -178,7 +183,7 @@ if(needs_layer(ND_gpkg, divides_nd)){ if(debug) { protoHRUs <- divides_lu %>% group_by("POI_ID") %>% - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% sfheaders::sf_remove_holes(.) %>% st_make_valid() write_sf(protoHRUs, ND_gpkg, HRU_layer) @@ -206,11 +211,11 @@ if(!"HUC_12" %in% unique(divides_lu$aggStep)){ } divides_lu <- assign_HUC12(divides_lu, xwalk_nhd_wbd, - filter(full_nhd, FTYPE != "Coastline"), CAC_thresh) + filter(vpu_nhd, FTYPE != "Coastline"), CAC_thresh) # Determine if any unaggregated catchments match HUC10 footprint (mostly in the west) divides_lu <- assign_HUC10(divides_lu, xwalk_nhd_wbd, - filter(full_nhd, FTYPE != "Coastline"), CAC_thresh) %>% + filter(vpu_nhd, FTYPE != "Coastline"), CAC_thresh) %>% select(-comid_char) if(debug) { @@ -218,7 +223,7 @@ if(!"HUC_12" %in% unique(divides_lu$aggStep)){ # HRU layer protoHRUs <- divides_lu %>% group_by(POI_ID) %>% - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% sfheaders::sf_remove_holes(.) %>% st_make_valid() write_sf(protoHRUs, ND_gpkg, HRU_layer) @@ -240,27 +245,28 @@ Note that in this case it can involve ugly multi-part HRUs which we may want to ```{r Coastal Cats} # aggregate frontal hucs of same type -if("Coastline" %in% unique(full_nhd$FTYPE)){ +if("Coastline" %in% unique(vpu_nhd$FTYPE)){ print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID")) # Only run through if coastal catchments are present if(!"coastal" %in% unique(divides_lu$aggStep)){ # Function to create coastal catchments by HUC12 and catch inflows. - divides_lu <- coastal_cats(divides_lu, full_nhd, vpu_WBD) + divides_lu <- coastal_cats(divides_lu, vpu_nhd, vpu_WBD) if(debug) { # DEBUG: # HRU layer protoHRUs <- divides_lu %>% group_by(POI_ID) %>% - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% sfheaders::sf_remove_holes(.) %>% st_make_valid() # Write out updates write_sf(protoHRUs, ND_gpkg, HRU_layer) rm(protoHRUs) } + # Update missing_cats write_sf(divides_lu, ND_gpkg, divides_nd) @@ -270,7 +276,7 @@ if("Coastline" %in% unique(full_nhd$FTYPE)){ } else { print ("No Coastlines in VPU") } -full_nhd <- filter(full_nhd, FTYPE != "Coastline") +vpu_nhd <- filter(vpu_nhd, FTYPE != "Coastline") ``` @@ -279,9 +285,8 @@ This is where alot of the more complex steps happen. For most VPUs, there are no 1, Determine where network terminals below the drainage area threshold used to derive the hyRefactor outlets contribute to (__assignable_cats__ in __non_dend.R__) 2. Assign other cats by the basis of HUC12 and HUC12 of bordering catchments (__assignable_cats__) 3. Uses the FAC grid to assign the rest ( __assign_holes__ ) - ```{r FixingtheHoles} -if(needs_layer(ND_gpkg, "missing_cats")){ +if(needs_layer(ND_gpkg, missing_cats)){ print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID")) if(sum(is.na(divides_lu$POI_ID)) > 0){ @@ -293,8 +298,7 @@ if(needs_layer(ND_gpkg, "missing_cats")){ mutate(row_ID = row_number()) %>% arrange(row_ID) - - HUC_sinks <- NHD_sinks(divides_lu, area_thresh = min_da_km/2, + HUC_sinks <- NHD_sinks(divides_lu, area_thresh = min_da_km_terminal/2, HUC12_table =file.path(data_paths$nhdplus_dir, "HUC12.rds"), NHD_sinks = read_sf(data_paths$nhdplus_gdb, "Sink")) @@ -303,8 +307,38 @@ if(needs_layer(ND_gpkg, "missing_cats")){ sinks_table <- HUC_sinks$sink_cats_table } + # Scan for terminals that may have been refactored + missing_ds <- filter(divides_lu, is.na(POI_ID)) + term_refactored <- lapply(missing_ds$member_COMID, function(x){ + ds_ref <- get_DM(vpu_nhd,x, include = F) + if(length(ds_ref) == 0){ + return(filter(missing_ds, member_COMID == x)) + } + + lookup_test <- filter(lookup, NHDPlusV2_COMID %in% ds_ref) + divides_test <- filter(divides_lu, ID %in% lookup_test$reconciled_ID) %>% + filter(!is.na(POI_ID)) + + print(unique(divides_test$POI_ID)) + + if(length(unique(divides_test$POI_ID)) == 1){ + return(filter(missing_ds, member_COMID == x) %>% + mutate(POI_ID = unique(divides_test$POI_ID))) + } else { + return(filter(missing_ds, member_COMID == x)) + } + + }) + + term_refactored <- data.table::rbindlist(term_refactored[lengths(term_refactored) > 1], + fill = TRUE) %>% + st_as_sf() + + divides_lu <- filter(divides_lu, !member_COMID %in% term_refactored$member_COMID) %>% + rbind(term_refactored) + if(sum(is.na(divides_lu$POI_ID)) > 0) { - divides_dem <- miss_term_assign(divides_lu, full_nhd, elev) + divides_dem <- miss_term_assign(divides_lu, vpu_nhd, elev) divides_lu <- divides_lu %>% left_join(select(divides_dem, member_COMID, agg_ID), @@ -314,16 +348,12 @@ if(needs_layer(ND_gpkg, "missing_cats")){ select(-agg_ID) if(exists("sinks_table")){ - sinks_table_fin <- filter(sinks_table, !member_COMID %in% divides_dem$member_COMID) %>% - rbind(divides_dem) + sinks_table_fin <- filter(sinks_table, !member_COMID %in% divides_dem$member_COMID) + sinks_table_fin <- data.table::rbindlist(list(sinks_table_fin, + divides_dem), fill = TRUE) } else { sinks_table_fin <- divides_dem } - - - if(sum(is.na(divides_lu$POI_ID)) > 0){ - write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats) - } write_sf(sinks_table_fin, ND_gpkg, ND_table) } @@ -332,10 +362,15 @@ if(needs_layer(ND_gpkg, "missing_cats")){ print ("all unaggregated catchments assigned") } + divides_lu <- dissolve_holes(divides_lu) + if(sum(is.na(divides_lu$POI_ID)) > 0){ + write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats) + } + # Prob HRU - filter(all_hrus, POI_ID == 140402000209) all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>% group_by(POI_ID) %>% - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% sfheaders::sf_remove_holes(.) %>% nhdplusTools:::check_valid(.) @@ -346,3 +381,32 @@ if(needs_layer(ND_gpkg, "missing_cats")){ all_hrus <- read_sf(ND_gpkg, HRU_layer) } ``` + +```{r Manual checks} +noagg_divides <- read_sf(ND_gpkg, missing_cats) %>% + select(row_ID, POI_ID_noagg = POI_ID) %>% + st_drop_geometry() + +if(all(!is.na(noagg_divides$POI_ID_noagg))){ + print ("all unaggregated divides accounted for") + + divides_lu <- divides_lu %>% + left_join(noagg_divides, by = "row_ID") %>% + mutate(POI_ID = ifelse(!is.na(POI_ID_noagg), POI_ID_noagg, POI_ID), + aggStep = ifelse(!is.na(aggStep), "manual", aggStep)) %>% + select(-POI_ID_noagg) + + # Prob HRU - filter(all_hrus, POI_ID == 140402000209) + all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>% + group_by(POI_ID) %>% + summarize(do_union = TRUE) %>% + sfheaders::sf_remove_holes(.) %>% + nhdplusTools:::check_valid(.) + + write_sf(divides_lu, ND_gpkg, divides_nd) + write_sf(all_hrus, ND_gpkg, HRU_layer) +} else { + print ("account for unaggregated divides") +} + +``` diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R index 7a4d3f9eeba2ce1070e23da6c68f2cff17f9fa45..3aaacd5b0eba95c911e763521fcbb0c90ef06a30 100644 --- a/workspace/R/NHD_navigate.R +++ b/workspace/R/NHD_navigate.R @@ -740,7 +740,7 @@ WB_event <- function(WBs, nhd_wb, type){ arrange(desc(Hydroseq)) %>% group_by(LevelPathI) %>% # union together - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% st_cast("LINESTRING") WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE") @@ -1288,7 +1288,7 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){ gage_reach <- gage_ds %>% group_by(REACHCODE) %>% - summarize(do_union = T, + summarize(do_union = TRUE, total_length = sum(LENGTHKM)) #print(nrow(gage_reach)) @@ -1433,7 +1433,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){ gage_reach <- gage_ds %>% group_by(REACHCODE) %>% - summarize(do_union = T, + summarize(do_union = TRUE, total_length = sum(LENGTHKM)) #print(nrow(gage_reach)) @@ -1719,4 +1719,4 @@ cs_group <- function(a, athres) { result = c(result, group) } return (result) -} \ No newline at end of file +} diff --git a/workspace/R/non_dend.R b/workspace/R/non_dend.R index 365052f38a90cfe09c41692578dfbc704c98fc88..090819d7284d014d08f4040cf0af5424040537cc 100644 --- a/workspace/R/non_dend.R +++ b/workspace/R/non_dend.R @@ -10,7 +10,8 @@ dissolve_holes <- function(divides_poi){ # HRU layer protoHRUs <- filter(divides_poi, !is.na(POI_ID)) %>% group_by(POI_ID) %>% - summarize(do_union = T) + summarize(do_union = TRUE) %>% + clean_geometry(., ID = "POI_ID") agg_noholes <- sfheaders::sf_remove_holes(protoHRUs) %>% st_cast("MULTIPOLYGON") @@ -467,7 +468,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){ filter(HUC_10 %in% divides_poi$HUC_10) %>% filter(HU_10_TYPE == "C" | HU_12_TYPE == "C") - # Determine which HUC12s have only one or less aggregated catcment + # Determine which HUC12s have only one or less aggregated catchment HUC12_poi <- filter(divides_poi, !is.na(POI_ID)) %>% distinct(POI_ID, HUC_12_int) %>% group_by(HUC_12_int) %>% @@ -501,7 +502,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){ # Aggretate sink cats to HUC12 and add area sink_cats_agg <- filter(sink_cats, is.na(POI_ID)) %>% group_by(HUC_12_int) %>% - summarize(do_union = T) %>% + summarize(do_union = TRUE) %>% mutate(agg_id = row_number()) m_per_km <- 1000 @@ -546,7 +547,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){ sink_cats_nhd <- filter(divides_poi, is.na(POI_ID)) %>% filter(grepl("-", .$member_COMID)) - agg_cats_small <- filter(agg_cats_poly, area < min_da_km) + agg_cats_small <- filter(agg_cats_poly, area < area_thresh) if (nrow(sink_cats_small) > 0){ # Get pertinent boundaries @@ -586,21 +587,27 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){ mutate(POI_ID = ifelse(!is.na(POI_ID_new), POI_ID_new, POI_ID), aggStep = ifelse(!is.na(POI_ID_new), "HUC12_sinks, small", aggStep)) %>% dplyr::select(-POI_ID_new) + + small_sinks <- select(sink_cats_nhd_table, member_COMID, HUC_12 = HUC_12_int, neighbor_COMID = bound_cat, + neighbor_HUC12 = bound_cat_HUC12, agg_ID = bound_cat_POI_ID, agg_id, aggStep) %>% + filter(member_COMID %in% sink_cats_small$member_COMID, !member_COMID %in% sink_cats_final$member_COMID) %>% + select(-agg_id) %>% + filter(HUC_12 == neighbor_HUC12) %>% + mutate(aggStep = "HUC12_sinks, small") %>% + st_drop_geometry() } } - large_sinks <- select(sink_cats, member_COMID, HUC_12 = HUC_12_int, agg_ID = POI_ID_new, aggStep) %>% + large_sinks <- select(sink_cats, member_COMID, HUC_12 = HUC_12_int, POI_ID, aggStep) %>% st_drop_geometry() - small_sinks <- select(sink_cats_nhd_table, member_COMID, HUC_12 = HUC_12_int, neighbor_COMID = bound_cat, - neighbor_HUC12 = bound_cat_HUC12, agg_ID = bound_cat_POI_ID, agg_id, aggStep) %>% - filter(member_COMID %in% sink_cats_small$member_COMID, !member_COMID %in% sink_cats_final$member_COMID) %>% - select(-agg_id) %>% - filter(HUC_12 == neighbor_HUC12) %>% - mutate(aggStep = "HUC12_sinks, small") %>% - st_drop_geometry() + if(exists("small_sinks")){ + sinks_table <- data.table::rbindlist(list(large_sinks, small_sinks), fill = TRUE) + } else { + sinks_table <- large_sinks + } + - sinks_table <- data.table::rbindlist(list(large_sinks, small_sinks), fill = TRUE) } return(list(divides_poi = divides_poi, sink_cats_table = sinks_table)) } diff --git a/workspace/R/poi_move.R b/workspace/R/poi_move.R new file mode 100644 index 0000000000000000000000000000000000000000..ebfdd7a06530d5af3f1fa65318f8be33ebbed02d --- /dev/null +++ b/workspace/R/poi_move.R @@ -0,0 +1,119 @@ +## TODO: document this file! + +# Load data +if(collapse){ + + # Move HUC12 to other POIs + moved_pois <- poi_move(final_POIs_prec, "Type_HUC12", poi_dar_move, + poi_distance_move) + + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois$moved_points %>% + mutate(move_type = "huc12 to other") + } else { + final_POIs <- moved_POIs + } + + # Gages to confluences, terminals + moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, + poi_distance_move, c("Type_Conf", "Type_Term")) + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "gages to conf")) + } else { + final_POIs <- moved_POIs + } + + # Gages to waterbody inlets + moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move, + poi_distance_move, c("Type_WBIn", "Type_WBOut")) + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "gages to wbin")) + } else { + final_POIs <- moved_pois + } + + # Waterbody inlet to confluence + moved_pois <- poi_move(final_POIs, "Type_WBIn", poi_dar_move/2, + poi_distance_move*0.4, "Type_Conf") + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "gages to term")) + } else { + final_POIs <- moved_pois + } + + # # Waterbody inlet to confluence + # moved_pois <- poi_move(final_POIs, "Type_WBOut", poi_dar_move/2, + # poi_distance_move*0.4, "Type_WBIn") + # if(!is.data.frame(moved_pois)){ + # final_POIs <- moved_pois$final_pois + # moved_pois_table <- moved_pois_table %>% + # rbind(moved_pois$moved_points %>% + # mutate(move_type = "gages to term")) + # } else { + # final_POIs <- moved_pois + # } + + # Waterbody inlet to confluence + moved_pois <- poi_move(final_POIs, c("Type_WBIn", "Type_HUC12"), poi_dar_move/2, + poi_distance_move*0.4, "Type_Conf") + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "gages to term")) + } else { + final_POIs <- moved_pois + } + + # NID to waterbody outlet + moved_pois <- poi_move(final_POIs, "Type_hilarri", poi_dar_move/2, + poi_distance_move * 0.4, c("Type_WBOut", "Type_TE")) + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "nid to wb_out")) + } else { + final_POIs <- moved_pois + } + + # NID to waterbody outlet + moved_pois <- poi_move(final_POIs, "Type_DA", poi_dar_move, + poi_distance_move) + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "nid to wb_out")) + } else { + final_POIs <- moved_pois + } + + if("Type_elev" %in% names(final_POIs)){ + # NID to waterbody outlet + moved_pois <- poi_move(final_POIs, "Type_elev", poi_dar_move, + poi_distance_move) + if(!is.data.frame(moved_pois)){ + final_POIs <- moved_pois$final_pois + moved_pois_table <- moved_pois_table %>% + rbind(moved_pois$moved_points %>% + mutate(move_type = "nid to wb_out")) + } else { + final_POIs <- moved_pois + } + } + + + write_sf(final_POIs, nav_gpkg, pois_all_layer) + write_sf(moved_pois_table, temp_gpkg, "pois_collapsed") +} \ No newline at end of file diff --git a/workspace/R/user_vars.R b/workspace/R/user_vars.R index b6488a7437f2ca05cda0511b4c8aa620edd4e847..980c4ffd0bbca918465f872daae640227d1c2abc 100644 --- a/workspace/R/user_vars.R +++ b/workspace/R/user_vars.R @@ -1,10 +1,19 @@ if(Sys.getenv("sb_user") == "CHANGEME") stop("must set sb_user") # Path mods -path_mod <- TRUE +path_mod <- FALSE if (path_mod){ + # 7-zip old_path <- Sys.getenv("PATH") Sys.setenv(PATH = paste(old_path, "C:\\Program Files\\7-zip", sep = ";")) + + # Mapshaper + old_path <- Sys.getenv("PATH") + Sys.setenv(PATH = paste(old_path, paste0(Sys.getenv("USERPROFILE"), "AppData\\Roaming\\npm"), sep = ";")) + + # Whitebox Tools + library(whitebox) + wbt_init(exe_path = paste0(Sys.getenv("USERPROFILE"), "AppData\\Roaming\\R\\data\\R\\whitebox\\WBT\\whitebox_tools.exe")) } # CRS for projected and geographic reference systems @@ -21,7 +30,8 @@ min_obs <- 90 # POI thresholds # Minimum drainage area for gages (square kilometers) -min_da_km_gages <- 5 +min_da_km_gages <- min_da_km <- 5 +sqkm_sqmi <- .386 min_da_km_terminal <- 10 # Variables used during 02_Navigate wb_area_thresh <- 1 # Waterbody size for inclusion as POI @@ -43,7 +53,7 @@ poi_distance_move <- 2.5 split_meters <- 10000 combine_meters <- 1000 reach_meas_thresh <- 10 -aggregate_da_thresh_sqkm <- 3 +aggregate_da_thresh_sqkm <- 5 # parallel factor for catchment reconciliation. para_reconcile <- 2 diff --git a/workspace/R/utils.R b/workspace/R/utils.R index 574020063d907f9cdfb984c10303b1da2520b427..7e495b2f0ae3c464ef06f750ed105876da8c7b3b 100644 --- a/workspace/R/utils.R +++ b/workspace/R/utils.R @@ -307,6 +307,10 @@ merge_refactor <- function(rpus, if("COMID" %in% names(df[[n]]) && is.numeric(df[[n]]$COMID)) df[[n]]$COMID <- as.character(df[[n]]$COMID) df[[n]] + + if("event_id" %in% names(df[[n]]) && !is.character(df[[n]]$event_id)) + df[[n]]$event_id <- as.character(df[[n]]$event_id) + df[[n]] }, n = x)) }, out = out), names(out[[1]]))