Handle NLDI gages better now

workspace/01_Gage_Selection.Rmd

@@ -48,13 +48,20 @@ if(file.exists("cache/dropped_gages.csv")){
   file.remove("cache/dropped_gages.csv")
 }
 
-# Update location index information from other sources where necessary from NLDI reference gage set
-ref_gages <- jsonlite::fromJSON("https://github.com/internetofwater/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson", 
-                                flatten = TRUE)$features %>%
-  select(-geometry.type, -geometry.coordinates, ref_site_no = properties.id, 
-         ref_RC = properties.nhdpv2_REACHCODE, ref_RM = properties.nhdpv2_REACH_measure, 
-         ref_COM = properties.nhdpv2_COMID) %>%
-  filter(!is.na(ref_COM) & ref_COM > 0)
+if(!file.exists("data/ref_gages.gpkg")){
+  # Update location index information from other sources where necessary from NLDI reference gage set
+  ref_gages <- jsonlite::fromJSON("https://github.com/internetofwater/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson", 
+                                  flatten = TRUE)$features %>%
+    select(-geometry.type, -geometry.coordinates, ref_site_no = properties.id, 
+           ref_RC = properties.nhdpv2_REACHCODE, ref_RM = properties.nhdpv2_REACH_measure, 
+           ref_COM = properties.nhdpv2_COMID) %>%
+    filter(!is.na(ref_COM) & ref_COM > 0)
+  
+  write_sf(ref_gages, "data/ref_gages.gpkg")
+} else {
+  ref_gages <- read_sf("data/ref_gages.gpkg")
+}
+
 ```
 
 ```{r GAGESIII-SWIM}
@@ -66,7 +73,7 @@ gagesIII_db_full <- read.dbf(file = file.path(data_paths$data_dir, "gages3GFinfo
 gagesIII_sf_full <- read_sf(data_paths$SWIM_points_path, "SWIM_gage_loc") %>%
   st_zm(drop = TRUE) %>%
   select(site_no = Gage_no, comid = COMID, reachcode = REACHCODE, reach_meas = REACH_meas) %>%
-  left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+  left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
   # Update index information with NLDI gages
   mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
   select(-c(ref_COM, ref_RC, ref_RM))
@@ -166,7 +173,7 @@ gageloc_sf_full <- read_sf(data_paths$nhdplus_gdb, "Gage") %>%
   # Disregard locations already removed in GAGESIII
   filter(!site_no %in% gagesIII_sf_full$site_no) %>%
   # Update index information with NLDI gages
-  left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+  left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
   mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
   select(-c(ref_COM, ref_RC, ref_RM))
 
@@ -269,7 +276,7 @@ gfv11_sf <- gfv11_sf_full %>%
          reach_meas = ifelse(is.na(reach_meas), FromMeas, reach_meas)) %>%
   select(-c(NHDPlusID, REACHCODE, FromMeas)) %>%
   # Update index information with NLDI gages
-  left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+  left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
   mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
   select(-c(ref_COM, ref_RC, ref_RM)) %>%
   st_sf()