diff --git a/examples/update_model_water_demand.ipynb b/examples/update_model_water_demand.ipynb index d8552257..5b379dee 100644 --- a/examples/update_model_water_demand.ipynb +++ b/examples/update_model_water_demand.ipynb @@ -44,6 +44,7 @@ "source": [ "Here, you see we will be adding the required information to simulate water demand and allocation in wflow, by running three different setup functions:\n", "- **setup_allocation_areas**: Adds a map that defines the different regions that will be used to allocate water. By default this is a mix between administrative boundaries and catchment boundaries\n", + "- **setup_allocation_surfacewaterfrac**: prepare the fraction of surface water used for allocation (can be reduced if groundwater or non conventional water sources are also present in the basin).\n", "- **setup_lulcmaps_with_paddy** (or **setup_lulcmaps**): update the landuse to include new parameters (crop coefficient kc and soil water pressure heads h) and add the paddies (rice fields). To allow for water to pool on the surface of the rice fields, soil parameters will also be updated to include an additional thin layer with limited vertical conductivity.\n", "- **setup_domestic_demand**: Add domestic water demands (gross and net) from gridded data and downscaled using high resolution population map.\n", "- **setup_other_demand**: Adds maps to the wflow schematization that describe how much water is demanded (gross and net amounts) by different sources: domestic (dom), industry (ind), and livestock (lsk). In our case, as we downscale domestic with population, we will here add industry and livestock.\n", @@ -73,7 +74,7 @@ "outputs": [], "source": [ "# NOTE: copy this line (without !) to your shell for more direct feedback\n", - "! hydromt build wflow \"./wflow_piave_water_demand\" -r \"{'basin': [12.2051, 45.8331], 'bounds': [11.70, 45.35, 12.95, 46.70]}\" -i wflow_build.yml -d \"artifact_data==v0.0.6\" -vv" + "! hydromt build wflow \"./wflow_piave_water_demand\" -r \"{'basin': [12.2051, 45.8331], 'bounds': [11.70, 45.35, 12.95, 46.70]}\" -i wflow_build.yml -d artifact_data -vv" ] }, { @@ -89,7 +90,7 @@ "metadata": {}, "outputs": [], "source": [ - "! hydromt update wflow wflow_piave_water_demand -i wflow_update_water_demand.yml -d artifact_data -d \"https://github.com/Deltares/hydromt_wflow/releases/download/v0.5.0/wflow_artifacts.yml\" -v" + "! hydromt update wflow wflow_piave_water_demand -i wflow_update_water_demand.yml -d artifact_data -d \"https://github.com/Deltares/hydromt_wflow/releases/download/v0.5.0/wflow_artifacts.yml\" -d ../tests/data/demand/data_catalog.yml -v" ] }, { @@ -110,7 +111,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -397,6 +398,34 @@ "source": [ "Here we see a couple of distinct administrative boundaries (black lines), some of which already follow the catchment boundaries. When the catchment crosses an administrative boundary, that region receives a different but unique identifier. Note that we are using level 2 boundaries here, which might not be the most realistic for a region of this size." ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Surface water frac used for water allocation\n", + "\n", + "By default, Wflow will allocate all water demands with water from the surface water (frac_sw_used=1). However, if in a certain areas, groundwater or other non conventional sources can be used, the frac_sw_used of each allocation area can be reduced and prepared using the **setup_allocation_surfacewaterfrac** method.\n", + "\n", + "Here we used global data from GLOFAS (Lisflood) for the fraction of grounwater used, presence of groundwater bodies and non conventional sources (0 is Piave). The water allocations are the ones we prepared in the previous step in order to match better the wflow model basin delineation.\n", + "\n", + "Let's have a look at the resulting map:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1)\n", + "\n", + "ax.set_title(\"Allocation surface water frac used\")\n", + "\n", + "mod.grid[\"frac_sw_used\"].raster.mask_nodata().plot(ax=ax, add_labels=False)\n", + "mod.geoms[\"basins\"].plot(ax=ax, facecolor=\"none\", edgecolor=\"black\")\n", + "mod.geoms[\"rivers\"].plot(ax=ax)" + ] } ], "metadata": { diff --git a/examples/wflow_update_water_demand.yml b/examples/wflow_update_water_demand.yml index d1d8fd22..66b5a6c7 100644 --- a/examples/wflow_update_water_demand.yml +++ b/examples/wflow_update_water_demand.yml @@ -2,6 +2,11 @@ setup_allocation_areas: min_area: 30 admin_bounds_fn: gadm_level2 +setup_allocation_surfacewaterfrac: + gwfrac_fn: lisflood_gwfrac + gwbodies_fn: lisflood_gwbodies + ncfrac_fn: lisflood_ncfrac + # Update model with GLCNMO landuse data in order to add the paddies setup_lulcmaps_with_paddy: lulc_fn: glcnmo diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index edbfbd47..64fac811 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -3381,47 +3381,57 @@ def setup_allocation_areas( self.set_grid(alloc, name="allocation_areas") # Update the settings toml - self.set_config("input.vertical.waterallocation.areas", "allocation_areas") + self.set_config("input.vertical.allocation.areas", "allocation_areas") # Add alloc to geoms self.set_geoms(alloc.raster.vectorize(), name="allocation_areas") - def setup_surfacewaterfrac( + def setup_allocation_surfacewaterfrac( self, gwfrac_fn: Union[str, xr.DataArray], - gwbodies_fn: Union[str, xr.DataArray], - ncfrac_fn: Union[str, xr.DataArray], - waterareas_fn: Union[str, xr.DataArray], + waterareas_fn: Optional[Union[str, xr.DataArray]] = None, + gwbodies_fn: Optional[Union[str, xr.DataArray]] = None, + ncfrac_fn: Optional[Union[str, xr.DataArray]] = None, + interpolate_nodata: bool = False, ): - """Create water demand surface water consumption fraction. + """Create the fraction of water allocated from surface water. - This fraction entails the division of the water demand between surface water - and ground water. + This fraction entails the division of the water demand between surface water, + ground water (aquifers) and non conventional sources (e.g. desalination plants). The surface water fraction is based on the raw groundwater fraction, if - groudwater bodies are present (these are absent in e.g. mountainous regions), + groundwater bodies are present (these are absent in e.g. mountainous regions), a fraction of water consumed that is obtained by non-conventional means and the water source areas. Non-conventional water could e.g. be water acquired by desalination of ocean or other brackish water. + Adds model layer: + + * **frac_sw_used**: fraction of water allocated from surface water [0-1] + Parameters ---------- gwfrac_fn : Union[str, xr.DataArray] The raw groundwater fraction per grid cell. The values of these cells need to be between 0 and 1. - gwbodies_fn : Union[str, xr.DataArray] - The presence of groundwater bodies per grid cell. The values are ought to - be binary (either 0 or 1). - ncfrac_fn : Union[str, xr.DataArray] - The non-conventional fraction. Same types of values apply as for - `gwfrac_fn`. waterareas_fn : Union[str, xr.DataArray] The areas over which the water has to be distributed. This may either be - a global (or more local map) or it may be defined as `allocation_areas`. - When `allocation_areas` is provided, the source areas created by the - `setup_allocation_areas` will be used for this setup method. + a global (or more local map). If not provided, the source areas created by + the `setup_allocation_areas` will be used. + gwbodies_fn : Union[str, xr.DataArray], optional + The presence of groundwater bodies per grid cell. The values are ought to + be binary (either 0 or 1). If they are not provided, we assume groundwater + bodies are present where gwfrac is more than 0. + ncfrac_fn : Union[str, xr.DataArray], optional + The non-conventional fraction. Same types of values apply as for + `gwfrac_fn`. If not provided, we assume no non-conventional sources are + used. + interpolate_nodata : bool, optional + If True, nodata values in the resulting frac_sw_used map will be linearly + interpolated. Else a default value of 1 will be used for nodata values + (default). """ self.logger.info("Preparing surface water fraction map.") # Load the data @@ -3430,21 +3440,33 @@ def setup_surfacewaterfrac( geom=self.region, buffer=2, ) - gwbodies = self.data_catalog.get_rasterdataset( - gwbodies_fn, - geom=self.region, - buffer=2, - ) - ncfrac = self.data_catalog.get_rasterdataset( - ncfrac_fn, - geom=self.region, - buffer=2, - ) + if gwbodies_fn is not None: + gwbodies = self.data_catalog.get_rasterdataset( + gwbodies_fn, + geom=self.region, + buffer=2, + ) + else: + gwbodies = None + if ncfrac_fn is not None: + ncfrac = self.data_catalog.get_rasterdataset( + ncfrac_fn, + geom=self.region, + buffer=2, + ) + else: + ncfrac = None # check wether to use the models own allocation areas - if waterareas_fn == "allocation_areas": + if waterareas_fn is None: self.logger.info("Using wflow model allocation areas.") - waterareas = self.grid[waterareas_fn] + if "allocation_areas" not in self.grid: + self.logger.error( + "No allocation areas found. Run setup_allocation_areas first " + "or provide a waterareas_fn." + ) + return + waterareas = self.grid["allocation_areas"] else: waterareas = self.data_catalog.get_rasterdataset( waterareas_fn, @@ -3453,22 +3475,23 @@ def setup_surfacewaterfrac( ) # Call the workflow - w_frac = workflows.demand.surfacewaterfrac( - self.grid["wflow_dem"], + w_frac = workflows.demand.surfacewaterfrac_used( gwfrac_raw=gwfrac_raw, + da_like=self.grid["wflow_dem"], + waterareas=waterareas, gwbodies=gwbodies, ncfrac=ncfrac, - waterareas=waterareas, + interpolate=interpolate_nodata, ) # Update the settings toml self.set_config( - "input.vertical.waterallocation.frac_sw_used", - "SurfaceWaterFrac", + "input.vertical.allocation.frac_sw_used", + "frac_sw_used", ) # Set the dataarray to the wflow grid - self.set_grid(w_frac, name="SurfaceWaterFrac") + self.set_grid(w_frac, name="frac_sw_used") def setup_domestic_demand( self, diff --git a/hydromt_wflow/workflows/demand.py b/hydromt_wflow/workflows/demand.py index 94ffee8f..f709497d 100644 --- a/hydromt_wflow/workflows/demand.py +++ b/hydromt_wflow/workflows/demand.py @@ -8,6 +8,7 @@ import numpy as np import xarray as xr from hydromt import raster +from hydromt.workflows.grid import grid_from_constant from scipy.ndimage import convolve logger = logging.getLogger(__name__) @@ -16,7 +17,7 @@ "allocation_areas", "domestic", "other_demand", - "surfacewaterfrac", + "surfacewaterfrac_used", "irrigation", ] @@ -387,28 +388,28 @@ def allocation_areas( return alloc -def surfacewaterfrac( - da_like: xr.DataArray, +def surfacewaterfrac_used( gwfrac_raw: xr.DataArray, - gwbodies: xr.DataArray, - ncfrac: xr.DataArray, + da_like: xr.DataArray, waterareas: xr.DataArray, - interpolate: bool = True, + gwbodies: Optional[xr.DataArray] = None, + ncfrac: Optional[xr.DataArray] = None, + interpolate: bool = False, ) -> xr.DataArray: """Create surface water fraction map. Parameters ---------- - da_like : xr.DataArray - Wflow model like grid. gwfrac_raw : xr.DataArray Raw groundwater fraction data map. - gwbodies : xr.DataArray - Groundwater bodies map. Either 0 or 1. - ncfrac : xr.DataArray - Non-conventional water fraction. + da_like : xr.DataArray + Wflow model like grid. waterareas : xr.DataArray Water source areas. + gwbodies : xr.DataArray, optional + Groundwater bodies map. Either 0 or 1. If None, assumes 1. + ncfrac : xr.DataArray, optional + Non-conventional water fraction. If None, assumes 0. interpolate : bool Interpolate missing data values within wflow model domain. @@ -417,41 +418,55 @@ def surfacewaterfrac( xr.DataArray Surface water fraction. """ - # Mask to int value for later use - da_like = da_like.raster.mask_nodata(-9999) - # Resample the data to model resolution gwfrac_raw_mr = gwfrac_raw.raster.reproject_like( da_like, - method="nearest", - ) - x, y = np.where(~np.isnan(gwfrac_raw_mr)) - gwbodies_mr = gwbodies.raster.reproject_like( - da_like, - method="nearest", - ) - ncfrac_mr = ncfrac.raster.reproject_like( - da_like, - method="nearest", - ) - # Prepare the nodata a little as lisflood data has none set by default - waterareas = waterareas.raster.mask_nodata(-9999) - if "_FillValue" not in waterareas.attrs: - waterareas = waterareas.assign_attrs({"_FillValue": -9999}) - # Reproject - waterareas_mr = waterareas.raster.reproject_like( - da_like, - method="nearest", + method="average", ) - # waterareas_mr = waterareas_mr.raster.mask_nodata(-9999) + # Only reproject waterareas if needed + if not waterareas.raster.identical_grid(da_like): + waterareas_mr = waterareas.raster.reproject_like( + da_like, + method="mode", + ) + else: + waterareas_mr = waterareas # Set nodata values to zeros waterareas_mr = waterareas_mr.where( waterareas_mr != waterareas_mr.raster.nodata, 0, ) + if gwbodies is not None: + gwbodies_mr = gwbodies.raster.reproject_like( + da_like, + method="mode", + ) + else: + gwbodies_mr = grid_from_constant( + grid_like=da_like, + constant=1, + name="gwbodies", + nodata=-9999, + dtype=np.int32, + ) + if ncfrac is not None: + ncfrac_mr = ncfrac.raster.reproject_like( + da_like, + method="average", + ) + else: + ncfrac_mr = grid_from_constant( + grid_like=da_like, + constant=0, + name="ncfrac", + nodata=-9999, + dtype=np.float32, + ) + # Get the fractions based on area count - w_pixels = np.take( + x, y = np.where(~np.isnan(gwfrac_raw_mr)) + gw_pixels = np.take( np.bincount( waterareas_mr.values[x, y], weights=gwbodies_mr.values[x, y], @@ -468,13 +483,12 @@ def surfacewaterfrac( # Determine the groundwater fraction gwfrac_val = np.minimum( - gwfrac_raw_mr.values[x, y] * (a_pixels / (w_pixels + 0.01)), + gwfrac_raw_mr.values[x, y] * (a_pixels / (gw_pixels + 0.01)), 1 - ncfrac_mr.values[x, y], ) gwfrac_val[np.where(gwbodies_mr.values[x, y] == 0)] = 0 # invert to get surface water frac - # TODO fix with line from listflood - gwfrac_val = 1 - gwfrac_val + swfrac_val = np.maximum(np.minimum(1 - gwfrac_val - ncfrac_mr.values[x, y], 1), 0) # create the dataarray for the fraction swfrac = xr.full_like( @@ -482,17 +496,18 @@ def surfacewaterfrac( fill_value=np.nan, dtype=np.float32, ).load() - swfrac.name = "SurfaceWaterFrac" + swfrac.name = "frac_sw_used" swfrac.attrs = {"_FillValue": -9999} - swfrac = swfrac.copy() + swfrac.raster.set_crs(da_like.raster.crs) # Set and interpolate the values - swfrac.values[x, y] = gwfrac_val + swfrac = swfrac.copy() # to avoid read-only error + swfrac.values[x, y] = swfrac_val if interpolate: - swfrac = swfrac.interpolate_na(dim=swfrac.raster.x_dim, method="linear") - swfrac = swfrac.interpolate_na( - dim=swfrac.raster.x_dim, method="linear", fill_value="extrapolate" - ) + swfrac = swfrac.raster.interpolate_na(method="linear", extrapolate=True) + else: + # Fill na with 1 (no groundwater or nc sources used) + swfrac = swfrac.fillna(1.0) # Set the nodata values based on the dem of the model (da_like) swfrac = swfrac.where(da_like != da_like.raster.nodata, -9999) diff --git a/tests/data/demand/data_catalog.yml b/tests/data/demand/data_catalog.yml new file mode 100644 index 00000000..e9bb80c7 --- /dev/null +++ b/tests/data/demand/data_catalog.yml @@ -0,0 +1,33 @@ +lisflood_gwbodies: + data_type: RasterDataset + path: lisflood_gwbodies.tif + driver: raster + filesystem: local + meta: + source_url: https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/CEMS-GLOFAS/LISFLOOD_static_and_parameter_maps_for_GloFAS/Water_demand/ + crs: 4326 +lisflood_gwfrac: + data_type: RasterDataset + path: lisflood_gwfrac.tif + driver: raster + filesystem: local + meta: + source_url: https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/CEMS-GLOFAS/LISFLOOD_static_and_parameter_maps_for_GloFAS/Water_demand/ + crs: 4326 +lisflood_ncfrac: + data_type: RasterDataset + path: lisflood_ncfrac.tif + driver: raster + filesystem: local + meta: + source_url: https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/CEMS-GLOFAS/LISFLOOD_static_and_parameter_maps_for_GloFAS/Water_demand/ + crs: 4326 +lisflood_waterregions: + data_type: RasterDataset + path: lisflood_waterregions.tif + driver: raster + filesystem: local + nodata: -9999 + meta: + source_url: https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/CEMS-GLOFAS/LISFLOOD_static_and_parameter_maps_for_GloFAS/Water_demand/ + crs: 4326 diff --git a/tests/data/demand/lisflood_gwbodies.tif b/tests/data/demand/lisflood_gwbodies.tif new file mode 100644 index 00000000..35d5daa8 Binary files /dev/null and b/tests/data/demand/lisflood_gwbodies.tif differ diff --git a/tests/data/demand/lisflood_gwfrac.tif b/tests/data/demand/lisflood_gwfrac.tif new file mode 100644 index 00000000..73f9c763 Binary files /dev/null and b/tests/data/demand/lisflood_gwfrac.tif differ diff --git a/tests/data/demand/lisflood_ncfrac.tif b/tests/data/demand/lisflood_ncfrac.tif new file mode 100644 index 00000000..85e17d6a Binary files /dev/null and b/tests/data/demand/lisflood_ncfrac.tif differ diff --git a/tests/data/demand/lisflood_waterregions.tif b/tests/data/demand/lisflood_waterregions.tif new file mode 100644 index 00000000..14e021d7 Binary files /dev/null and b/tests/data/demand/lisflood_waterregions.tif differ diff --git a/tests/test_model_methods.py b/tests/test_model_methods.py index 6ddc658e..d3ae24c4 100644 --- a/tests/test_model_methods.py +++ b/tests/test_model_methods.py @@ -926,6 +926,58 @@ def test_setup_allocation_areas(example_wflow_model, tmpdir): assert np.all(np.sort(uni) == [36, 37, 57]) +def test_setup_allocation_surfacewaterfrac(example_wflow_model, tmpdir): + # Read the data and set new root + example_wflow_model.read() + example_wflow_model.set_root( + Path( + tmpdir, + ), + mode="w", + ) + # Add lisflood data from test + lisflood_yml = join(TESTDATADIR, "demand", "data_catalog.yml") + example_wflow_model.data_catalog = example_wflow_model.data_catalog.from_yml( + lisflood_yml + ) + + # Use the method fully with lisflood + example_wflow_model.setup_allocation_surfacewaterfrac( + gwfrac_fn="lisflood_gwfrac", + waterareas_fn="lisflood_waterregions", + gwbodies_fn="lisflood_gwbodies", + ncfrac_fn="lisflood_ncfrac", + interpolate_nodata=False, + ) + + # Assert entries + assert "frac_sw_used" in example_wflow_model.grid + assert np.isclose( + example_wflow_model.grid["frac_sw_used"].raster.mask_nodata().mean().values, + 0.9411998, + ) + + # Use the method without gwbodies and ncfrac and waterareas from wflow + example_wflow_model.setup_allocation_areas( + admin_bounds_fn="gadm_level2", + min_area=30, + ) + example_wflow_model.setup_allocation_surfacewaterfrac( + gwfrac_fn="lisflood_gwfrac", + waterareas_fn=None, + gwbodies_fn=None, + ncfrac_fn=None, + interpolate_nodata=False, + ) + + # Assert entries + assert "frac_sw_used" in example_wflow_model.grid + assert np.isclose( + example_wflow_model.grid["frac_sw_used"].raster.mask_nodata().mean().values, + 0.9865037, + ) + + def test_setup_non_irrigation(example_wflow_model, tmpdir): # Read the data example_wflow_model.read()