update test to use net1. add user manual documentation. figure and da…

…ta added for raster example

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/elevation.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.
@@ -822,3 +829,88 @@ 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.
+
+Find the intersect between geometries
+--------------------------------------
+
+The :class:`~wntr.gis.sample_raster` function can be used to query a raster file at point geometries,
+such as the nodes of a water network.
+
+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')
+
+Assign elevations to nodes
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Elevation is an essential attribute for accurate simulation of pressure in a water network. Elevation data is
+commonly provided in GeoTIFF files, such as that provided by the 
+`USGS 3D Elevation Program <https://www.usgs.gov/3d-elevation-program>_`. 
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> junctions = wn_gis.junctions
+    >>> print(junctions["elevation"])
+    name
+    10    216.408
+    11    216.408
+    12    213.360
+    13    211.836
+    21    213.360
+    22    211.836
+    23    210.312
+    31    213.360
+    32    216.408
+    Name: elevation, dtype: float64
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> elevation_data_path = 'data/elevation.tif' # doctest: +SKIP
+    >>> junctions["elevation"] = wntr.gis.sample_raster(junctions, elevation_data_path)
+    >>> print(junctions["elevation"])
+    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
+    Name: elevation, dtype: float64
+
+The assigned elevations can be plotted as follows. The 
+resulting :numref:`fig-assingn-elevations` illustrates Net1 with the elevations queried from the raster file.
+Note that to use these elevations in a simulation, they would need to be added to the water network object directly.
+Tanks, in addition to junctions, would need their elevations updated.
+
+.. doctest::
+    :skipif: gpd is None or rasterio is None
+
+    >>> ax = wntr.graphics.plot_network(wn, node_attribute=junctions["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('assign_elevations.png', dpi=300)
+    >>> plt.close()
+
+.. _fig-assign-elevations:
+.. figure:: figures/assign_elevations.png
+   :width: 640
+   :alt: Net1 junctions with elevations from raster.
+
+    Net1 junctions with elevations assigned from raster.

wntr/gis/geospatial.py

@@ -21,10 +21,10 @@
     has_geopandas = False
 
 try:
-    import rasterio as rio
+    import rasterio
     has_rasterio = True
 except ModuleNotFoundError:
-    rio = None
+    rasterio = None
     has_rasterio = False
 
 
@@ -324,7 +324,7 @@ def sample_raster(A, filepath, bands=1):
     assert isinstance(filepath, str)
     assert isinstance(bands, (int, list))
 
-    with rio.open(filepath) as raster:
+    with rasterio.open(filepath) as raster:
         xys = zip(A.geometry.x, A.geometry.y)
 
         values = np.array(

wntr/tests/test_gis.py

@@ -24,10 +24,10 @@
     has_geopandas = False
 
 try:
-    import rasterio as rio
+    import rasterio
     has_rasterio = True
 except ModuleNotFoundError:
-    rio = None
+    rasterio = None
     has_rasterio = False
 
 testdir = dirname(abspath(str(__file__)))
@@ -323,36 +323,29 @@ def test_snap_points_to_lines(self):
 class TestRaster(unittest.TestCase):
     @classmethod
     def setUpClass(self):
-
-        # raster testing
-        points = [
-            (-120.5, 38.5),
-            (-120.6, 38.6),
-            (-120.55, 38.65),
-            (-120.65, 38.55),
-            (-120.7, 38.7)
-        ]
-        point_geometries = [Point(xy) for xy in points]
-        points = gpd.GeoDataFrame(geometry=point_geometries, crs="EPSG:4326")
-        points.index = ["A", "B", "C", "D", "E"]
+        # use net1 junctions as example points
+        inp_file = join(ex_datadir, "Net1.inp")
+        wn = wntr.network.WaterNetworkModel(inp_file)
+        wn_gis = wn.to_gis(crs="EPSG:4326")
+        points = wn_gis.junctions
         self.points = points
 
-        # create example raster
-        minx, miny, maxx, maxy = points.total_bounds 
-        raster_width = 100
-        raster_height = 100 
+        # create test raster file
+        min_lon, min_lat = -180, -90  
+        max_lon, max_lat = 180, 90
 
-        x = np.linspace(0, 1, raster_width)
-        y = np.linspace(0, 1, raster_height)
-        raster_data = np.cos(y)[:, np.newaxis] * np.sin(x) # arbitrary values
+        resolution = 1.0
+
+        lon_values = np.arange(min_lon, max_lon, resolution)
+        lat_values = np.arange(max_lat, min_lat, -resolution)  # Decreasing order for latitudes
+        raster_data = np.outer(lon_values,lat_values) # value is product of coordinate
 
-        transform = rio.transform.from_bounds(minx, miny, maxx, maxy, raster_width, raster_height)
-        self.transform = transform
+        transform = rasterio.transform.from_bounds(min_lon, min_lat, max_lon, max_lat, raster_data.shape[0], raster_data.shape[1])
 
-        with rio.open(
-            "test_raster.tif", "w", driver="GTiff", height=raster_height, width=raster_width, 
-            count=1, dtype=raster_data.dtype, crs="EPSG:4326", transform=transform) as dst:
-            dst.write(raster_data, 1) 
+        with rasterio.open(
+            "test_raster.tif", "w", driver="GTiff", height=raster_data.shape[1], width=raster_data.shape[0], 
+            count=1, dtype=raster_data.dtype, crs="EPSG:4326", transform=transform) as file:
+            file.write(raster_data, 1) 
 
     @classmethod
     def tearDownClass(self):
@@ -363,7 +356,8 @@ def test_sample_raster(self):
         assert (raster_values.index == self.points.index).all()
 
         # self.points.plot(column=values, legend=True)
-        expected_values = np.array([0.000000, 0.423443, 0.665369, 0.174402, 0.000000])
+        # values should be product of coordinates
+        expected_values = self.points.apply(lambda row: row.geometry.x * row.geometry.y, axis=1)
         assert np.isclose(raster_values.values, expected_values, atol=1e-5).all()