Merge remote-tracking branch 'usepa/main' into demands-in-gis

.github/workflows/build_deploy_pages.yml

@@ -1,11 +1,12 @@
-# This is a basic workflow to help you get started with Actions
-
 name: docs
 
 on:
   push:
     branches: [main]
+  pull_request:
+    branches: [main]
   workflow_dispatch:
+
 jobs:
   build:
     name: Build the documentation with Sphinx
@@ -36,6 +37,7 @@ jobs:
   deploy:
     name: Deploy documentation to GitHub Pages
     needs: build
+    if: github.event_name == 'push'
     permissions:
       contents: read
       pages: write      # to deploy to Pages

.github/workflows/quick_check.yml

@@ -16,6 +16,7 @@ jobs:
     strategy:
       matrix:
         python-version: ['3.9', '3.11']
+      fail-fast: false
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v2
@@ -33,7 +34,8 @@ jobs:
         python -m pip install -e .
     - name: Run tests and coverage (unittests plus doctests)
       run: |
-        coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" -m pytest -m "not time_consuming" --doctest-modules --doctest-glob="*.rst" wntr        
+        coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" -m pytest -m "not time_consuming" --doctest-modules --doctest-glob="*.rst" wntr
+        coverage run --source=wntr --omit="*/tests/*","*/sim/network_isolation/network_isolation.py","*/sim/aml/evaluator.py" --append -m pytest --doctest-glob="*.rst" documentation        
         coverage report --fail-under=70        
     # coverage run --source=wntr --omit="*/tests/*" --append -m pytest --doctest-glob="*.rst" documentation 
 

.gitignore

@@ -20,6 +20,7 @@ temp*
 
 examples/*.inp
 wntr/tests/*.png
+wntr/tests/*.tif
 
 documentation/_local
 documentation/apidoc

documentation/gis.rst

@@ -31,6 +31,12 @@
     >>> hydrant_data = gpd.read_file(examples_dir+'/data/Net1_hydrant_data.geojson')
     >>> valve_data = gpd.read_file(examples_dir+'/data/Net1_valve_data.geojson')
 
+.. doctest::
+    :hide:
+    :skipif: gpd is None or rasterio is None
+
+    >>> elevation_data_path = examples_dir+'/data/Net1_elevation_data.tif'
+
 .. _geospatial:
 
 Geospatial capabilities
@@ -47,7 +53,8 @@ The following section describes capabilities in WTNR that use GeoPandas GeoDataF
 
 .. note:: 
    Functions that use GeoDataFrames require the Python package **geopandas** :cite:p:`jvfm21` 
-   and **rtree** :cite:p:`rtree`. Both are optional dependencies of WNTR.
+   and **rtree** :cite:p:`rtree`, and functions that use raster files require the
+   Python package **rasterio**. All three are optional dependencies of WNTR.
    Note that **shapely** is installed with geopandas.
 
 The following examples use a water network generated from Net1.inp.
@@ -112,13 +119,13 @@ For example, the junctions GeoDataFrame contains the following information:
     :skipif: gpd is None
 
     >>> print(wn_gis.junctions.head())
-         node_type  base_demand demand_pattern  elevation  initial_quality                   geometry
-    name                                                                                             
-    10    Junction        0.000              1    216.408        5.000e-04  POINT (20.00000 70.00000)
-    11    Junction        0.009              1    216.408        5.000e-04  POINT (30.00000 70.00000)
-    12    Junction        0.009              1    213.360        5.000e-04  POINT (50.00000 70.00000)
-    13    Junction        0.006              1    211.836        5.000e-04  POINT (70.00000 70.00000)
-    21    Junction        0.009              1    213.360        5.000e-04  POINT (30.00000 40.00000)
+          elevation  initial_quality                   geometry
+    name                                                       
+    10      216.408        5.000e-04  POINT (20.00000 70.00000)
+    11      216.408        5.000e-04  POINT (30.00000 70.00000)
+    12      213.360        5.000e-04  POINT (50.00000 70.00000)
+    13      211.836        5.000e-04  POINT (70.00000 70.00000)
+    21      213.360        5.000e-04  POINT (30.00000 40.00000)
 
 Each GeoDataFrame contains attributes and geometry:
 
@@ -334,23 +341,23 @@ and then translates the GeoDataFrames coordinates to EPSG:3857.
 
     >>> wn_gis = wntr.network.to_gis(wn, crs='EPSG:4326')
     >>> print(wn_gis.junctions.head())
-         node_type  base_demand demand_pattern  elevation  initial_quality                   geometry
-    name                                                                                             
-    10    Junction        0.000              1    216.408        5.000e-04  POINT (20.00000 70.00000)
-    11    Junction        0.009              1    216.408        5.000e-04  POINT (30.00000 70.00000)
-    12    Junction        0.009              1    213.360        5.000e-04  POINT (50.00000 70.00000)
-    13    Junction        0.006              1    211.836        5.000e-04  POINT (70.00000 70.00000)
-    21    Junction        0.009              1    213.360        5.000e-04  POINT (30.00000 40.00000)
+          elevation  initial_quality                   geometry
+    name                                                       
+    10      216.408        5.000e-04  POINT (20.00000 70.00000)
+    11      216.408        5.000e-04  POINT (30.00000 70.00000)
+    12      213.360        5.000e-04  POINT (50.00000 70.00000)
+    13      211.836        5.000e-04  POINT (70.00000 70.00000)
+    21      213.360        5.000e-04  POINT (30.00000 40.00000)
 
     >>> wn_gis.to_crs('EPSG:3857')
     >>> print(wn_gis.junctions.head())
-         node_type  base_demand demand_pattern  elevation  initial_quality                          geometry
-    name                                                                                                    
-    10    Junction        0.000              1    216.408        5.000e-04  POINT (2226389.816 11068715.659)
-    11    Junction        0.009              1    216.408        5.000e-04  POINT (3339584.724 11068715.659)
-    12    Junction        0.009              1    213.360        5.000e-04  POINT (5565974.540 11068715.659)
-    13    Junction        0.006              1    211.836        5.000e-04  POINT (7792364.356 11068715.659)
-    21    Junction        0.009              1    213.360        5.000e-04   POINT (3339584.724 4865942.280)
+          elevation  initial_quality                          geometry
+    name                                                              
+    10      216.408        5.000e-04  POINT (2226389.816 11068715.659)
+    11      216.408        5.000e-04  POINT (3339584.724 11068715.659)
+    12      213.360        5.000e-04  POINT (5565974.540 11068715.659)
+    13      211.836        5.000e-04  POINT (7792364.356 11068715.659)
+    21      213.360        5.000e-04   POINT (3339584.724 4865942.280)
 
 Snap point geometries to the nearest point or line
 ----------------------------------------------------
@@ -822,3 +829,101 @@ the census tracts (polygons) is different than the junction and pipe attributes.
    :alt: Intersection of junctions and pipes with mean income demographic data in EPANET example Net1
 
    Net1 with mean income demographic data intersected with junctions and pipes.
+
+Sample raster at points geometries
+--------------------------------------
+
+The :class:`~wntr.gis.sample_raster` function can be used to sample a raster file at point geometries,
+such as the nodes of a water network. A common use case for this function is to assign elevation to the 
+nodes of a water network model, however other geospatial information such as climate or hazard data could be sampled 
+using this function.
+
+The network file, Net1.inp, in EPSG:4326 CRS is used in the example below. 
+The raster data in the GeoTIFF format is also in EPSG:4326 CRS.
+See :ref:`crs` for more information.
+
+.. doctest::
+    :skipif: gpd is None
+
+    >>> wn = wntr.network.WaterNetworkModel('networks/Net1.inp') # doctest: +SKIP
+    >>> wn_gis = wntr.network.to_gis(wn, crs='EPSG:4326')
+
+Sample elevations at junctions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Elevation is an essential attribute for accurate simulation of pressure in a water network and is
+commonly provided in GeoTIFF (.tif) files. The following example shows how such files can be sampled 
+and assigned to the junctions and tanks of a network. Note that elevation data generally needs 
+to be adjusted to account for buried pipes.
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> elevation_data_path = 'data/Net1_elevation_data.tif' # doctest: +SKIP
+    >>> junctions = wn_gis.junctions
+    >>> junction_elevations = wntr.gis.sample_raster(junctions, elevation_data_path)
+    >>> print(junction_elevations)
+    name
+    10    1400.0
+    11    2100.0
+    12    3500.0
+    13    4900.0
+    21    1200.0
+    22    2000.0
+    23    2800.0
+    31     300.0
+    32     500.0
+    dtype: float64
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> tanks = wn_gis.tanks
+    >>> tank_elevations = wntr.gis.sample_raster(tanks, elevation_data_path)
+    >>> print(tank_elevations)
+    name
+    2    4500.0
+    dtype: float64
+
+To use these elevations for hydraulic simulations, 
+they need to be added to the water network object.
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> for junction_name in wn.junction_name_list:
+    ...     junction = wn.get_node(junction_name)
+    ...     junction.elevation = junction_elevations[junction_name]
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> for tank_name in wn.tank_name_list:
+    ...     tank = wn.get_node(tank_name)
+    ...     tank.elevation = tank_elevations[tank_name]
+
+The sampled elevations can be plotted as follows. The 
+resulting :numref:`fig-sample-elevations` illustrates Net1 with the elevations 
+sampled from the raster file.
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> ax = wntr.graphics.plot_network(wn, node_attribute="elevation", link_width=1.5, 
+    ...     node_size=40, node_colorbar_label='Raster Elevation')
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+    :hide:
+
+    >>> bounds = ax.axis('equal')
+    >>> plt.tight_layout()
+    >>> plt.savefig('sample_elevations.png', dpi=300)
+    >>> plt.close()
+
+.. _fig-sample-elevations:
+.. figure:: figures/sample_elevations.png
+   :width: 640
+   :alt: Net1 with elevations sampled from raster.
+
+   Net1 with elevations sampled from raster.

documentation/installation.rst

@@ -272,6 +272,8 @@ The following Python packages are optional:
   https://pypi.org/project/utm/
 * geopandas :cite:p:`jvfm21`: used to work with geospatial data,
   https://geopandas.org/
+* rasterio :cite:p:`rasterio`: used to work with raster data,
+  https://rasterio.readthedocs.io/
 * rtree :cite:p:`rtree`: used for overlay operations in geopandas,
   https://rtree.readthedocs.io/
 * openpyxl :cite:p:`gacl18`: used to read/write to Microsoft® Excel® spreadsheets,

documentation/model_io.rst

@@ -206,27 +206,29 @@ GeoJSON files
 GeoJSON files are commonly used to store geographic data structures. 
 More information on GeoJSON files can be found at https://geojson.org.
 
-To use GeoJSON files in WNTR, a set of valid base column names are required.
-Valid base GeoJSON column names can be obtained using the 
-:class:`~wntr.network.io.valid_gis_names` function.
-The following example returns valid base GeoJSON column names for junctions.
+When reading GeoJSON files into WNTR, only a set of valid column names can be used.
+Valid GeoJSON column names can be obtained using the 
+:class:`~wntr.network.io.valid_gis_names` function. By default, the function
+returns all column names, both required and optional.
+The following example returns valid GeoJSON column names for junctions.
 
 .. doctest::
     :skipif: gpd is None
 
     >>> geojson_column_names = wntr.network.io.valid_gis_names()
     >>> print(geojson_column_names['junctions'])
-    ['name', 'elevation', 'coordinates', 'base_demand', 'demand_pattern', 'demand_category', 'emitter_coefficient', 'initial_quality', 'minimum_pressure', 'required_pressure', 'pressure_exponent', 'tag']
+    ['name', 'elevation', 'geometry', 'emitter_coefficient', 'initial_quality', 'minimum_pressure', 'required_pressure', 'pressure_exponent', 'tag']
 
-A minimal list of valid column names can also be obtained by setting ``complete_list`` to False.
-Column names that are optional (i.e., ``initial_quality``) and not included in the GeoJSON file are defined using default values.
+A minimal list of required column names can also be obtained by setting ``complete_list`` to False.
+Column names that are optional (i.e., ``initial_quality``) and not included in the GeoJSON file are 
+defined using default values.
 
 .. doctest::
     :skipif: gpd is None
 
     >>> geojson_column_names = wntr.network.io.valid_gis_names(complete_list=False)
     >>> print(geojson_column_names['junctions'])
-    ['name', 'elevation', 'coordinates']
+    ['name', 'elevation', 'geometry']
 
 Note that GeoJSON files can contain additional custom column names that are assigned to WaterNetworkModel objects.
 
@@ -253,7 +255,7 @@ Note that patterns, curves, sources, controls, and options are not stored in the
 
 The :class:`~wntr.network.io.read_geojson` function creates a WaterNetworkModel from a 
 dictionary of GeoJSON files.  
-Valid base column names and additional custom attributes are added to the model.
+Valid column names and additional custom attributes are added to the model.
 The function can also be used to append information from GeoJSON files into an existing WaterNetworkModel.
 
 .. doctest::
@@ -300,20 +302,21 @@ To use Esri Shapefiles in WNTR, several formatting requirements are enforced:
   assumed that the first 10 characters of each attribute are unique.  
 
 * To create WaterNetworkModel from Shapefiles, a set of valid field names are required.
-  Valid base Shapefiles field names can be obtained using the 
-  :class:`~wntr.network.io.valid_gis_names` function.
-  For Shapefiles, the `truncate` input parameter should be set to 10 (characters).
-  The following example returns valid base Shapefile field names for junctions.
-  Note that attributes like ``base_demand`` are truncated to ``base_deman``. 
+  Valid Shapefiles field names can be obtained using the 
+  :class:`~wntr.network.io.valid_gis_names` function. By default, the function
+  returns all column names, both required and optional.
+  For Shapefiles, the `truncate_names` input parameter should be set to 10 (characters).
+  The following example returns valid Shapefile field names for junctions.
+  Note that attributes like ``minimum_pressure`` are truncated to ``minimum_pr``. 
 
   .. doctest::
       :skipif: gpd is None
 
       >>> shapefile_field_names = wntr.network.io.valid_gis_names(truncate_names=10)
       >>> print(shapefile_field_names['junctions'])
-      ['name', 'elevation', 'coordinate', 'base_deman', 'demand_pat', 'demand_cat', 'emitter_co', 'initial_qu', 'minimum_pr', 'required_p', 'pressure_e', 'tag']
+      ['name', 'elevation', 'geometry', 'emitter_co', 'initial_qu', 'minimum_pr', 'required_p', 'pressure_e', 'tag']
 
-  A minimal list of valid field names can also be obtained by setting ``complete_list`` to False.
+  A minimal list of required field names can also be obtained by setting ``complete_list`` to False.
   Field names that are optional (i.e., ``initial_quality``) and not included in the Shapefile are defined using default values.
 
   .. doctest::
@@ -322,7 +325,7 @@ To use Esri Shapefiles in WNTR, several formatting requirements are enforced:
       >>> shapefile_field_names = wntr.network.io.valid_gis_names(complete_list=False, 
       ...    truncate_names=10)
       >>> print(shapefile_field_names['junctions'])
-      ['name', 'elevation', 'coordinate']
+      ['name', 'elevation', 'geometry']
 
 * Shapefiles can contain additional custom field names that are assigned to WaterNetworkModel objects.
 
@@ -349,7 +352,7 @@ Note that patterns, curves, sources, controls, and options are not stored in the
 
 The :class:`~wntr.network.io.read_shapefile` function creates a WaterNetworkModel from a dictionary of
 Shapefile directories.
-Valid base field names and additional custom field names are added to the model.
+Valid field names and additional custom field names are added to the model.
 The function can also be used to append information from Shapefiles into an existing WaterNetworkModel.
 
 .. doctest::

documentation/references.bib

@@ -128,6 +128,14 @@ @misc{jvfm21
     year = "2021"
 }
 
+@misc{rasterio,
+  author =    {Sean Gillies and others},
+  organization = {Mapbox},
+  title =     {{Rasterio: geospatial raster I/O for {Python} programmers}},
+  year =      {2013--},
+  url = {"https://github.com/rasterio/rasterio"}
+}
+
 @article{jcmg11,
     author = "Joyner, David and \v{C}ert\'{i}k, Ond\v{r}ej and Meurer, Aaron and Granger, Brian E.",
     address = "New York, NY, USA",

requirements.txt

@@ -13,6 +13,7 @@ utm
 openpyxl
 geopandas<1.0
 fiona<1.10
+rasterio
 rtree
 
 # Documentation

wntr/epanet/io.py

@@ -674,6 +674,8 @@ def _write_tanks(self, f, wn, version=2.2):
         f.write('\n'.encode(sys_default_enc))
 
     def _read_pipes(self):
+        darcy_weisbach = self.wn.options.hydraulic.headloss == "D-W"
+
         for lnum, line in self.sections['[PIPES]']:
             line = line.split(';')[0]
             current = line.split()
@@ -701,7 +703,8 @@ def _read_pipes(self):
                         current[2],
                         to_si(self.flow_units, float(current[3]), HydParam.Length),
                         to_si(self.flow_units, float(current[4]), HydParam.PipeDiameter),
-                        float(current[5]),
+                        to_si(self.flow_units, float(current[5]), HydParam.RoughnessCoeff, 
+                              darcy_weisbach=darcy_weisbach),
                         minor_loss,
                         link_status,
                         check_valve)
@@ -711,6 +714,8 @@ def _read_pipes(self):
                 raise ENValueError(211, str(e.args[0]), line_num=lnum) from e
 
     def _write_pipes(self, f, wn):
+        darcy_weisbach = wn.options.hydraulic.headloss == "D-W"
+
         f.write('[PIPES]\n'.encode(sys_default_enc))
         f.write(_PIPE_LABEL.format(';ID', 'Node1', 'Node2', 'Length', 'Diameter',
                                    'Roughness', 'Minor Loss', 'Status').encode(sys_default_enc))
@@ -723,7 +728,9 @@ def _write_pipes(self, f, wn):
                  'node2': pipe.end_node_name,
                  'len': from_si(self.flow_units, pipe.length, HydParam.Length),
                  'diam': from_si(self.flow_units, pipe.diameter, HydParam.PipeDiameter),
-                 'rough': pipe.roughness,
+                 'rough': from_si(self.flow_units, pipe.roughness, 
+                                  HydParam.RoughnessCoeff, 
+                                  darcy_weisbach=darcy_weisbach),
                  'mloss': pipe.minor_loss,
                  'status': str(pipe.initial_status),
                  'com': ';'}