Water demand #226
Water demand #226
Conversation
code also useful for rootingdepth, but should require additional information before being useful for wflow
|
Thanks for adding irrigation and crop factor maps @JoostBuitink ! For a test version I think it would be good to also add a soil profile for irrigated paddy areas: |
There was a problem hiding this comment.
Nice developments @dalmijn and @JoostBuitink ! And a good start to get demand data from pcr_globwb or irrigation/paddy areas!
I reviewed mainly the functions from @dalmijn, as I don't think setup_irrigation is fully done yet.
There's a few comments and suggestions here and there. Maybe one function I did not fully understand is how the allocation areas are derived so would be happy to discuss further!
hydromt_wflow/wflow.py
Outdated
| def setup_allocation( | ||
| self, | ||
| min_area: float | int = 0, | ||
| admin_bounds_fn: str = "gadm", |
There was a problem hiding this comment.
We are cleaning this up still for wflow but we do not want to assume a default dataset value anymore, they are just mandatory arguments :) the example from artifact_data can be mentioned in the examples. Which might be good to add here! (update water demands example to showcase the new methods)
hydromt_wflow/wflow.py
Outdated
| self, | ||
| min_area: float | int = 0, | ||
| admin_bounds_fn: str = "gadm", | ||
| admin_level: int = None, |
There was a problem hiding this comment.
This is too specific to gadm and people should be able to use their own geometry. I think in the data catalog, it's actually gadm_level1 that you can call so just use this in admin_bounds_fn (in testing / examples).
Possibly you can add **kwargs here to the get_geodataframe method so that a user is able to select a specific layer in their data via the driver_kwargs (if they did not do that in their own data catalog)
hydromt_wflow/wflow.py
Outdated
| @@ -2494,6 +2498,292 @@ def setup_rootzoneclim( | |||
| # update config | |||
| self.set_config("input.vertical.rootingdepth", update_toml_rootingdepth) | |||
|
|
|||
| def setup_allocation( | |||
There was a problem hiding this comment.
| def setup_allocation( | |
| def setup_allocation_areas( |
To match better wflow names. Makes it also clearer for me that you are only preparing areas here not demand or actual allocation data.
hydromt_wflow/wflow.py
Outdated
|
|
||
| def setup_non_irigation( | ||
| self, | ||
| non_irigation_fn: str = "pcr_globwb", |
There was a problem hiding this comment.
Either None if optional or required argument
| non_irigation_fn: str = "pcr_globwb", | |
| non_irigation_fn: Union[str, Path, xr.Dataset], |
hydromt_wflow/workflows/demand.py
Outdated
| sub_basins = sub_basins.dissolve("uid", sort=False, as_index=False) | ||
| _count += 1 | ||
|
|
||
| alloc = full_like(ds_like["dem_subgrid"], nodata=-9999, lazy=True).astype(int) |
There was a problem hiding this comment.
Why dem_subgrid? If kinematic wave models, I'm not event sure the map is available... "wflow_dem" or "wflow_ldd" or "wflow_subcatch" are safer as they should always be there (especially wflow_subcatch that you should actually use rather than self.basins actually as self.basins is unfortunately basins at merit_hydro resolution rather than wflow model resolution...)
There was a problem hiding this comment.
Updated the code to no longer rely on dem_subgrid. Other remark about self.basins vs. wflow_subcatch is still open I guess
hydromt_wflow/workflows/demand.py
Outdated
| ds_like.raster.bounds, | ||
| ) | ||
|
|
||
| # Create dummy dataset from this info |
There was a problem hiding this comment.
I guess this is to ensure better mass conservation? Ie you are not creating or removing m3 just because or up or downscaling? Then I think it's nice to apply to not just domestic demand.
If you make a nice enough I would also move this to hydromt core. It's useful to have this kind of mass conservation reprojection for water demands or population or livestock or any count data :)
There was a problem hiding this comment.
Not quite sure anymore what this is referring to.
some todo's still open, but contains all required functionality
They both point to the same map currently
|
ToDo's still open:
|
There was a problem hiding this comment.
First comments from the review, mainly for allocation_areas and surfacewaterfrac and docs. I will continue tomorrow or later this week for the other functions.
I also ran the example notebook and it looks really good!
I have one comment already on setup_irrigation from that notebook: to simplify inputs, I think it would be good to have consistent maps for the irrigated areas. I assume anyway wflow internally will not apply crop_factor or irrigation_trigger for non irrigated lands but would be good to also then get consistent maps from hydromt on this to better understand what really happens in Wflow.jl
-
crop_factor: I would assume only valid for irrigation areas so should be 1 if no irrigated areas (paddy or non paddy)
-
Same for irrigation_trigger, I would assume always 0 for non-irrigation areas
docs/changelog.rst
Outdated
There was a problem hiding this comment.
Changes from this PR should move to the unreleased section
docs/user_guide/wflow_update.rst
Outdated
There was a problem hiding this comment.
Also need to update the table in docs/user_guide/wflow_model_setup.rst to include the new setup methods
| # Get the fractions based on area count | ||
| w_pixels = np.take( | ||
| np.bincount( | ||
| waterareas_mr.values[x, y], | ||
| weights=gwbodies_mr.values[x, y], | ||
| ), | ||
| waterareas_mr.values[x, y], | ||
| ) |
There was a problem hiding this comment.
I'm a little reluctant to use numpy instead of xarray and if I understand bincount well if you have poorly named waterareas you may get bins from a very wide range of numbers (eg two basins with id 1 and 1000 or more).
Could you maybe use a groupby and count? https://docs.xarray.dev/en/stable/generated/xarray.DataArray.groupby.html
If too complex this is okay for now
| swfrac = xr.full_like( | ||
| gwfrac_raw_mr, | ||
| fill_value=np.nan, | ||
| dtype=np.float32, | ||
| ).load() | ||
| swfrac.name = "SurfaceWaterFrac" | ||
| swfrac.attrs = {"_FillValue": -9999} | ||
| swfrac = swfrac.copy() |
There was a problem hiding this comment.
maybe you can use hydromt.workflows.grid.grid_from_constant? Also good to have the crs property in there
| 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" | ||
| ) |
There was a problem hiding this comment.
I think we discussed this but why not hydromt.raster.interpolate_na with extrapolate=True?
Or else use keep_attrs=True else you will loose again nodata and crs attributes
There was a problem hiding this comment.
It had a more specific use case in that instance we discussed, Ill look at it for this instance.
hydromt_wflow/workflows/demand.py
Outdated
| gwfrac_val = np.minimum( | ||
| gwfrac_raw_mr.values[x, y] * (a_pixels / (w_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 | ||
| gwfrac_val = 1 - gwfrac_val |
There was a problem hiding this comment.
Maybe I need to read the paper but I would say sw_frac = 1 - gwfrac - ncfrac no?
But here let's say gwfrac is 0.4 and ncfrac 0.2 then swfrac is 0.6 and not 0.4?
hydromt_wflow/wflow.py
Outdated
| if "time" in data.dims: | ||
| tname = "time" | ||
| time_axes = { | ||
| k: v for k, v in dict(self.grid.dims).items() if k.startswith("time") | ||
| } | ||
| if data["time"].size not in time_axes.values(): | ||
| tname = f"time_{data['time'].size}" if "time" in time_axes else tname | ||
| else: | ||
| k = list( | ||
| filter(lambda x: time_axes[x] == data["time"].size, time_axes) | ||
| )[0] | ||
| tname = k | ||
|
|
||
| if tname != "time": | ||
| data = data.rename_dims({"time": tname}) |
There was a problem hiding this comment.
What is this??
Does Wflow supports new types of cyclic inputs?
Is this really needed or can't we standardize in the workflows?
There was a problem hiding this comment.
Well review is complete!
This is a massive amount of work and very important so well done to you both :)
I (of course) have some comments from the "hydromt" and user side so let's discuss the main feedbacks from the review together.
A recap here of the main TODOS / things to discuss:
- Update changelog
- Update setup methods table for wflow in the docs
- Replace key
[input.vertical.waterallocation]with[input.vertical.allocation] - Add setup_surfacewater to example
- Use of _MAPS dictionnary
- All setup methods: document required variables and units in the input datasets
- Setup_irrigation: Vertical.paddy h_min, h_opt, h_max do not appear in toml
- Discuss renaming of the following methods:
setup_allocation_areas,setup_surfacewaterfrac,setup_non_irrigation - Discuss the implementation of
setup_allocation_areas: intersection is good but would like to discuss deriving the subbasins after that - Discuss the implementation of
setup_surfacewaterfrac: default values and mistake in the calculation method? + take out the preparation of the allocation areas map? - Discuss the implementation of
setup_non_irrigation: resampling method and use of an extra low res grid (res_factor !=0) - Discuss the implementation of
setup_irrigation: Split for paddy and non-paddy? Simplify arguments (and code)? Harmonise crop_factor and irrigation_trigger - Discuss the implementation of
set_grid: multiple time dims in wflow staticmaps? - Optional ie something from the Wflow.jl side or just a new hydromt feature that becomes more than nice to have for v1 of wflow: new TOML arguments are even more "subdivided"
Finally, when all comments are addressed:
- Do a check that the produced model does indeed run with the latest version of the wflow code
| _count = 0 | ||
|
|
||
| # Create touched and not touched by rivers datasets | ||
| while True: |
There was a problem hiding this comment.
And why not the outflow_idxs method above instead of your long while loop? You have to trust the user to use an uparea threshold more than their river but it's already the case for the setup_floodplains method so we can actually do a check on this in the setup method first :)
Unless you would like to allow for exchanges between subbasins in terms of water allocation areas? I guess this would mainly be for small coastal basins as in your example above?
Then you can maybe decide to interpolate/extrapolate with nearest neighbor?
| nodata = -9999 | ||
| # Split based on admin bounds | ||
| sub_basins = basins.copy() | ||
| sub_basins["uid"] = range(len(sub_basins)) |
There was a problem hiding this comment.
I'm not sure how well 0 is supported, you may want to use arange to start at 1
Co-authored-by: hboisgon <[email protected]>
Co-authored-by: hboisgon <[email protected]>
Co-authored-by: hboisgon <[email protected]>
Now includes fix when layers have a different dimension, and includes a better docstring
Issue addressed
Fixes #87
Explanation
Explain how you addressed the bug/feature request, what choices you made and why.
Checklist
mainAdditional Notes (optional)
Add any additional notes or information that may be helpful.