History of nhdplus_v2
The changes to the source code of this entry in the repository, newest first, each with the differences it made to the entry's files. They are also listed on GitLab.
[entries] Lay out the spatial networks with hfd_layout() instead of cfd_layout() cc3bcb0
Changed: none
diff --git a/entries/nhdplus_v2/__init__.py b/entries/nhdplus_v2/__init__.py
index a8d3220..8a795b0 100644
--- a/entries/nhdplus_v2/__init__.py
+++ b/entries/nhdplus_v2/__init__.py
@@ -293,6 +293,6 @@ def parse(alts=None):
yield None, g
def get_pos(g, sbm_fit):
- # Spatial network: laid out with the default parameters of cfd_layout(),
+ # Spatial network: laid out with the default parameters of hfd_layout(),
# without the partition of the SBM fit
- return cfd_layout(g)
+ return hfd_layout(g)
[entries] Add the date when each entry was added, and remove icon_hash beff7c9
The dates are those of the commits that first added the entries' files (2020-08-11 for the entries of the initial commit); renamed entries keep the date of their first name. Changed: none
diff --git a/entries/nhdplus_v2/__init__.py b/entries/nhdplus_v2/__init__.py
index 5c307ed..a8d3220 100644
--- a/entries/nhdplus_v2/__init__.py
+++ b/entries/nhdplus_v2/__init__.py
@@ -54,7 +54,7 @@ citation = [('R. B. Moore, T. G. Dewald, "The Road to NHDPlus — Advancements i
'https://doi.org/10.1111/1752-1688.12389'),
('L. McKay, T. Bondelid, T. Dewald, J. Johnston, R. Moore, A. Rea, "NHDPlus Version 2: User Guide." US Environmental Protection Agency (2012)',
'https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus')]
-icon_hash = None
+date_added = '2026-10-31'
upstream_license = ("Public domain (US government work)", "https://www.usa.gov/government-works")
# The value-added attributes (with the flow table), the NHD snapshot (with the
[entries] Prefix the properties in the descriptions with (Node), (Edge) or (Graph) bc3dbbd
The lists of properties at the end of the descriptions marked edge properties with "Edge" but left node properties unmarked, or grouped them in "Nodes: `a`, ...; `b`, ..." items, so the kind of a property was often unclear. Every item now starts with (Node), (Edge) or (Graph), followed by the property names, a colon and the explanation, one item per property or group of properties; qualifications (kinds of nodes, networks) follow the names in parentheses. Properties that were only named got a short explanation, and a few items that named the wrong kind of property were corrected.
diff --git a/entries/nhdplus_v2/__init__.py b/entries/nhdplus_v2/__init__.py
index 7861151..5c307ed 100644
--- a/entries/nhdplus_v2/__init__.py
+++ b/entries/nhdplus_v2/__init__.py
@@ -26,27 +26,27 @@ Nodes are the flowlines (identified by their COMID), and a directed edge goes fr
The network has the following properties (missing values are -1, NaN or empty strings):
-- `comid`: the identifier of the flowline (COMID).
-- `name`, `gnis_id`: the name of the river or stream and its identifier in the Geographic Names Information System (GNIS).
-- `reachcode`: the 14-digit reach code of the NHD (whose first 8 digits are the hydrologic unit).
-- `ftype`, `fcode`: the type of the feature (e.g. `StreamRiver`, `ArtificialPath`, `CanalDitch`, `Connector`, `Coastline`, `Pipeline`) and its detailed code.
-- `vpu`: the vector processing unit of the flowline (e.g. `01`, `03N`, `10L`).
-- `in_network`: whether the flowline is part of the network of NHDPlus (has value-added attributes).
-- `length`: the length of the flowline, in km.
-- `stream_order`: the Strahler stream order (modified for divergences); `stream_calc`: the stream order calculated without the minor paths of divergences (0 for these).
-- `stream_level`: the level of the stream (1 for the main paths that end at the coast or the Great Lakes, increasing for tributaries; -9 or 0 for coastlines).
-- `hydroseq`: the hydrologic sequence number (sorting the flowlines from downstream to upstream); `level_path`: the hydroseq of the most downstream flowline of the main path (level path) of the flowline; `terminal_path`: the hydroseq of the terminal flowline of its drainage network.
-- `path_length`: the distance from the end of the flowline to the end of its drainage network along the main paths, in km; `arbolate_sum`: the total length of the flowlines upstream, including this one, in km.
-- `divergence`: 0 for flowlines without divergence, 1 for the main path of a divergence, 2 for minor paths (9 for most coastlines).
-- `start_flag`, `terminal_flag`: whether the flowline is a headwater, or ends the drainage network.
-- `area`, `total_area`, `divergence_routed_area`: the area of the catchment of the flowline, its total drainage area, and its drainage area with the flow split at divergences, in km$^2$.
-- `tidal`: whether the flowline is tidal.
-- `mean_flow`, `mean_velocity`: the estimated mean annual flow (gage-adjusted, in cubic feet per second) and velocity (in feet per second) of the Enhanced Unit Runoff Method (EROM, not available for the Pacific islands).
-- `slope`: the slope of the flowline (m/m); `max_elevation`, `min_elevation`: the smoothed elevations of its ends, in m.
-- `latitude`, `longitude`: the position of the downstream end of the flowline.
-- Edge `direction`: the code of the connection in the flow table (709 for the connections between flowlines, 714 for the connections into coastlines, from rivers or between consecutive coastlines).
-- Edge `delta_level`: the difference of the stream levels of the flowlines.
-- Edge `gap_distance`: the distance across a gap between the flowlines, in km (0 for most)."""
+- (Node) `comid`: the identifier of the flowline (COMID).
+- (Node) `name`, `gnis_id`: the name of the river or stream and its identifier in the Geographic Names Information System (GNIS).
+- (Node) `reachcode`: the 14-digit reach code of the NHD (whose first 8 digits are the hydrologic unit).
+- (Node) `ftype`, `fcode`: the type of the feature (e.g. `StreamRiver`, `ArtificialPath`, `CanalDitch`, `Connector`, `Coastline`, `Pipeline`) and its detailed code.
+- (Node) `vpu`: the vector processing unit of the flowline (e.g. `01`, `03N`, `10L`).
+- (Node) `in_network`: whether the flowline is part of the network of NHDPlus (has value-added attributes).
+- (Node) `length`: the length of the flowline, in km.
+- (Node) `stream_order`: the Strahler stream order (modified for divergences); `stream_calc`: the stream order calculated without the minor paths of divergences (0 for these).
+- (Node) `stream_level`: the level of the stream (1 for the main paths that end at the coast or the Great Lakes, increasing for tributaries; -9 or 0 for coastlines).
+- (Node) `hydroseq`: the hydrologic sequence number (sorting the flowlines from downstream to upstream); `level_path`: the hydroseq of the most downstream flowline of the main path (level path) of the flowline; `terminal_path`: the hydroseq of the terminal flowline of its drainage network.
+- (Node) `path_length`: the distance from the end of the flowline to the end of its drainage network along the main paths, in km; `arbolate_sum`: the total length of the flowlines upstream, including this one, in km.
+- (Node) `divergence`: 0 for flowlines without divergence, 1 for the main path of a divergence, 2 for minor paths (9 for most coastlines).
+- (Node) `start_flag`, `terminal_flag`: whether the flowline is a headwater, or ends the drainage network.
+- (Node) `area`, `total_area`, `divergence_routed_area`: the area of the catchment of the flowline, its total drainage area, and its drainage area with the flow split at divergences, in km$^2$.
+- (Node) `tidal`: whether the flowline is tidal.
+- (Node) `mean_flow`, `mean_velocity`: the estimated mean annual flow (gage-adjusted, in cubic feet per second) and velocity (in feet per second) of the Enhanced Unit Runoff Method (EROM, not available for the Pacific islands).
+- (Node) `slope`: the slope of the flowline (m/m); `max_elevation`, `min_elevation`: the smoothed elevations of its ends, in m.
+- (Node) `latitude`, `longitude`: the position of the downstream end of the flowline.
+- (Edge) `direction`: the code of the connection in the flow table (709 for the connections between flowlines, 714 for the connections into coastlines, from rivers or between consecutive coastlines).
+- (Edge) `delta_level`: the difference of the stream levels of the flowlines.
+- (Edge) `gap_distance`: the distance across a gap between the flowlines, in km (0 for most)."""
tags = ['Hydrology', 'Spatial', 'Unweighted', 'Metadata']
url = ['https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus',
'https://dmap-data-commons-ow.s3.amazonaws.com/index.html#NHDPlusV21/']
[entries] New entry: nhdplus_v2 f737b29
diff --git a/entries/nhdplus_v2/__init__.py b/entries/nhdplus_v2/__init__.py
new file mode 100644
index 0000000..7861151
--- /dev/null
+++ b/entries/nhdplus_v2/__init__.py
@@ -0,0 +1,298 @@
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2026 Tiago de Paula Peixoto <tiago@skewed.de>
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Affero General Public License along
+# with this program. If not, see <http://www.gnu.org/licenses/>.
+
+from .. import *
+import struct
+
+title = "US river network (NHDPlus version 2.1)"
+description = r"""The rivers, streams, canals and other flowlines of the United States as represented by NHDPlus version 2.1 (the National Hydrography Dataset Plus) of the US Environmental Protection Agency and the US Geological Survey, which combines the medium-resolution (1:100,000) National Hydrography Dataset (NHD) with the National Elevation Dataset and the Watershed Boundary Dataset into a network of flowlines with flow directions and value-added attributes, such as the stream order, the drainage areas and estimates of the mean annual flow. It covers the conterminous United States, Hawaii, Puerto Rico and the US Virgin Islands, and the Pacific islands of American Samoa, Guam and the Northern Mariana Islands (not Alaska), and is distributed in 26 vector processing units (VPUs), which are combined here into a single network.
+
+Nodes are the flowlines (identified by their COMID), and a directed edge goes from a flowline to each flowline downstream of it, as given by the flow table (`PlusFlow`) of NHDPlus, which is the basis of the navigation in the network (the connections across the boundaries of the processing units are included once). The network is not a forest of trees, since a flowline can have several flowlines downstream (divergences, e.g. in braided channels, deltas and canals; the main path is marked by `divergence`). Besides the flowlines of the network, the NHD contains flowlines that are not part of it (e.g. coastlines and some isolated or uninitialized features, without value-added attributes), which are included as nodes, marked by `in_network` (311,005 of the 3,035,617 flowlines, mostly canals and ditches and isolated streams; all of them but one are isolated nodes). The positions of the nodes are the downstream ends of the flowlines (whose geometries are digitized in the direction of flow), in NAD83 coordinates (which differ from WGS84 by about 1 m). The geometries, the catchments and the monthly flow estimates are not included.
+
+The network has the following properties (missing values are -1, NaN or empty strings):
+
+- `comid`: the identifier of the flowline (COMID).
+- `name`, `gnis_id`: the name of the river or stream and its identifier in the Geographic Names Information System (GNIS).
+- `reachcode`: the 14-digit reach code of the NHD (whose first 8 digits are the hydrologic unit).
+- `ftype`, `fcode`: the type of the feature (e.g. `StreamRiver`, `ArtificialPath`, `CanalDitch`, `Connector`, `Coastline`, `Pipeline`) and its detailed code.
+- `vpu`: the vector processing unit of the flowline (e.g. `01`, `03N`, `10L`).
+- `in_network`: whether the flowline is part of the network of NHDPlus (has value-added attributes).
+- `length`: the length of the flowline, in km.
+- `stream_order`: the Strahler stream order (modified for divergences); `stream_calc`: the stream order calculated without the minor paths of divergences (0 for these).
+- `stream_level`: the level of the stream (1 for the main paths that end at the coast or the Great Lakes, increasing for tributaries; -9 or 0 for coastlines).
+- `hydroseq`: the hydrologic sequence number (sorting the flowlines from downstream to upstream); `level_path`: the hydroseq of the most downstream flowline of the main path (level path) of the flowline; `terminal_path`: the hydroseq of the terminal flowline of its drainage network.
+- `path_length`: the distance from the end of the flowline to the end of its drainage network along the main paths, in km; `arbolate_sum`: the total length of the flowlines upstream, including this one, in km.
+- `divergence`: 0 for flowlines without divergence, 1 for the main path of a divergence, 2 for minor paths (9 for most coastlines).
+- `start_flag`, `terminal_flag`: whether the flowline is a headwater, or ends the drainage network.
+- `area`, `total_area`, `divergence_routed_area`: the area of the catchment of the flowline, its total drainage area, and its drainage area with the flow split at divergences, in km$^2$.
+- `tidal`: whether the flowline is tidal.
+- `mean_flow`, `mean_velocity`: the estimated mean annual flow (gage-adjusted, in cubic feet per second) and velocity (in feet per second) of the Enhanced Unit Runoff Method (EROM, not available for the Pacific islands).
+- `slope`: the slope of the flowline (m/m); `max_elevation`, `min_elevation`: the smoothed elevations of its ends, in m.
+- `latitude`, `longitude`: the position of the downstream end of the flowline.
+- Edge `direction`: the code of the connection in the flow table (709 for the connections between flowlines, 714 for the connections into coastlines, from rivers or between consecutive coastlines).
+- Edge `delta_level`: the difference of the stream levels of the flowlines.
+- Edge `gap_distance`: the distance across a gap between the flowlines, in km (0 for most)."""
+tags = ['Hydrology', 'Spatial', 'Unweighted', 'Metadata']
+url = ['https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus',
+ 'https://dmap-data-commons-ow.s3.amazonaws.com/index.html#NHDPlusV21/']
+citation = [('R. B. Moore, T. G. Dewald, "The Road to NHDPlus — Advancements in Digital Stream Networks and Associated Catchments." Journal of the American Water Resources Association 52, 890-900 (2016)',
+ 'https://doi.org/10.1111/1752-1688.12389'),
+ ('L. McKay, T. Bondelid, T. Dewald, J. Johnston, R. Moore, A. Rea, "NHDPlus Version 2: User Guide." US Environmental Protection Agency (2012)',
+ 'https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus')]
+icon_hash = None
+upstream_license = ("Public domain (US government work)", "https://www.usa.gov/government-works")
+
+# The value-added attributes (with the flow table), the NHD snapshot (with the
+# flowlines) and the EROM flow estimates of each processing unit; the file
+# names carry their version numbers. About 3 GB in all.
+upstream_prefix = 'https://dmap-data-commons-ow.s3.amazonaws.com/NHDPlusV21/Data'
+# (drainage area, VPU, directory, versions of the NHDPlusAttributes,
+# NHDSnapshot and EROMExtension archives)
+vpus = [("NE", "01", "NHDPlusNE", "09", "04", "06"),
+ ("MA", "02", "NHDPlusMA", "09", "04", "06"),
+ ("SA", "03N", "NHDPlusSA/NHDPlus03N", "07", "04", "05"),
+ ("SA", "03S", "NHDPlusSA/NHDPlus03S", "07", "06", "05"),
+ ("SA", "03W", "NHDPlusSA/NHDPlus03W", "07", "04", "05"),
+ ("GL", "04", "NHDPlusGL", "14", "08", "11"),
+ ("MS", "05", "NHDPlusMS/NHDPlus05", "09", "06", "06"),
+ ("MS", "06", "NHDPlusMS/NHDPlus06", "10", "09", "07"),
+ ("MS", "07", "NHDPlusMS/NHDPlus07", "10", "08", "07"),
+ ("MS", "08", "NHDPlusMS/NHDPlus08", "09", "07", "06"),
+ ("SR", "09", "NHDPlusSR", "07", "07", "05"),
+ ("MS", "10U", "NHDPlusMS/NHDPlus10U", "10", "07", "07"),
+ ("MS", "10L", "NHDPlusMS/NHDPlus10L", "12", "06", "08"),
+ ("MS", "11", "NHDPlusMS/NHDPlus11", "08", "06", "04"),
+ ("TX", "12", "NHDPlusTX", "09", "05", "06"),
+ ("RG", "13", "NHDPlusRG", "07", "05", "05"),
+ ("CO", "14", "NHDPlusCO/NHDPlus14", "10", "07", "05"),
+ ("CO", "15", "NHDPlusCO/NHDPlus15", "09", "04", "07"),
+ ("GB", "16", "NHDPlusGB", "06", "06", "04"),
+ ("PN", "17", "NHDPlusPN", "10", "08", "09"),
+ ("CA", "18", "NHDPlusCA", "08", "05", "07"),
+ ("HI", "20", "NHDPlusHI", "02", "02", "03"),
+ ("CI", "21", "NHDPlusCI", "02", "02", "03"),
+ ("PI", "22AS", "NHDPlusPI/NHDPlus22AS", "02", "01", None),
+ ("PI", "22GU", "NHDPlusPI/NHDPlus22GU", "02", "01", None),
+ ("PI", "22MP", "NHDPlusPI/NHDPlus22MP", "02", "01", None)]
+
+def _archives(da, vpu, d, va, vs, ve):
+ a = [f"{d}/NHDPlusV21_{da}_{vpu}_NHDPlusAttributes_{va}.7z",
+ f"{d}/NHDPlusV21_{da}_{vpu}_NHDSnapshot_{vs}.7z"]
+ if ve is not None:
+ a.append(f"{d}/NHDPlusV21_{da}_{vpu}_EROMExtension_{ve}.7z")
+ return a
+
+files = [(tuple(a for v in vpus for a in _archives(*v)), None, None)]
+
+# the members of the archives (whose paths differ in case between units)
+def _members(da, vpu, d, va, vs, ve):
+ root = f"NHDPlus{da}/NHDPlus{vpu}"
+ att, snap, erom = _archives(da, vpu, d, va, vs, ve) + [None] * (ve is None)
+ hyd = "hydrography" if vpu in ("20", "21", "22AS", "22MP") else "Hydrography"
+ m = {"flow": f"{att}:{root}/NHDPlusAttributes/PlusFlow.dbf",
+ "vaa": f"{att}:{root}/NHDPlusAttributes/PlusFlowlineVAA.dbf",
+ "elev": f"{att}:{root}/NHDPlusAttributes/{'Elevslope' if vpu == '20' else 'elevslope'}.dbf",
+ "dbf": f"{snap}:{root}/NHDSnapshot/{hyd}/NHDFlowline.dbf",
+ "shp": f"{snap}:{root}/NHDSnapshot/{hyd}/NHDFlowline.shp"}
+ if erom is not None:
+ m["erom"] = f"{erom}:{'NHDPlusCI' if vpu == '21' else root}/EROMExtension/EROM_MA0001.DBF"
+ return m
+
+members = {v[1]: _members(*v) for v in vpus}
+
+files = [(tuple(f for m in members.values() for f in m.values()), None, None)]
+
+def fetch_upstream(force=False):
+ # (the archives are checked by name, not by extracting their members)
+ return fetch_upstream_files(__name__.split(".")[-1], upstream_prefix, files,
+ force, check=False)
+
+def read_dbf(f, fields):
+ """Read the fields `fields` (a dict of name -> type, "float", "int" or
+ "string"; names in upper case) of a dBase (DBF) file, as numpy arrays
+ (strings as object arrays); missing numbers are NaN or -1, deleted
+ records are left out."""
+ head = f.read(32)
+ n, hlen, rlen = struct.unpack("<IHH", head[4:12])
+ rest = f.read(hlen - 32)
+ spec = {}
+ pos = 0
+ off = 1
+ while rest[pos] != 0x0D:
+ d = rest[pos:pos + 32]
+ spec[d[:11].split(b"\0")[0].decode().upper()] = (off, d[16])
+ off += d[16]
+ pos += 32
+ data = numpy.frombuffer(f.read(n * rlen), dtype=numpy.uint8)
+ if len(data) < n * rlen:
+ raise ValueError("truncated DBF file")
+ data = data.reshape(n, rlen)
+ keep = data[:, 0] != ord("*")
+ cols = {}
+ for k, t in fields.items():
+ off, l = spec[k]
+ x = numpy.char.strip(data[keep, off:off + l].copy().view(f"S{l}").ravel())
+ if t == "string":
+ cols[k] = numpy.array([s.decode("utf8", errors="replace") for s in x],
+ dtype=object)
+ else:
+ x = numpy.where(x == b"", b"nan", x).astype("double")
+ if t == "int":
+ x = numpy.where(numpy.isnan(x), -1, x).astype("int64")
+ cols[k] = x
+ return cols
+
+def iter_shp_last(f):
+ """Iterate over the records of a shapefile of (multi)lines, read as a
+ stream, yielding the last point of each (NaN for null shapes)."""
+ f.read(100)
+ while True:
+ h = f.read(8)
+ if len(h) < 8:
+ return
+ n, length = struct.unpack(">ii", h)
+ c = f.read(2 * length)
+ stype = struct.unpack("<i", c[:4])[0]
+ if stype == 0: # null shape
+ yield (numpy.nan, numpy.nan)
+ continue
+ if stype not in (3, 13, 23): # PolyLine, PolyLineZ, PolyLineM
+ raise ValueError(f"unexpected shape type: {stype}")
+ nparts, npoints = struct.unpack("<ii", c[36:44])
+ o = 44 + 4 * nparts + 16 * (npoints - 1)
+ yield struct.unpack("<dd", c[o:o + 16])
+
+# The attributes kept: (file, column, name, type)
+FLOWLINE = [("GNIS_ID", "gnis_id", "string"), ("GNIS_NAME", "name", "string"),
+ ("REACHCODE", "reachcode", "string"), ("FTYPE", "ftype", "string"),
+ ("FCODE", "fcode", "int"), ("LENGTHKM", "length", "float")]
+VAA = [("STREAMORDE", "stream_order", "int"), ("STREAMCALC", "stream_calc", "int"),
+ ("STREAMLEVE", "stream_level", "int"), ("HYDROSEQ", "hydroseq", "int"),
+ ("LEVELPATHI", "level_path", "int"), ("TERMINALPA", "terminal_path", "int"),
+ ("PATHLENGTH", "path_length", "float"), ("ARBOLATESU", "arbolate_sum", "float"),
+ ("DIVERGENCE", "divergence", "int"), ("STARTFLAG", "start_flag", "int"),
+ ("TERMINALFL", "terminal_flag", "int"), ("AREASQKM", "area", "float"),
+ ("TOTDASQKM", "total_area", "float"), ("DIVDASQKM", "divergence_routed_area", "float"),
+ ("TIDAL", "tidal", "int")]
+EROM = [("Q0001E", "mean_flow", "float"), ("V0001E", "mean_velocity", "float")]
+ELEV = [("SLOPE", "slope", "float"), ("MAXELEVSMO", "max_elevation", "float"),
+ ("MINELEVSMO", "min_elevation", "float")]
+
+def _join(ids, other_ids, vals, missing):
+ "The values `vals` of `other_ids` aligned with `ids` (`missing` if absent)."
+ out = numpy.full(len(ids), missing, dtype=vals.dtype)
+ if len(other_ids) == 0:
+ return out
+ order = numpy.argsort(other_ids)
+ j = numpy.searchsorted(other_ids[order], ids)
+ j = numpy.minimum(j, len(order) - 1)
+ found = other_ids[order][j] == ids
+ out[found] = vals[order][j[found]]
+ return out
+
+@cache_network()
+@coerce_props(vforce={"reachcode": "string", "gnis_id": "string"})
+@annotate()
+def parse(alts=None):
+ name = __name__.split(".")[-1]
+ vcols = defaultdict(list)
+ flows = defaultdict(list)
+ for vpu, m in members.items():
+ with open_upstream_file(name, m["dbf"], "rb") as f:
+ fl = read_dbf(f, {"COMID": "int", **{c: t for c, n, t in FLOWLINE}})
+ with open_upstream_file(name, m["shp"], "rb") as f:
+ pos = numpy.array(list(iter_shp_last(f)))
+ if len(pos) != len(fl["COMID"]):
+ raise ValueError(f"{vpu}: {len(pos)} shapes for {len(fl['COMID'])} records")
+ ids = fl["COMID"]
+ vcols["comid"].append(ids)
+ vcols["vpu"].append(numpy.array([vpu] * len(ids), dtype=object))
+ for c, n, t in FLOWLINE:
+ vcols[n].append(fl[c])
+ vcols["longitude"].append(pos[:, 0])
+ vcols["latitude"].append(pos[:, 1])
+ with open_upstream_file(name, m["vaa"], "rb") as f:
+ vaa = read_dbf(f, {"COMID": "int", **{c: t for c, n, t in VAA}})
+ vcols["in_network"].append(numpy.isin(ids, vaa["COMID"]))
+ for c, n, t in VAA:
+ vcols[n].append(_join(ids, vaa["COMID"], vaa[c], -1 if t == "int" else numpy.nan))
+ if "erom" in m:
+ with open_upstream_file(name, m["erom"], "rb") as f:
+ er = read_dbf(f, {"COMID": "int", **{c: t for c, n, t in EROM}})
+ else:
+ er = {"COMID": numpy.zeros(0, dtype="int64"),
+ **{c: numpy.zeros(0) for c, n, t in EROM}}
+ for c, n, t in EROM:
+ vcols[n].append(_join(ids, er["COMID"], er[c], numpy.nan))
+ with open_upstream_file(name, m["elev"], "rb") as f:
+ el = read_dbf(f, {"COMID": "int", **{c: t for c, n, t in ELEV}})
+ for c, n, t in ELEV:
+ x = el[c]
+ x = numpy.where(x <= -9998, numpy.nan, x) # missing values
+ if c != "SLOPE":
+ x = x / 100 # elevations in cm
+ vcols[n].append(_join(ids, el["COMID"], x, numpy.nan))
+ with open_upstream_file(name, m["flow"], "rb") as f:
+ fw = read_dbf(f, {"FROMCOMID": "int", "TOCOMID": "int",
+ "DIRECTION": "int", "DELTALEVEL": "int",
+ "GAPDISTKM": "float"})
+ for k, v in fw.items():
+ flows[k].append(v)
+ vcols = {k: numpy.concatenate(v) for k, v in vcols.items()}
+ flows = {k: numpy.concatenate(v) for k, v in flows.items()}
+ comid = vcols["comid"]
+ order = numpy.argsort(comid, kind="stable")
+ if numpy.any(numpy.diff(comid[order]) == 0):
+ raise ValueError("duplicate COMIDs")
+ # the connections between two flowlines (the others mark the starts and
+ # ends of the networks), counted once if listed by two units
+ sel = (flows["FROMCOMID"] > 0) & (flows["TOCOMID"] > 0)
+ flows = {k: v[sel] for k, v in flows.items()}
+ _, first = numpy.unique(numpy.column_stack((flows["FROMCOMID"], flows["TOCOMID"])),
+ axis=0, return_index=True)
+ first.sort()
+ flows = {k: v[first] for k, v in flows.items()}
+ def index(x):
+ j = numpy.minimum(numpy.searchsorted(comid[order], x), len(order) - 1)
+ if numpy.any(comid[order][j] != x):
+ raise ValueError(f"unknown COMIDs: {numpy.unique(x[comid[order][j] != x])[:10]}")
+ return order[j]
+ g = Graph(directed=True)
+ g.add_vertex(len(comid))
+ g.add_edge_list(numpy.column_stack((index(flows["FROMCOMID"]), index(flows["TOCOMID"]))))
+ for k, v in vcols.items():
+ if v.dtype == object:
+ t = "string"
+ elif v.dtype == bool:
+ t = "bool"
+ elif v.dtype.kind == "i":
+ t = "int64_t"
+ else:
+ t = "double"
+ g.vp[k] = g.new_vp(t, vals=v)
+ for k in ("tidal", "start_flag", "terminal_flag"):
+ g.vp[k] = g.new_vp("bool", vals=g.vp[k].a == 1)
+ g.ep.direction = g.new_ep("int", vals=flows["DIRECTION"])
+ g.ep.delta_level = g.new_ep("int", vals=flows["DELTALEVEL"])
+ g.ep.gap_distance = g.new_ep("double", vals=flows["GAPDISTKM"])
+ yield None, g
+
+def get_pos(g, sbm_fit):
+ # Spatial network: laid out with the default parameters of cfd_layout(),
+ # without the partition of the SBM fit
+ return cfd_layout(g)