Skip to contents

Network Manipulation

Mike Johnson

Lynker, NOAA-Affiliate

David Blodgett

USGS WMA
Source: vignettes/03-processing-deep-dive.Rmd
03-processing-deep-dive.Rmd

Network Manipulation: Geoprocessing

Getting the reference fabric

All reference and (precomputed) refactored data products live on ScienceBase. They can be accessed with the web interface or can be downloaded programatically. The hydrofab::get_hydrofabric() utility will download the most current geofabric for a Vector Processing Unit (VPU). Options include downloading the “refactored” (default) or “reference” data. If the requested file already exists, the file path will be returned.

Example

  • To exemplify this process we do the following:
    1. Define a terminal location of interest
    2. Subset the upstream reference fabric using hydrofab::subset_reference 
library(hydrofabric)

# Define starting feature by source and ID
(nldi_feature = list(featureSource = "nwis", featureID = "06752260"))

# Use this feature to build a reference subset
  # This function uses the NLDI to identify an upstream tributary
  # and SQL to extract a subset GPKG
subset_reference(nldi_feature = nldi_feature, 
                 gpkg = get_hydrofabric(VPU = "10L",  type = "reference", dir = "cihro-data"),
                 export_gpkg = "cihro-data/reference.gpkg")

1) Running the Process: Refactoring

Concept

  • Refactoring is a geoprocessing workflow that seeks to

    1. Split large or long catchments into a more uniform catchment size distribution and
    2. collapse catchment topology to eliminate small catchments

  • The fundamental key is that no network resolution is lost! That means the total path length of the network going in, is what comes out.

  • The workflow is can be parameterized using three primary values:

Parameter Purpose Elected Value
split_flines_meters the maximum length flowpath desired in the output. 10,000
collapse_flines_meters the minimum length of inter-confluence flowpath desired in the output. 1,000
collapse_flines_main_meters the minimum length of between-confluence flowpaths. 1,000

You might also have areas where you want to avoid, or, enforce a splitting event. These can be defined with the following values:

Parameter Purpose
exclude_cats integer vector of COMIDs to be excluded from collapse modifications.
events data.frame containing events as generated by get_flowline_index()

  • POI definition and selection is model application specific. Here, we will ignore this aspect.

  • With only information about the network, refactor can only refactor the flowpath network.

  • In order to reconcile the catchment network,a set of flow accumulation (FAC) and flow direction (FDR) grids must be provided.

  • For the reference fabric (e.g. NHDPlusV2), we supply a national VRT for each of these that can be accessed at: s3://nextgen-hydrofabric/DEM-products/{product}.vrt

  • These (and other gridded products) can be found here

Example 

refactor = refactor(gpkg = "cihro-data/reference.gpkg",
                    split_flines_meters = 10000, 
                    collapse_flines_meters = 1000, 
                    collapse_flines_main_meters = 1000,
                    fac = '/vsis3/nextgen-hydrofabric/gridded-resources/fac.vrt',
                    fdr = '/vsis3/nextgen-hydrofabric/gridded-resources/fdr.vrt',
                    outfile = "cihro-data/refactor.gpkg")

Outputs

To get a high level understanding of what happened with this “refactor”, we can look at the length distributions:

And the area distributions:

Lastly, we look at the feature count of the network:

  • Finally, we can zoom into a layer of this network to see what changes exist.

  • In the figure below, the white edges represent the reference catchment network, while the black edges represent the refactored network

  • Since refactoring requires the preservation of the flowpath network, the blue lines are representative of both the reference and refactored network with the caveat they are broken and different places.

2. Running the Process: Aggregating

  • Aggregation is a primarily a divide oriented workflow. It collapses the network to provide a new discretization.

  • Two aggregation methods:

    • To POIs - you define network outlets (NHM and SPARROW)
    • To a statistical distribution, with or w/o enforced POIs (NextGen).

Parameter Purpose Elected Value
ideal_size_sqkm the maximum length flowpath desired in the output. 10
min_length_km the minimum length of inter-confluence flowpath desired in the output. 1
min_area_sqkm the minimum length of between-confluence flowpaths. 3

Here, a hydrolocation POINT layer can be passed to help direct the aggregation, but for simplicity is ignored here:

Example 

aggregate_to_distribution(gpkg = "cihro-data/refactor.gpkg",
                          ideal_size_sqkm = 10, 
                          min_length_km = 1, 
                          min_area_sqkm = 3, 
                          outfile = "cihro-data/aggregate.gpkg")

Outputs

To get a high level understanding of what happens with this “refactor”, we can look at the length distributions:

And the area distributions:

  • Finally, we can zoom into a layer of this network to see what changes exist.

  • In the left-hand figure below, the white edges represent the reference catchment network, the black edges represent the refactored network, and the red represnet the aggregated network

  • In the right-hand figure, we can see the blue flowlines (from reference and refactor) that were prunned in the aggregation process.

Lastly, we can look at the cumulative network traits of each fabric:

The refactored output is shared under the Refactored Parent item of the above ScienceBase resource (available here). It is also avaliable in the Lynker s3 account here

3) NextGen

Once we have a network aggregated to a scale matching the desired hydrologic processes we need to turn it into something NextGen can use (modeling task)

Divergent Topology

Nextgen operates on a [flowpath --> nexus] vs [flowpath --> flowpath] topology

This is due to the HY Features conceptula catchment that has 1 inflow draining to one outflow.

And example of this can be seen below:

captioncaption

caption

Character based indentification

Nextgen requires integer based identification, like described in the data model, prefixed with a character string defining what the feature is

Right now, the following prefixs are used to distiguish between type of model features.

Parameter Purpose Elected Value
nexus_prefix the maximum length flowpath desired in the output. “nex-”
terminal_nexus_prefix the minimum length of inter-confluence flowpath desired in the output. “tnx-”
coastal_nexus_prefix the minimum length of between-confluence flowpaths. “cnx-”
internal_nexus_prefix the maximum length flowpath desired in the output. “inx-”
catchment_prefix the minimum length of inter-confluence flowpath desired in the output. “cat-”
waterbody_prefix the minimum length of between-confluence flowpaths. “wb-”

The following function (1) identifies nexus locations, (2) moves them when needed and (3) applies the above schema to the features.

apply_nexus_topology(gpkg = "cihro-data/aggregate.gpkg", export_gpkg = "cihro-data/nextgen.gpkg")

We can see the results of this by opening the hydrofabric and adding it to a map!

nexus = read_sf("cihro-data/nextgen.gpkg", "nexus") 

read_hydrofabric("cihro-data/nextgen.gpkg") %>% 
  mapview() + 
  nexus

Extending NWM attributes

The core utilities provide a series of flowpath, divide, and nexuses.

sf::st_layers("cihro-data/nextgen.gpkg")
#> Driver: GPKG 
#> Available layers:
#>            layer_name geometry_type features fields             crs_name
#> 1             divides       Polygon      832      9 NAD83 / Conus Albers
#> 2               nexus         Point      411      4 NAD83 / Conus Albers
#> 3           flowpaths   Line String      832     10 NAD83 / Conus Albers
#> 4             network            NA     1540     15                 <NA>
#> 5 flowpath_attributes            NA      832     16                 <NA>

However, other information is needed to run some/all NextGen formulations. These include the following:

Lake Attributes

  • WBOut Hydrolocations are mapped to the NHDPlusWaterBody COMIDs used in the NWM.

Flowpath Attributes

  • Flowpath attributes are extracted from the a Routelink file
  • The values are length averaged by the portion of length each makes up in the refactored/aggregated network

For example, if a 75m flowline has a roughness of 0.05 and a 25m flowline with a roughness of 0.2

(n = (.75 * .05) + (.25 * .2))
#> [1] 0.0875

Flowpaths attributes and lake parameters can be added by pointing to a set of NWM domain files like those found here

add_flowpath_attributes("cihro-data/nextgen.gpkg", 
                        rl_path = "/Volumes/MyBook/nextgen/RouteLink_nwm_v2_2_3.fst") %>%
   add_lake_attributes(lake_path = '/Volumes/Transcend/nwmCONUS-v216/LAKEPARM_CONUS.nc')

CFE/NOAH-OWP attributes

Catchment attribute data is needed to run some formulations like CFE and NOAH-OWP. These can be computed using the following logic.

 add_cfe_noahowp_attributes(gpkg        = 'cihro-data/nextgen.gpkg"',
                            nwm_dir     = "/Volumes/Transcend/nwmCONUS-v216/",
                            outfile     = 'cihro-data/cfe_noahowp.parquet')

By default, this code take the gridded data supplied by the soilproperties_CONUS_FullRouting.nc, wrfinput_CONUS, and GWBUCKPARM_CONUS_FullRouting.nc and computes the following summaries:

Tables Description Layer(s) Summary Function Source
bexp Beta Parameter 4 mode soilproperties_CONUS_FullRouting.nc
IVGTYP Dominant category 1 mode wrfinput_CONUS.nc
ISLTYP Dominant category 1 mode wrfinput_CONUS.nc
dksat Saturated Soil Connectivity 4 geometric mean soilproperties_CONUS_FullRouting.nc
psisat Saturated soil matric potential 4 geometric mean soilproperties_CONUS_FullRouting.nc
slope Slope Index 1 mean soilproperties_CONUS_FullRouting.nc
smcmax Saturated value of soil moisture [volumetric] 4 mean soilproperties_CONUS_FullRouting.nc
smcwlt Wilting point soil moisture [volumetric] 4 mean soilproperties_CONUS_FullRouting.nc
refkdt Parameter in the surface runoff parameterization 1 mean soilproperties_CONUS_FullRouting.nc
cwpvt Empirical canopy wind parameter 1 mean soilproperties_CONUS_FullRouting.nc
vcmx25 Maximum rate of carboxylation at 25 C [ umol CO2/m2/s] 1 mean soilproperties_CONUS_FullRouting.nc
mp Slope of Conductance to photosynthesis relationship 1 mean soilproperties_CONUS_FullRouting.nc
mfsno Snowmelt m parameter 1 mean soilproperties_CONUS_FullRouting.nc
Coef Coefficient 1 mean GWBUCKPARM_CONUS_FullRouting.nc
Zmax Zmax 1 mean GWBUCKPARM_CONUS_FullRouting.nc
Expon Exponent 1 mode GWBUCKPARM_CONUS_FullRouting.nc
Quartz Soil Quartz 4 mean soilproperties_CONUS_FullRouting.nc

Themeing

As an option, QGIS QML theming files can be added to a gpkg. Some default themeing files come with hydrofabric and can be specified/added with the append_styles utility on the desired layer_names

append_style("cihro-data/nextgen.gpkg", layer_names = c("nexus", "hydrolocations", "flowpaths", "divides", "lakes"))

In QGIS, double clicking the gpkg file will allow you to select which layers to load.

Full Run-through 


reference_gpkg = get_hydrofabric(VPU = "10L",  type = "reference", dir = "cihro-data")

# Define starting feature by source and ID
list(featureSource = "nwis", featureID = "06752260") %>% 
  # Subset Reference
  subset_reference(gpkg = reference_gpkg,
                   export_gpkg = "cihro-data/reference.gpkg")


  refactor(gpkg = "cihro-data/reference.gpkg",
           split_flines_meters = 10000, 
           collapse_flines_meters = 1000, 
           collapse_flines_main_meters = 1000,
           fac = '/vsis3/nextgen-hydrofabric/gridded-resources/fac.vrt',
           fdr = '/vsis3/nextgen-hydrofabric/gridded-resources/fdr.vrt',
           outfile = "cihro-data/refactor.gpkg") %>% 
  aggregate_to_distribution(ideal_size_sqkm = 10, 
                          min_length_km = 1, 
                          min_area_sqkm = 3, 
                          outfile = "cihro-data/aggregate.gpkg") %>% 
  apply_nexus_topology(export_gpkg = "cihro-data/nextgen.gpkg") %>% 
  add_flowpath_attributes(rl_path = "/Volumes/MyBook/nextgen/RouteLink_nwm_v2_2_3.fst") %>%
  add_lake_attributes(lake_path = '/Volumes/Transcend/nwmCONUS-v216/LAKEPARM_CONUS.nc') %>%
  append_style(layer_names = c("nexus", "hydrolocations", "flowpaths", "divides", "lakes"))

add_cfe_noahowp_attributes(gpkg        = 'cihro-data/nextgen.gpkg"',
                           nwm_dir     = "nwm-data",
                           outfile     = 'cihro-data/cfe_noahowp.parquet')

nhdplusTools:

With all network manipulations, fundamental network traits change. This requires the utilities to rapidly and efficiently recopute key network metric. The nhdplusTools package provides the option to regenerate all or some of these on the fly using graph algroithms and logic.

We will return to the get_sorted() utility in the subsetting section

nhdplusTools::add_plus_network_attributes()
nhdplusTools::get_streamorder()
nhdplusTools::calculate_total_drainage_area()
nhdplusTools::get_sorted()