Accessing Deltares global water availability on Azure¶
Deltares has produced a hydrological model approach to simulate historical daily reservoir variations for 3,236 locations across the globe for the period 1970-2020 using the distributed wflow_sbm model. The model outputs long-term daily information on reservoir volume, inflow and outflow dynamics, as well as information on upstream hydrological forcing.
They hydrological model was forced with 5 different precipitation products. Two products (ERA5 and CHIRPS) are available at the global scale, while for Europe, USA and Australia a regional product was use (i.e. EOBS, NLDAS and BOM, respectively). Using these different precipitation products, it becomes possible to assess the impact of uncertainty in the model forcing. A different number of basins upstream of reservoirs are simulated, given the spatial coverage of each precipitation product.
Environment setup¶
In [1]:
import warnings
warnings.simplefilter("ignore")
import fsspec
import xarray as xr
import pystac_client
import planetary_computer
import requests
import geopandas
import contextily
import dask.distributed
Create a local Dask cluster¶
Enable parallel reads and processing of data using Dask and xarray.
In [2]:
client = dask.distributed.Client(processes=False)
print(client.dashboard_link)
/user/taugspurger@microsoft.com/proxy/8787/status
Data access¶
The collection is made up of five netCDF files, each corresponding to a different precipitation product for hydrological model forcing. Each variation covers a different spatial extent, with different number of reservoirs modeled. We'll choose one of the smaller files using a BOM dataset covering Australia here. The delatres:reservoir key can be used to filter to a specific dataset.
The other files available from the same root store are CHIRPS, EOBS, ERA5, and NLDAS.
We'll open the dataset and look at the variables and dimensions it contains.
In [3]:
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)
In [4]:
item = catalog.get_collection("deltares-water-availability").get_item("BOM")
item
Out[4]:
Item: BOM
| id: BOM |
| bbox: [115.95417023, -43.04583359, 153.3374939, -17.17083359] |
| datetime: 1982-01-01T00:00:00Z |
| end_datetime: 2020-12-31T00:00:00Z |
| cube:variables: {'P': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Average precipitation upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average precipitation upstream of reservoir'}, 'ETa': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Average simulated actual evapotransporation upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average simulated actual evapotransporation upstream of reservoir'}, 'PET': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Average potential evapotranspiration upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average potential evapotranspiration upstream of reservoir'}, 'Melt': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Average simulated snow melt upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average simulated snow melt upstream of reservoir'}, 'Snow': {'type': 'data', 'unit': 'mm', 'attrs': {'units': 'mm', 'description': 'Average simulated snow depth upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average simulated snow depth upstream of reservoir'}, 'Temp': {'type': 'data', 'unit': 'degrees C', 'attrs': {'units': 'degrees C', 'description': 'Average surface temperature upstream of reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Average surface temperature upstream of reservoir'}, 'P_res': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Precipitation reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Precipitation reservoir'}, 'S_res': {'type': 'data', 'unit': 'm3', 'attrs': {'units': 'm3', 'description': 'Simulated reservoir volume'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Simulated reservoir volume'}, 'Ea_res': {'type': 'data', 'unit': 'mm per day', 'attrs': {'units': 'mm per day', 'description': 'Simulated actual evaporation reservoir'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Simulated actual evaporation reservoir'}, 'Qin_res': {'type': 'data', 'unit': 'm3 per s', 'attrs': {'units': 'm3 per s', 'description': 'Simulated reservoir inflow (surface+subsurface)'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Simulated reservoir inflow (surface+subsurface)'}, 'FracFull': {'type': 'data', 'unit': 'm3', 'attrs': {'units': 'm3', 'description': 'Simulated reservoir volume'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Simulated reservoir volume'}, 'Qout_res': {'type': 'data', 'unit': 'm3 per s', 'attrs': {'units': 'm3 per s', 'description': 'Simulated reservoir outflow'}, 'shape': [14245, 116, 5], 'dimensions': ['time', 'GrandID', 'ksathorfrac'], 'description': 'Simulated reservoir outflow'}, 'latitude': {'type': 'data', 'unit': 'degrees', 'attrs': {'units': 'degrees', 'description': 'Latitude of reservoir'}, 'shape': [116], 'dimensions': ['GrandID'], 'description': 'Latitude of reservoir'}, 'longitude': {'type': 'data', 'unit': 'degrees', 'attrs': {'units': 'degrees', 'description': 'Longitude of reservoir'}, 'shape': [116], 'dimensions': ['GrandID'], 'description': 'Longitude of reservoir'}} |
| start_datetime: 1982-01-01T00:00:00Z |
| cube:dimensions: {'time': {'step': 'P1DT0H0M0S', 'type': 'temporal', 'extent': ['1982-01-01T00:00:00Z', '2020-12-31T00:00:00Z']}, 'GrandID': {'type': 'identifier', 'extent': [5822, 6752], 'description': 'GrandID number of the reservoir of interest'}, 'ksathorfrac': {'type': 'level', 'values': [5, 20, 50, 100, 250], 'description': 'Five different value lateral anisotropy values used'}} |
| deltares:reservoir: BOM |
STAC Extensions
| https://stac-extensions.github.io/datacube/v2.0.0/schema.json |
Assets
Asset: Flood Map
| href: https://deltaresreservoirssa.blob.core.windows.net/reservoirs/v2021.12/reservoirs_BOM.nc?st=2023-10-10T11%3A43%3A01Z&se=2023-10-18T11%3A43%3A02Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-10-11T11%3A43%3A00Z&ske=2023-10-18T11%3A43%3A00Z&sks=b&skv=2021-06-08&sig=PLo049qGIsJeYKHbZfTY4Pi4%2BysZgKz7qIh0l0Dy9CE%3D |
| type: application/x-netcdf |
| title: Flood Map |
| description: Inundation maps of flood depth using a model that takes into account water level attenuation and is forced by sea level. |
| roles: ['data'] |
| owner: BOM |
Asset: Index file
| href: https://deltaresreservoirssa.blob.core.windows.net/references/reservoirs/BOM.json?st=2023-10-10T11%3A43%3A02Z&se=2023-10-18T11%3A43%3A02Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-10-11T11%3A43%3A01Z&ske=2023-10-18T11%3A43%3A01Z&sks=b&skv=2021-06-08&sig=Uo3bwXYV/PY1xAUKYZQBLmnGUzGjzTcO/olsV/aiSOU%3D |
| type: application/json |
| title: Index file |
| description: Kerchunk index file. |
| roles: ['index'] |
| owner: BOM |
Links
Link:
| rel: self |
| href: https://planetarycomputer.microsoft.com/api/stac/v1/collections/deltares-water-availability/items/BOM |
| type: application/json |
Link:
| rel: collection |
| href: https://planetarycomputer.microsoft.com/api/stac/v1/collections/deltares-water-availability |
| type: application/json |
Link:
| rel: parent |
| href: https://planetarycomputer.microsoft.com/api/stac/v1/collections/deltares-water-availability |
| type: application/json |
Link: Deltares Global Water Availability
| rel: root |
| href: https://planetarycomputer.microsoft.com/api/stac/v1/collections/deltares-water-availability |
| type: application/json |
| title: Deltares Global Water Availability |
Each item has two assets: data pointing to the NetCDF file, and index pointing to an index file enabling cloud-optimzed access. When accessing the data from Azure, we recommend using the index file.
In [5]:
index = planetary_computer.sign(requests.get(item.assets["index"].href).json())
store = fsspec.get_mapper("reference://", fo=index)
ds = xr.open_dataset(store, engine="zarr", consolidated=False, chunks={})
ds
Out[5]:
<xarray.Dataset>
Dimensions: (time: 14245, GrandID: 116, ksathorfrac: 5)
Coordinates:
* GrandID (GrandID) float64 6.675e+03 6.666e+03 ... 6.594e+03 6.706e+03
* ksathorfrac (ksathorfrac) float64 5.0 20.0 50.0 100.0 250.0
* time (time) datetime64[ns] 1982-01-01 1982-01-02 ... 2020-12-31
Data variables: (12/14)
ETa (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
Ea_res (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
FracFull (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
Melt (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
P (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
PET (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
... ...
Qout_res (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
S_res (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
Snow (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
Temp (time, GrandID, ksathorfrac) float32 dask.array<chunksize=(1, 116, 5), meta=np.ndarray>
latitude (GrandID) float32 dask.array<chunksize=(116,), meta=np.ndarray>
longitude (GrandID) float32 dask.array<chunksize=(116,), meta=np.ndarray>Inspect the reservoirs¶
First, we'll plot a simple map of where all the modeled reservoirs are located in Australia.
In [6]:
# Just interested in the GrandID dimension
df = ds.drop_dims(["time", "ksathorfrac"]).to_pandas()
In [10]:
geometry = geopandas.GeoSeries(
geopandas.points_from_xy(df.longitude, df.latitude), crs="WGS84"
)
ax = geometry.plot(figsize=(12, 12))
contextily.add_basemap(
ax, crs=geometry.crs.to_string(), source=contextily.providers.Esri.NatGeoWorldMap
)
ax.set_axis_off()
ax.set(ylim=(-44, -9), xlim=(112, 156));
Mean precipitation upstream of all reservoirs, single year¶
Now we'll pluck a single year in the series and plot the mean precipitation upstream of all the reservoir basins.
In [11]:
precip_yr = ds.P.sel(time="1989").load()
In [12]:
g = precip_yr.groupby("time")
p = g.mean(...).plot(figsize=(15, 10))
Inspect a single reservoir¶
We'll select the first reservoir by Id and a KsatHorFrac of 50.
In [13]:
single_res = ds.sel(GrandID=ds.GrandID[0], ksathorfrac=50)
single_res
Out[13]:
<xarray.Dataset>
Dimensions: (time: 14245)
Coordinates:
GrandID float64 6.675e+03
ksathorfrac float64 50.0
* time (time) datetime64[ns] 1982-01-01 1982-01-02 ... 2020-12-31
Data variables: (12/14)
ETa (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
Ea_res (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
FracFull (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
Melt (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
P (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
PET (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
... ...
Qout_res (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
S_res (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
Snow (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
Temp (time) float32 dask.array<chunksize=(1,), meta=np.ndarray>
latitude float32 dask.array<chunksize=(), meta=np.ndarray>
longitude float32 dask.array<chunksize=(), meta=np.ndarray>Observe the impact of KsatHorFrac¶
Here we'll select a new reservoir and plot how the different KsatHorFrac values affect volume in the model. Because this touches every file in the timeseries, you might want a to run it with a cluster.
single_res_ksats = ds.sel(GrandID=6662).groupby("time.year").mean()
p = single_res_ksats.FracFull.plot.line(figsize=(12, 8), hue="ksathorfrac", x="year")
Inspect all reservoirs¶
We can also plot out all the reservoirs' FracFull variable and see how the different KsatHorFrac values affect volume in the model.
y = ds.FracFull.groupby("time.year").mean().compute()
p = y.plot.line(col="GrandID", x="year", col_wrap=6)