Merge branch 'main' into 'main'

Minor mods to NonDend

See merge request wma/nhgf/gfv2.0!131

workspace/07-2_NonDend.Rmd

@@ -38,6 +38,8 @@ if (vpu == "20"){
    full_nhd <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline.rds"))
 }
 
+elev <- data_paths$elev_cm[grepl(paste0("Ned", substr(vpu, 1, 2)), data_paths$elev_cm, ignore.case = TRUE)]
+
 # HUC extraction for specific NHD+ vpus
 if(vpu == "02"){
   grep_exp <-"^02|^04"
@@ -45,6 +47,7 @@ if(vpu == "02"){
   grep_exp <- "^03|^08"
 } else {
   grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
+  elev <- append(elev, list(rpu_03g = data_paths$elev_cm$rpu_03g))
 }
 
 cat_rpu_table <- readRDS(data_paths$fullcats_table)
@@ -57,8 +60,6 @@ full_nhd <- full_nhd %>%
 vpu_WBD <- readRDS(file.path(data_paths$nhdplus_dir, "HUC12.rds")) %>%
   filter(grepl(paste0("^", grep_exp, ".*"), .data$HUC_12))
 
-elev <- data_paths$elev_cm[grepl(paste0("Ned", vpu), data_paths$elev_cm, ignore.case = TRUE)]
-
 nhd <- st_transform(read_sf(gf_gpkg, nhd_flowline), crs)
 cats <- st_transform(read_sf(gf_gpkg, nhd_catchment), crs)
 divides <- st_transform(read_sf(gf_gpkg, divide_layer), crs)
@@ -98,7 +99,7 @@ if(needs_layer(ND_gpkg, xwalk_layer)){
   
   # Bring over divides/HUC12 intersection information into divides layer
   xwalk_nhd_wbd <- st_drop_geometry(nhd_wbd_int$cats_HUC12) %>%
-    select(-c(ACRES, SOURCEFC, Shape_Length, Shape_Area, HU_12_MOD))
+    select(-c(ACRES, HU_12_MOD))
 
   rm(nhd_wbd_int)
   
@@ -350,7 +351,7 @@ if(needs_layer(ND_gpkg,  "missing_cats")){
     st_make_valid()
   
   # Prob HRU - filter(all_hrus, POI_ID == 140402000209)
-  all_hrus <- divides_lu %>%
+  all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>%
     group_by(POI_ID) %>%
     summarize(do_union = T) %>%
     sfheaders::sf_remove_holes(.) %>%

workspace/R/non_dend.R

@@ -82,7 +82,7 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
     ungroup() %>%
     group_by(FEATUREID, HUC_12_int) %>%
     summarize(n_NHD = n(),
-              NHD_area = mean(AreaSqKM),
+              NHD_area = mean(AREASQKM),
               HUC_12_area = sum(AreaHUC12),
               intarea = sum(intArea),
               CAC = intarea / ((NHD_area - intarea) + (HUC_12_area - intarea) + intarea)) %>%
@@ -118,7 +118,7 @@ assign_HUC12 <- function(divides_poi, HUC12_xwalk, nhd, CAC_num){
     ungroup() %>%
     group_by(HUC_12_int) %>%
     summarize(n_HUC12 = n(),
-              NHD_area = sum(AreaSqKM),
+              NHD_area = sum(AREASQKM),
               HUC_12_area = mean(AreaHUC12),
               intarea = sum(intArea),
               CAC = intarea / ((NHD_area - intarea) + (HUC_12_area - intarea) + intarea)) %>%
@@ -233,7 +233,7 @@ assign_HUC10 <- function(divides, HUC12_xwalk, nhd, CAC_num){
     ungroup() %>%
     group_by(HUC_10) %>%
     mutate(n_NHD = n(),
-           NHD_area = mean(AreaSqKM),
+           NHD_area = mean(AREASQKM),
            HUC_10_area = mean(AreaHUC10),
            intarea = sum(intArea),
            CAC = intarea / ((NHD_area - intarea) + (HUC_10_area - intarea) + intarea)) %>%
@@ -849,35 +849,46 @@ miss_term_assign <- function(term_outlets, divides_poi, nhd, elev){
       st_intersection(divides_lu_poi) %>%
       filter(!member_COMID == incomid)
     
-    # convert cat buf to format for extract input
-    buf <- st_as_sf(cats_buff)
-    
-    #**************************************************************
-    # Then do a zonal statitics type operation of the DEM within each 
-    #      buffered polygon to get the min cell value, this is a hard coded path to where the dems reside
-    
-    elev <- as.character(elev)
-    elev <- elev[grepl(unique(missing_cats$rpu), elev)]
-    
-    dem <- terra::rast(elev)
-    ex <- terra::extract(dem, terra::vect(buf), min)
-    buf_df <- as.data.frame(ex)[,2, drop = FALSE]
-    
-    # Attribute with values
-    buf_df$outlet_COMID <- incomid
-    buf_df$neighbor_COMID <- as.character(cats_buff$member_COMID)
-    buf_df$POI_ID = as.character(cats_buff$POI_ID)
-    
-    colnames(buf_df) <- c("Elev","outlet_COMID", "neighbor_COMID", "POI_ID")
-    buf_df2 <- buf_df[which.min(buf_df$Elev), ]
-    
-    return(buf_df2)
+    if(nrow(cats_buff) == 0){
+      miss_cat <- dplyr::select(st_drop_geometry(missing_cats), 
+                         outlet_COMID = member_COMID) %>%
+        dplyr::mutate(Elev = NA, neighbor_COMID = NA, POI_ID = NA) %>%
+        dplyr::select(Elev, outlet_COMID, neighbor_COMID, POI_ID)
+      return(miss_cat)
+    } else {
+      # convert cat buf to format for extract input
+      buf <- st_as_sf(cats_buff)
+      
+      #**************************************************************
+      # Then do a zonal statitics type operation of the DEM within each 
+      #      buffered polygon to get the min cell value, this is a hard coded path to where the dems reside
+      
+      elev <- as.character(elev)
+      elev <- elev[grepl(unique(missing_cats$rpu), elev)]
+      
+      dem <- terra::rast(elev)
+      ex <- terra::extract(dem, terra::vect(buf), min)
+      buf_df <- as.data.frame(ex)[,2, drop = FALSE]
+      
+      # Attribute with values
+      buf_df$outlet_COMID <- incomid
+      buf_df$neighbor_COMID <- as.character(cats_buff$member_COMID)
+      buf_df$POI_ID = as.character(cats_buff$POI_ID)
+      
+      colnames(buf_df) <- c("Elev","outlet_COMID", "neighbor_COMID", "POI_ID")
+      buf_df2 <- buf_df[which.min(buf_df$Elev), ]
+      
+      return(buf_df2)
+    }
   }
 
-  out_df <- do.call(rbind, 
-                    pbapply::pblapply(unique(term_outlets$outlet_COMID), 
+  out_df <- do.call(rbind,
+                    pbapply::pblapply(unique(term_outlets$outlet_COMID),
                                       assign_func, divides_poi, nhd, elev, cl = NULL))
-
+  
+  # out_df <- lapply(unique(term_outlets$outlet_COMID), 
+  #                                    assign_func, divides_poi, nhd, elev)
+  
   out_df <- out_df %>%
     left_join(dplyr::select(st_drop_geometry(divides_poi), HUC12_neighbor = HUC_12_int, member_COMID), 
               by = c("neighbor_COMID" = "member_COMID")) %>%