diff --git a/documentation/figures/assign_elevations.png b/documentation/figures/assign_elevations.png new file mode 100644 index 000000000..847da7f1a Binary files /dev/null and b/documentation/figures/assign_elevations.png differ diff --git a/documentation/gis.rst b/documentation/gis.rst index e57d6886d..ae79f614f 100644 --- a/documentation/gis.rst +++ b/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 _`. + +.. 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. diff --git a/examples/data/elevation.tif b/examples/data/elevation.tif new file mode 100644 index 000000000..6a40d473c Binary files /dev/null and b/examples/data/elevation.tif differ diff --git a/wntr/gis/geospatial.py b/wntr/gis/geospatial.py index bcb7c2f0f..8b2989bc2 100644 --- a/wntr/gis/geospatial.py +++ b/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( diff --git a/wntr/tests/test_gis.py b/wntr/tests/test_gis.py index d7bf91135..c9857f053 100644 --- a/wntr/tests/test_gis.py +++ b/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()