Missing Catchments in GF_01 reference layers, odd ID assignment in refactored output

I got blocked on adding code to the original post so here's the example:

gpkg = "/Users/mjohnson/GF_01.gpkg"

comids = c(6747170, 6745232, 6747172)

(fl <- sf::read_sf(gpkg, 'nhd_flowline')  |>
  sf::st_transform(5070) |> 
  dplyr::filter(COMID %in% comids))
#> Simple feature collection with 3 features and 36 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2000352 ymin: 2428754 xmax: 2002178 ymax: 2430478
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 3 × 37
#>     COMID VPUID RPUID FTYPE      FromNode ToNode StartFlag StreamCalc Divergence
#> *   <int> <chr> <chr> <chr>         <dbl>  <dbl>     <dbl>      <dbl>      <dbl>
#> 1 6747170 01    01a   StreamRiv…   1.50e8 1.50e8         0          1          0
#> 2 6745232 01    01a   StreamRiv…   1.50e8 1.50e8         1          1          0
#> 3 6747172 01    01a   StreamRiv…   1.50e8 1.50e8         1          1          0
#> # … with 28 more variables: WBAREACOMI <int>, DnMinorHyd <dbl>, toCOMID <int>,
#> #   FCODE <int>, LENGTHKM <dbl>, REACHCODE <chr>, FromMeas <dbl>, ToMeas <dbl>,
#> #   AreaSqKM <dbl>, ArbolateSu <dbl>, TerminalPa <int>, Hydroseq <int>,
#> #   LevelPathI <int>, Pathlength <dbl>, DnLevelPat <int>, DnHydroseq <int>,
#> #   TotDASqKM <dbl>, TerminalFl <int>, streamleve <int>, StreamOrde <int>,
#> #   dend <dbl>, poi <dbl>, minNet <dbl>, WB <dbl>, struct_POI <dbl>,
#> #   struct_Net <dbl>, POI_ID <dbl>, geom <LINESTRING [m]>

(catchments = sf::read_sf(gpkg, 'nhd_catchment') |>
  sf::st_set_crs(5070) |> 
  dplyr::filter(FEATUREID %in% comids))
#> Simple feature collection with 1 feature and 6 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1999305 ymin: 2426385 xmax: 2002785 ymax: 2431215
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 1 × 7
#>   FEATUREID AREASQKM VPUID RPUID hy_cats full_cats                          geom
#> *     <dbl>    <dbl> <chr> <chr>   <dbl>     <dbl>                 <POLYGON [m]>
#> 1   6747172     7.15 01    01a         1         1 ((1999815 2427285, 1999995 2…


plot(catchments$geom, col = "gray")
plot(fl['COMID'], add = TRUE, lwd = 3)


ids = c(10023610, 10023647)

(fl <- sf::read_sf(gpkg, 'reconciled')  |>
    sf::st_transform(5070) |> 
    dplyr::filter(ID %in% ids))
#> Simple feature collection with 2 features and 7 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 2000352 ymin: 2427877 xmax: 2002471 ymax: 2430478
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 2 × 8
#>         ID     toID LENGTHKM TotDASqKM member_COMID      LevelPathID refactor_ID
#> *    <dbl>    <dbl>    <dbl>     <dbl> <chr>                   <int>       <int>
#> 1 10023610 10024735    2.27      0.713 6747170,6745232,…       33946       23611
#> 2 10023647 10023646    0.642     7.14  6747172                 22081       23648
#> # … with 1 more variable: geom <MULTILINESTRING [m]>

(catchments = sf::read_sf(gpkg, 'divides') |>
    sf::st_set_crs(5070) |> 
    dplyr::filter(ID %in% ids))
#> Simple feature collection with 2 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1999305 ymin: 2426385 xmax: 2002785 ymax: 2431215
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 2 × 5
#>         ID areasqkm member_COMID rpu                                        geom
#> *    <dbl>    <dbl> <chr>        <chr>                             <POLYGON [m]>
#> 1 10023610    0.724 5866323      01a   ((2002215 2426535, 2002095 2426685, 2001…
#> 2 10023647    7.15  6747172      01a   ((1999815 2427285, 1999995 2427285, 2000…


plot(catchments$geom, col = "gray")
plot(fl['ID'], add = TRUE, lwd = 3)

^{Created on 2022-05-20 by the reprex package (v2.0.1)}

changed the description

I am thinking its because its digitized the wrong way?


outlets = fl |> 
  dplyr::mutate(geom = nhdplusTools::get_node(fl, "end")$geometry) 

plot(catchments$geom, col = "gray")
plot(fl['ID'], add = TRUE, lwd = 3)
plot(outlets$geom, add = TRUE, lwd = 3)

^{Created on 2022-05-20 by the reprex package (v2.0.1)}

This one's interesting.

6747170 and 6745232 are both marked Catchment: N in the burn line event table. This indicates that those lines will not have catchments created for them.

So, from a catchment perspective, they don't exist and no land-area will ever route to them.

If you search the divides for a member_COMID that includes those two orphan flow lines, you won't find any.

So, should we just drop little orphan networks like these when we start the process? They really aren't adding anything.

Looking more closely at the flow direction issue -- this can be fixed. This is a case where the flowline is being replaced by a burn line but, because the no-catchment flowlines connect to it, the direction is not being checked. Definitely a bug. I'll get this taken care of with an update to the enhd attributes. Going to have some fairly far-reaching impacts, but I'm glad we caught this.

OK -- I now see:

> c(6747170, 6745232, 6747172) %in% new_atts$comid
[1] FALSE FALSE  TRUE

Running the reference flowlines gpkg now.

I was wrong about the solution here. There's a much lighter hand approach. I'm now just culling flowlines that don't need to be included in prep nhd prior to generating the regional reference gpkgs. Merge request incoming and I'll work on getting the regional files updated tomorrow.

mentioned in commit 79fbc74b

closed with commit 613a07f0

mentioned in commit 9052cbae

mentioned in commit 12b2dac7

Missing Catchments in GF_01 reference layers, odd ID assignment in refactored output

Designs

Child items ...

Activity