nhdplusTools：get_split_catchment 未描绘倾点以上的整个流域集水盆地

Question

我正在尝试使用包

get_split_catchment()

中的函数

nhdplusTools

从倾泻点描绘分水岭集水盆地。该函数正确地描绘了与倾点关联的 COMID，但是，当它遇到下一个上游 COMID 时，它似乎停止运行。理想的输出是为所提供倾点的整个上游区域绘制一个盆地。

library(tidyverse)
library(nhdplusTools)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# 1. Define the starting point (pour point) ----
start_point <- st_sfc(st_point(c(-78.3361152, 38.9767759)), crs = 4269)

# Snap point to the nearest NHDPlus flowline
trace <- get_raindrop_trace(start_point)

snap_point <- st_sfc(st_point(trace$intersection_point[[1]]),
                         crs = 4326)

# 2. Retrieve NHDPlus COMID for the snapped point ----
start_comid <- discover_nhdplus_id(snap_point)

# 3. Navigate and extract upstream tributary flowlines ----
flowline <- navigate_nldi(list(featureSource = "comid", 
                               featureID = start_comid), 
                          mode = "upstreamTributaries", 
                          distance_km = 1000)

# 4. Download and subset NHDPlus data for the flowlines ----
subset_file <- tempfile(fileext = ".gpkg")
subset <- subset_nhdplus(comids = as.integer(flowline$UT$nhdplus_comid),
                         output_file = subset_file,
                         nhdplus_data = "download", 
                         flowline_only = FALSE,
                         return_data = TRUE, overwrite = TRUE)
#> All intersections performed in latitude/longitude.
#> Reading NHDFlowline_Network
#> Writing NHDFlowline_Network
#> Reading CatchmentSP
#> Writing CatchmentSP
flowline <- subset$NHDFlowline_Network

# Define the NLDI feature list for the starting COMID ----
nldi_nwis <- list(featureSource = "comid", featureID = start_comid)

# 5. Delineate the basin for the snap point ----
basin <- get_split_catchment(snap_point, upstream = TRUE)
basin <- basin[2,3]

# 6. Visualize watershed ----
ggplot() +
  geom_sf(data = flowline, color = "blue", size = 1) +
  geom_sf(data = basin, fill = NA, color = "red", size = 1) +
  geom_sf(data = start_point, color = "green", size = 3, shape = 21, fill = "green") +
  theme_bw() +
  theme(axis.text = element_text(size = 14, color = 'black'))

下图中，划定的盆地为红色区域，倾点为绿点，流线为蓝线。

sessionInfo()
#> R version 4.4.1 (2024-06-14 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 22631)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] sf_1.0-16          nhdplusTools_1.2.1 lubridate_1.9.3    forcats_1.0.0     
#>  [5] stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5       
#>  [9] tidyr_1.3.1        tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0   
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.5         xfun_0.45            tzdb_0.4.0          
#>  [4] vctrs_0.6.5          tools_4.4.1          generics_0.1.3      
#>  [7] curl_5.2.1           parallel_4.4.1       proxy_0.4-27        
#> [10] fansi_1.0.6          pkgconfig_2.0.3      KernSmooth_2.23-24  
#> [13] data.table_1.15.4    assertthat_0.2.1     lifecycle_1.0.4     
#> [16] farver_2.1.2         compiler_4.4.1       munsell_0.5.1       
#> [19] fst_0.9.8            htmltools_0.5.8.1    class_7.3-22        
#> [22] yaml_2.3.9           pillar_1.9.0         classInt_0.4-10     
#> [25] hydroloom_1.0.2      cachem_1.1.0         wk_0.9.2            
#> [28] zip_2.3.1            tidyselect_1.2.1     digest_0.6.36       
#> [31] stringi_1.8.4        arrow_16.1.0         fastmap_1.2.0       
#> [34] grid_4.4.1           colorspace_2.1-0     cli_3.6.3           
#> [37] magrittr_2.0.3       utf8_1.2.4           e1071_1.7-14        
#> [40] withr_3.0.0          scales_1.3.0         bit64_4.0.5         
#> [43] timechange_0.3.0     rmarkdown_2.27       httr_1.4.7          
#> [46] dataRetrieval_2.7.15 bit_4.0.5            RANN_2.6.1          
#> [49] hms_1.1.3            fstcore_0.9.18       memoise_2.0.1       
#> [52] pbapply_1.7-2        evaluate_0.24.0      knitr_1.48          
#> [55] s2_1.1.7             rlang_1.1.4          Rcpp_1.0.12         
#> [58] glue_1.7.0           DBI_1.2.3            xml2_1.3.6          
#> [61] reprex_2.1.1         rstudioapi_0.16.0    jsonlite_1.8.8      
#> [64] R6_2.5.1             fs_1.6.4             units_0.8-5

^{创建于 2024-07-31，使用 reprex v2.1.1}

Answer 1

事实证明这是一个罕见但已知的问题，可以通过稍微调整纬度/经度倾点来解决。

library(tidyverse)
library(nhdplusTools)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# 1. Define the starting point (pour point) ----
start_point <- st_sfc(st_point(c(-78.336, 38.976)), crs = 4326)

# Snap point to the nearest NHDPlus flowline
trace <- get_raindrop_trace(start_point)

# Use the intersection point as the pour point
snap_point <- st_sfc(st_point(trace$intersection_point[[1]]),
                         crs = 4326)

# 2. Retrieve NHDPlus COMID for the snapped point ----
start_comid <- discover_nhdplus_id(snap_point)

# 3. Navigate and extract upstream tributary flowlines ----
flowline <- navigate_nldi(list(featureSource = "comid", 
                               featureID = start_comid), 
                          mode = "upstreamTributaries", 
                          distance_km = 1000)

# 4. Download and subset NHDPlus data for the flowlines ----
subset_file <- tempfile(fileext = ".gpkg")
subset <- subset_nhdplus(comids = as.integer(flowline$UT$nhdplus_comid),
                         output_file = subset_file,
                         nhdplus_data = "download", 
                         flowline_only = FALSE,
                         return_data = TRUE, overwrite = TRUE)
#> All intersections performed in latitude/longitude.
#> Reading NHDFlowline_Network
#> Writing NHDFlowline_Network
#> Reading CatchmentSP
#> Writing CatchmentSP
flowline <- subset$NHDFlowline_Network

# Define the NLDI feature list for the starting COMID ----
nldi_nwis <- list(featureSource = "comid", featureID = start_comid)

# 5. Delineate the basin for the snap point ----
basin <- get_split_catchment(snap_point, upstream = TRUE)
basin <- basin[2,3]

# 6. Visualize watershed ----
ggplot() +
  geom_sf(data = flowline, color = "blue", size = 1) +
  geom_sf(data = basin, fill = NA, color = "red", size = 1) +
  geom_sf(data = start_point, color = "green", size = 3, shape = 21, fill = "green") +
  theme_bw() +
  theme(axis.text = element_text(size = 14, color = 'black'))

^{创建于 2024-07-31，使用 reprex v2.1.1}