Geographical analysis for hydrological modelling

This notebook shows how to delineate a catchment and extract properties from a digital elevation model (DEM) and a land-use data set. The processes are hosted on the Raven server, which in the background connects to a GeoServer instance to query watershed contours, DEM and land-use data.

We connect to Raven’s Web Processing Service interface using birdy’s WPSClient.

from birdy import WPSClient
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import json
import os

# Set environment variable WPS_URL to "http://localhost:9099" to run on the default local server
url = os.environ.get("WPS_URL", "https://pavics.ouranos.ca/twitcher/ows/proxy/raven/wps")

# Connect to the PAVICS-Hydro Raven WPS server
wps = WPSClient(url)

We first extract the watershed contour for the point of interest. The process looks into the HydroSheds databast to finds the watershed enclosing the given location. The location parameter identifies the outlet of the watershed, and aggregate_upstream determines whether or not we want the service to return the union of all upstream basins. Here we set it to False to reduce the size of the basin and speed-up computations.

The output of the hydrosheds-select process is a GeoJSON geometry.

r_select = wps.hydrobasins_select(location="-71.291660, 50.492758", aggregate_upstream=False)

# Get GeoJSON polygon of the delineated catchment.
# We can either get links to the files stored on the server, or get the data directly.
[feature_url, upstream_basins_url] = r_select.get(asobj=False)
[feature, upstream_basins] = r_select.get(asobj=True)

We can now plot the outline of the watershed by loading it into GeoPandas.

# df = gpd.read_file(feature_url)
df = gpd.GeoDataFrame.from_features([feature])
df.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x7f6bbe091940>

Now that we have delineated a catchment, lets find the zonal statistics and other properties of the catchment. We can pass the catchment outline either as a link to a file (e.g. feature_url), or pass the data directly as a string.

# Here we are using the geojson file created on the server in the last step as the input value to a process
# computing watershed properties.
#crs=4326
#projected_crs=32198
resp = wps.shape_properties(shape=json.dumps(feature))

Now, let’s extract the data from the WPS service response:

[properties, ]=resp.get(asobj=True)
prop = properties[0]
print(prop)

area = prop['area']/1000000.0
longitude = prop['centroid'][0]
latitude = prop['centroid'][1]
gravelius = prop['gravelius']
perimeter = prop['perimeter']

shapeProperties = {'area':area, 'longitude':longitude, 'latitude':latitude, 'gravelius':gravelius, 'perimeter':perimeter}
shapeProperties

{'id': '96929', 'gml_id': 'USGS_HydroBASINS_lake_na_lev12.96929', 'HYBAS_ID': 7120270182, 'NEXT_DOWN': 7120270181, 'NEXT_SINK': 7120034330, 'MAIN_BAS': 7120034330, 'DIST_SINK': 490.9, 'DIST_MAIN': 490.9, 'SUB_AREA': 29.0, 'UP_AREA': 9419.6, 'PFAF_ID': 724089370000, 'SIDE': 'R', 'LAKE': 0, 'ENDO': 0, 'COAST': 0, 'ORDER': 1, 'SORT': 96929, 'area': 28764849.46504176, 'centroid': [-71.3409808346398, 50.478880424555875], 'perimeter': 33017.42666655582, 'gravelius': 1.7366297521025682}

{'area': 28.76484946504176,
 'longitude': -71.3409808346398,
 'latitude': 50.478880424555875,
 'gravelius': 1.7366297521025682,
 'perimeter': 33017.42666655582}

Note that these properties are a mix of the properties of the original file where the shape is stored, and properties computed by the process (area, centroid, perimeter and gravelius). Note also that the computed area is in m², while the “SUB_AREA” property is in km², and that there are slight differences between the two values.

Now we’ll extract the land-use properties of the watershed. We pass the link to the watershed outline to a process using in the background the North American Land Change Monitoring System dataset, and retrieve properties over the region.

# Use the geoserver to extract the land cover over the appropriate bounding box (automatic)
resp = wps.nalcms_zonal_stats(shape=feature_url, select_all_touching=True, band=1, simple_categories=True)

# Note that geojson needs to be installed for this to work.
# $ pip install -r requirements_extra.txt
features, statistics  = resp.get(asobj=True)
print(features[0]['properties'], '\n\n')
print(statistics)

{'id': '96929', 'gml_id': 'USGS_HydroBASINS_lake_na_lev12.96929', 'HYBAS_ID': 7120270182, 'NEXT_DOWN': 7120270181, 'NEXT_SINK': 7120034330, 'MAIN_BAS': 7120034330, 'DIST_SINK': 490.9, 'DIST_MAIN': 490.9, 'SUB_AREA': 29.0, 'UP_AREA': 9419.6, 'PFAF_ID': 724089370000, 'SIDE': 'R', 'LAKE': 0, 'ENDO': 0, 'COAST': 0, 'ORDER': 1, 'SORT': 96929, '1': 13532, '2': 335, '5': 4367, '6': 847, '8': 10239, '10': 60, '12': 380, '14': 21, '16': 39, '18': 3248, 'count': 33068, 'nodata': 8.0, 'nan': 0, 'Ocean': 0, 'Forest': 19081, 'Shrubs': 10239, 'Grass': 479, 'Wetland': 21, 'Crops': 0, 'Urban': 0, 'Water': 3248, 'SnowIce': 0}


[{'Ocean': 0, 'Forest': 19081, 'Shrubs': 10239, 'Grass': 479, 'Wetland': 21, 'Crops': 0, 'Urban': 0, 'Water': 3248, 'SnowIce': 0}]

features

{"features": [{"geometry": {"coordinates": [[[[1981874.096748, 967399.091191], [1981560.497165, 967640.799616], [1981434.936988, 967882.780767], [1980961.263774, 969007.129417], [1980150.453239, 969613.136897], [1979557.788254, 970982.006618], [1979190.993435, 971263.662311], [1978968.917959, 971674.108713], [1979068.737029, 972221.498664], [1978911.276644, 972753.897266], [1978998.518432, 973199.540058], [1978148.536013, 975486.000658], [1978229.072674, 975929.06883], [1977903.742496, 976905.067086], [1977984.173254, 977348.115496], [1977736.398881, 978051.400822], [1977900.425286, 978054.503503], [1978726.767107, 977503.147592], [1979063.162155, 977488.410214], [1979585.086481, 977194.438111], [1979728.314495, 977383.335657], [1979905.626466, 977961.050647], [1980249.760066, 977888.82216], [1980777.020525, 976821.614008], [1981210.625424, 976808.412368], [1981599.965266, 976584.172313], [1981738.25575, 977013.547429], [1981974.993178, 977251.010851], [1982328.051065, 978903.0378], [1982681.072971, 978846.538719], [1982996.790791, 978315.034465], [1983354.336434, 977981.589513], [1983659.560071, 977324.816337], [1983830.052137, 976894.322922], [1983556.471154, 976787.837664], [1982897.423789, 976458.687678], [1982761.508473, 976042.300741], [1982257.567097, 975700.92043], [1981989.759672, 974591.183905], [1982132.826199, 973944.08598], [1981931.599976, 972702.854607], [1981635.845096, 972309.345787], [1981135.976984, 971194.624495], [1981172.131364, 971160.205611], [1981342.839416, 971032.658661], [1981721.483719, 970283.259585], [1982217.756236, 968683.141738], [1982076.182758, 967974.105719], [1981874.096748, 967399.091191]]]], "type": "MultiPolygon"}, "id": "0", "properties": {"1": 13532, "10": 60, "12": 380, "14": 21, "16": 39, "18": 3248, "2": 335, "5": 4367, "6": 847, "8": 10239, "COAST": 0, "Crops": 0, "DIST_MAIN": 490.9, "DIST_SINK": 490.9, "ENDO": 0, "Forest": 19081, "Grass": 479, "HYBAS_ID": 7120270182, "LAKE": 0, "MAIN_BAS": 7120034330, "NEXT_DOWN": 7120270181, "NEXT_SINK": 7120034330, "ORDER": 1, "Ocean": 0, "PFAF_ID": 724089370000, "SIDE": "R", "SORT": 96929, "SUB_AREA": 29.0, "Shrubs": 10239, "SnowIce": 0, "UP_AREA": 9419.6, "Urban": 0, "Water": 3248, "Wetland": 21, "count": 33068, "gml_id": "USGS_HydroBASINS_lake_na_lev12.96929", "id": "96929", "nan": 0, "nodata": 8.0}, "type": "Feature"}], "type": "FeatureCollection"}

We now have the statistics from the NALCMS zonal_stats toolbox regarding the land use, from which we calculate the ratio of each land-use component.

# total = sum(lu.values())
lu = statistics[0]
total = sum(lu.values())

landUse = {k: (v / total) for (k,v) in lu.items()}
landUse

{'Ocean': 0.0,
 'Forest': 0.5770231039071005,
 'Shrubs': 0.3096346921495101,
 'Grass': 0.014485303011975323,
 'Wetland': 0.0006350550381033022,
 'Crops': 0.0,
 'Urban': 0.0,
 'Water': 0.09822184589331075,
 'SnowIce': 0.0}

The next step will be to collect terrain data, such as elevation, slope and aspect. We will do this using the terrain_analysis WPS service:

resp = wps.terrain_analysis(shape=feature_url, select_all_touching=True, projected_crs=3978)

Now let’s extract the properties from the WPS response. Use asobj=True to have Birdy preprocess the data and return the data directly.

properties, dem = resp.get(asobj=True)

elevation=properties[0]['elevation']
slope=properties[0]['slope']
aspect=properties[0]['aspect']

terrain_data={'elevation':elevation, 'slope':slope,'aspect':aspect}

Finally, display all the extracted parameters for the user:

all_properties={**shapeProperties, **landUse, **terrain_data}
print(all_properties)

{'area': 28.76484946504176, 'longitude': -71.3409808346398, 'latitude': 50.478880424555875, 'gravelius': 1.7366297521025682, 'perimeter': 33017.42666655582, 'Ocean': 0.0, 'Forest': 0.5770231039071005, 'Shrubs': 0.3096346921495101, 'Grass': 0.014485303011975323, 'Wetland': 0.0006350550381033022, 'Crops': 0.0, 'Urban': 0.0, 'Water': 0.09822184589331075, 'SnowIce': 0.0, 'elevation': 490.04395604395603, 'slope': 3.9660612485567572, 'aspect': 116.79663053081183}

Note here that while the feature outline is defined above in terms of geographic coordinates (latitude, longitude), the DEM is projected onto a 2D cartesion coordinate system (here NAD83, the Canada Atlas Lambert projection). This is necessary to perform slope calculations.

For more information on this, see: https://en.wikipedia.org/wiki/Map_projection

# NBVAL_SKIP
import cartopy.crs as ccrs
import rasterio
from rasterio.plot import show
import xarray as xr
import rasterio
from rasterio.io import MemoryFile

with MemoryFile(dem) as mem:
    with mem.open(driver='gtiff') as src:
        crs = ccrs.LambertConformal(central_latitude=49, central_longitude=-95, standard_parallels=(49, 77))
        da = xr.open_rasterio(src)
        da.name = 'Elevation'
        da.attrs['units'] = 'm'
        ax = plt.subplot(projection=crs)
        da.where(da!=-32768).sel(band=1).plot.imshow(ax=ax, transform=crs)
        plt.show()