Add column that indicates if a flowline has a catchment polygon

[jds485]

I recently wanted to plot all the NHDPlus v2.1 flowlines that do not have catchment polygons, and needed to load in the flowline and catchment datasets to do that. I think it would be helpful to create a new column in the flowline data table with this information (no_catchment below). You could also list the polygon ID number instead of a binary TRUE/FALSE.

#flowlines
flines <- sf::st_read('NHDPlusV21_NationalData_Seamless_Geodatabase_Lower48_07/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb', 
                      layer = 'NHDFlowline_Network')

#catchments
catchments <- sf::st_read('NHDPlusV21_NationalData_Seamless_Geodatabase_Lower48_07/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb', 
                       layer = 'Catchment')

flines$no_catchment <- TRUE
flines[which(flines$COMID %in% catchments$FEATUREID), 'no_catchment'] <- FALSE

[dblodgett-usgs]

I've always used the drainage area attribute as a proxy for "doesn't have a catchment" -- did you check if that has a 1:1 correspondence with this?

[jds485]

Thanks, Dave. The only flowlines with AreaSqKM values that are 0 do not have a catchment polygon. There are very small AreaSqKM values for some flowlines that do have catchments (0.0009 km^2) but no 0s.

> head(sort(hascat))
[1] 9e-04 9e-04 9e-04 9e-04 9e-04 9e-04
> unique(nocat)
[1] 0

[dblodgett-usgs]

OK -- that was what I remembered, but I've stumbled on this and I know a lot of people do. Any ideas for how to make it more obvious that you can use that rule to determine which flowlines have catchment polygons?

[jds485]

Maybe change the AreaSqKM to NA for flowlines that do not have catchments, and then document in metadata what NA indicates. I think that would be better than saying 0 because we know the area is larger than 0, it just has not been estimated. I'd be interested in @lkoenig-usgs thoughts on this, too

[lkoenig-usgs]

I personally would not be in favor of swapping 0 for NA in the AreaSqKM attribute. My reasoning is that AreaSqKM = 0 is set in the nhdplus data model, and it would be confusing if {nhdplusTools} (the utility/toolkit) deviated or adopted a different schema from nhdplus (the source dataset).

@jds485 Would you find it helpful if this behavior is documented, either in a vignette or a note in the function documentation for get_nhdplus (or somewhere else)?

[jds485]

Okay, sounds good. Function documentation works for me. If the NHD network flowlines or catchments change over time, then having a function that tests if all flowlines with 0 for AreaSqKM do not have catchments, and all flowlines with non-zero AreaSqKM do have catchments would be helpful to support the documentation