-
Notifications
You must be signed in to change notification settings - Fork 13
/
wflow.py
3146 lines (2823 loc) · 138 KB
/
wflow.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""Implement wflow model class"""
# Implement model class following model API
import os
from os.path import join, dirname, basename, isfile, isdir
from pathlib import Path
from typing import Union, Optional, List
import glob
import numpy as np
import pandas as pd
import geopandas as gpd
import xarray as xr
from shapely.geometry import box
import pyproj
import toml
import codecs
import pyflwdir
from pyflwdir import core_d8, core_ldd, core_conversion
from dask.diagnostics import ProgressBar
import logging
import rioxarray # required for rio accessor
# from dask.distributed import LocalCluster, Client, performance_report
import hydromt
from hydromt.models.model_api import Model
from hydromt import flw
from hydromt.io import open_mfraster
from . import utils, workflows, DATADIR
__all__ = ["WflowModel"]
logger = logging.getLogger(__name__)
# specify pcraster map types
# NOTE non scalar (float) data types only
PCR_VS_MAP = {
"wflow_ldd": "ldd",
"wflow_river": "bool",
"wflow_streamorder": "ordinal",
"wflow_gauges": "nominal", # to avoid large memory usage in pcraster.aguila
"wflow_subcatch": "nominal", # idem.
"wflow_landuse": "nominal",
"wflow_soil": "nominal",
"wflow_reservoirareas": "nominal",
"wflow_reservoirlocs": "nominal",
"wflow_lakeareas": "nominal",
"wflow_lakelocs": "nominal",
"wflow_glacierareas": "nominal",
}
class WflowModel(Model):
"""This is the wflow model class"""
_NAME = "wflow"
_CONF = "wflow_sbm.toml"
_DATADIR = DATADIR
_GEOMS = {}
_MAPS = {
"flwdir": "wflow_ldd",
"elevtn": "wflow_dem",
"subelv": "dem_subgrid",
"uparea": "wflow_uparea",
"strord": "wflow_streamorder",
"basins": "wflow_subcatch",
"rivlen": "wflow_riverlength",
"rivmsk": "wflow_river",
"rivwth": "wflow_riverwidth",
"lndslp": "Slope",
"rivslp": "RiverSlope",
"rivdph": "RiverDepth",
"rivman": "N_River",
"gauges": "wflow_gauges",
"landuse": "wflow_landuse",
"resareas": "wflow_reservoirareas",
"reslocs": "wflow_reservoirlocs",
"lakeareas": "wflow_lakeareas",
"lakelocs": "wflow_lakelocs",
"glacareas": "wflow_glacierareas",
"glacfracs": "wflow_glacierfrac",
"glacstore": "wflow_glacierstore",
}
_FOLDERS = [
"staticgeoms",
"instate",
"run_default",
]
_CATALOGS = join(_DATADIR, "parameters_data.yml")
def __init__(
self,
root=None,
mode="w",
config_fn=None,
data_libs=None,
deltares_data=False,
logger=logger,
):
super().__init__(
root=root,
mode=mode,
config_fn=config_fn,
data_libs=data_libs,
deltares_data=deltares_data,
logger=logger,
)
# wflow specific
self._intbl = dict()
self._tables = dict()
self._flwdir = None
self.data_catalog.from_yml(self._CATALOGS)
# COMPONENTS
def setup_basemaps(
self,
region,
res=1 / 120.0,
hydrography_fn="merit_hydro",
basin_index_fn="merit_hydro_index",
upscale_method="ihu",
):
"""
This component sets the ``region`` of interest and ``res`` (resolution in degrees) of the
model. All DEM and flow direction related maps are then build.
If the model resolution is larger than the source data resolution,
the flow direction is upscaled using the ``upscale_method``, by default the
Iterative Hydrography Upscaling (IHU).
The default ``hydrography_fn`` is "merit_hydro" (`MERIT hydro <http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro/index.html>`_
at 3 arcsec resolution) Alternative sources include "merit_hydro_1k" at 30 arcsec resolution.
Users can also supply their own elevation and flow direction data in any CRS and not only EPSG:4326.
Note that in order to define the region, using points or bounding box, the coordinates of the points / bounding box
should be in the same CRS than the hydrography data. The wflow model will then also be in the same CRS than the
hydrography data in order to avoid assumptions and reprojection errors. If the user wishes to use a different CRS,
we recommend first to reproject the hydrography data seperately because calling hydromt build. You can find examples
on how to reproject or prepare hydrography data in the
`prepare flow directions example notebok <https://deltares.github.io/hydromt_wflow/latest/_examples/prepare_ldd.html>`.
Adds model layers:
* **wflow_ldd** map: flow direction in LDD format [-]
* **wflow_subcatch** map: basin ID map [-]
* **wflow_uparea** map: upstream area [km2]
* **wflow_streamorder** map: Strahler stream order [-]
* **wflow_dem** map: average elevation [m+REF]
* **dem_subgrid** map: subgrid outlet elevation [m+REF]
* **Slope** map: average land surface slope [m/m]
* **basins** geom: basins boundary vector
* **region** geom: region boundary vector
Parameters
----------
hydrography_fn : str
Name of data source for basemap parameters.
* Required variables: ['flwdir', 'uparea', 'basins', 'strord', 'elevtn']
* Optional variables: ['lndslp', 'mask']
basin_index_fn : str
Name of data source for basin_index data linked to hydrography_fn.
region : dict
Dictionary describing region of interest.
See :py:function:~basin_mask.parse_region for all options
res : float
Output model resolution
upscale_method : {'ihu', 'eam', 'dmm'}
Upscaling method for flow direction data, by default 'ihu'.
See Also
--------
hydromt.workflows.parse_region
hydromt.workflows.get_basin_geometry
workflows.hydrography
workflows.topography
"""
self.logger.info(f"Preparing base hydrography basemaps.")
# retrieve global data (lazy!)
ds_org = self.data_catalog.get_rasterdataset(hydrography_fn)
# get basin geometry and clip data
kind, region = hydromt.workflows.parse_region(region, logger=self.logger)
xy = None
if kind in ["basin", "subbasin", "outlet"]:
if basin_index_fn is not None:
bas_index = self.data_catalog[basin_index_fn]
else:
bas_index = None
geom, xy = hydromt.workflows.get_basin_geometry(
ds=ds_org,
kind=kind,
basin_index=bas_index,
logger=self.logger,
**region,
)
elif "bbox" in region:
bbox = region.get("bbox")
geom = gpd.GeoDataFrame(geometry=[box(*bbox)], crs=4326)
elif "geom" in region:
geom = region.get("geom")
else:
raise ValueError(f"wflow region argument not understood: {region}")
if geom is not None and geom.crs is None:
raise ValueError("wflow region geometry has no CRS")
ds_org = ds_org.raster.clip_geom(geom, align=res, buffer=10)
ds_org.coords["mask"] = ds_org.raster.geometry_mask(geom)
self.logger.debug("Adding basins vector to staticgeoms.")
self.set_staticgeoms(geom, name="basins")
# Check on resolution (degree vs meter) depending on ds_org res/crs
scale_ratio = int(np.round(res / ds_org.raster.res[0]))
if scale_ratio < 1:
raise ValueError(
f"The model resolution {res} should be larger than the {hydrography_fn} resolution {ds_org.raster.res[0]}"
)
if ds_org.raster.crs.is_geographic:
if res > 1: # 111 km
raise ValueError(
f"The model resolution {res} should be smaller than 1 degree (111km) for geographic coordinate systems. "
"Make sure you provided res in degree rather than in meters."
)
# setup hydrography maps and set staticmap attribute with renamed maps
ds_base, _ = workflows.hydrography(
ds=ds_org,
res=res,
xy=xy,
upscale_method=upscale_method,
logger=self.logger,
)
# Convert flow direction from d8 to ldd format
flwdir_data = ds_base["flwdir"].values.astype(np.uint8) # force dtype
# if d8 convert to ldd
if core_d8.isvalid(flwdir_data):
data = core_conversion.d8_to_ldd(flwdir_data)
da_flwdir = xr.DataArray(
name="flwdir",
data=data,
coords=ds_base.raster.coords,
dims=ds_base.raster.dims,
attrs=dict(
long_name="ldd flow direction",
_FillValue=core_ldd._mv,
),
)
ds_base["flwdir"] = da_flwdir
# Rename and add to staticmaps
rmdict = {k: v for k, v in self._MAPS.items() if k in ds_base.data_vars}
self.set_staticmaps(ds_base.rename(rmdict))
# setup topography maps
ds_topo = workflows.topography(
ds=ds_org, ds_like=self.staticmaps, method="average", logger=self.logger
)
ds_topo["lndslp"] = np.maximum(ds_topo["lndslp"], 0.0)
rmdict = {k: v for k, v in self._MAPS.items() if k in ds_topo.data_vars}
self.set_staticmaps(ds_topo.rename(rmdict))
# set basin geometry
self.logger.debug(f"Adding region vector to staticgeoms.")
self.set_staticgeoms(self.region, name="region")
# update toml for degree/meters if needed
if ds_base.raster.crs.is_projected:
self.set_config("model.sizeinmetres", True)
def setup_rivers(
self,
hydrography_fn,
river_geom_fn=None,
river_upa: float = 30,
rivdph_method: str = "powlaw",
slope_len: float = 2e3,
min_rivlen_ratio=0.0,
min_rivdph: float = 1,
min_rivwth: float = 30,
smooth_len: float = 5e3,
rivman_mapping_fn: str = "roughness_river_mapping_default",
elevtn_map: str = "wflow_dem",
river_routing: str = "kinematic-wave",
connectivity: int = 8,
**kwargs,
):
"""
This component sets the all river parameter maps.
The river mask is defined by all cells with a mimimum upstream area threshold
``river_upa`` [km2].
The river length is defined as the distance from the subgrid outlet pixel to
the next upstream subgrid outlet pixel. The ``min_rivlen_ratio`` is the minimum
global river length to avg. cell resolution ratio and is used as a threshold in
window based smoothing of river length.
The river slope is derived from the subgrid elevation difference between pixels at a
half distance ``slope_len`` [m] up- and downstream from the subgrid outlet pixel.
The river manning roughness coefficient is derived based on reclassification
of the streamorder map using a lookup table ``rivman_mapping_fn``.
The river width is derived from the nearest river segment in ``river_geom_fn``.
Data gaps are filled by the nearest valid upstream value and averaged along
the flow directions over a length ``smooth_len`` [m]
The river depth is calculated using the ``rivdph_method``, by default powlaw:
h = hc*Qbf**hp, which is based on qbankfull discharge from the nearest river
segment in ``river_geom_fn`` and takes optional arguments for the hc
(default = 0.27) and hp (default = 0.30) parameters. For other methods see
:py:meth:`hydromt.workflows.river_depth`.
If ``river_routing`` is set to "local-inertial", the bankfull elevantion map can be
conditioned based on the average cell elevation ("wflow_dem") or subgrid outlet pixel
elevation ("dem_subgrid"). The subgrid elevation might provide a better representation
of the river elevation profile, however in combination with local-inertial land routing
(see :py:meth:`setup_floodplains`) the subgrid elevation will likely overestimate the
floodplain storage capacity. Note that the same input elevation map should be used for
river bankfull elevation and land elevation when using local-inertial land routing.
Adds model layers:
* **wflow_river** map: river mask [-]
* **wflow_riverlength** map: river length [m]
* **wflow_riverwidth** map: river width [m]
* **RiverDepth** map: bankfull river depth [m]
* **RiverSlope** map: river slope [m/m]
* **N_River** map: Manning coefficient for river cells [s.m^1/3]
* **rivers** geom: river vector based on wflow_river mask
* **hydrodem** map: hydrologically conditioned elevation [m+REF]
Parameters
----------
hydrography_fn : str, Path
Name of data source for hydrography data.
Must be same as setup_basemaps for consistent results.
* Required variables: ['flwdir', 'uparea', 'elevtn']
* Optional variables: ['rivwth', 'qbankfull']
river_geom_fn : str, Path, optional
Name of data source for river data.
* Required variables: ['rivwth', 'qbankfull']
river_upa : float
minimum upstream area threshold for the river map [km2]
slope_len : float
length over which the river slope is calculated [km]
min_rivlen_ratio: float
Ratio of cell resolution used minimum length threshold in a moving
window based smoothing of river length, by default 0.0
The river length smoothing is skipped if min_riverlen_ratio = 0.
For details about the river length smoothing, see :py:meth:`pyflwdir.FlwdirRaster.smooth_rivlen`
rivdph_method : {'gvf', 'manning', 'powlaw'}
see py:meth:`hydromt.workflows.river_depth` for details, by default "powlaw"
river_routing : {'kinematic-wave', 'local-inertial'}
Routing methodology to be used, by default "kinematic-wave".
smooth_len : float, optional
Length [m] over which to smooth the output river width and depth, by default 5e3
min_rivdph : float, optional
Minimum river depth [m], by default 1.0
min_rivwth : float, optional
Minimum river width [m], by default 30.0
elevtn_map : str, optional
by default "wflow_dem",
See Also
--------
workflows.river_bathymetry
hydromt.workflows.river_depth
pyflwdir.FlwdirRaster.river_depth
setup_floodplains
"""
self.logger.info(f"Preparing river maps.")
rivdph_methods = ["gvf", "manning", "powlaw"]
if rivdph_method not in rivdph_methods:
raise ValueError(f'"{rivdph_method}" unknown. Select from {rivdph_methods}')
routing_options = ["kinematic-wave", "local-inertial"]
if river_routing not in routing_options:
raise ValueError(
f'river_routing="{river_routing}" unknown. Select from {routing_options}.'
)
# read data
ds_hydro = self.data_catalog.get_rasterdataset(
hydrography_fn, geom=self.region, buffer=10
)
ds_hydro.coords["mask"] = ds_hydro.raster.geometry_mask(self.region)
# get rivmsk, rivlen, rivslp
# read model maps and revert wflow to hydromt map names
inv_rename = {v: k for k, v in self._MAPS.items() if v in self.staticmaps}
ds_riv = workflows.river(
ds=ds_hydro,
ds_model=self.staticmaps.rename(inv_rename),
river_upa=river_upa,
slope_len=slope_len,
channel_dir="up",
min_rivlen_ratio=min_rivlen_ratio,
logger=self.logger,
)[0]
ds_riv["rivmsk"] = ds_riv["rivmsk"].assign_attrs(
river_upa=river_upa, slope_len=slope_len, min_rivlen_ratio=min_rivlen_ratio
)
dvars = ["rivmsk", "rivlen", "rivslp"]
rmdict = {k: self._MAPS.get(k, k) for k in dvars}
self.set_staticmaps(ds_riv[dvars].rename(rmdict))
# TODO make separate workflows.river_manning method
# Make N_River map from csv file with mapping between streamorder and N_River value
strord = self.staticmaps[self._MAPS["strord"]].copy()
df = self.data_catalog.get_dataframe(rivman_mapping_fn)
# max streamorder value above which values get the same N_River value
max_str = df.index[-2]
# if streamorder value larger than max_str, assign last value
strord = strord.where(strord <= max_str, max_str)
# handle missing value (last row of csv is mapping of missing values)
strord = strord.where(strord != strord.raster.nodata, -999)
strord.raster.set_nodata(-999)
ds_nriver = workflows.landuse(
da=strord,
ds_like=self.staticmaps,
df=df,
logger=self.logger,
)
self.set_staticmaps(ds_nriver)
# get rivdph, rivwth
# while we still have setup_riverwidth one can skip river_bathymetry here
# TODO make river_geom_fn required after removing setup_riverwidth
if river_geom_fn is not None:
gdf_riv = self.data_catalog.get_geodataframe(
river_geom_fn, geom=self.region
)
# reread model data to get river maps
inv_rename = {v: k for k, v in self._MAPS.items() if v in self.staticmaps}
ds_riv1 = workflows.river_bathymetry(
ds_model=self.staticmaps.rename(inv_rename),
gdf_riv=gdf_riv,
method=rivdph_method,
smooth_len=smooth_len,
min_rivdph=min_rivdph,
min_rivwth=min_rivwth,
logger=self.logger,
**kwargs,
)
rmdict = {k: v for k, v in self._MAPS.items() if k in ds_riv1.data_vars}
self.set_staticmaps(ds_riv1.rename(rmdict))
# update config
self.set_config("input.lateral.river.bankfull_depth", self._MAPS["rivdph"])
self.logger.debug(f"Adding rivers vector to staticgeoms.")
self.staticgeoms.pop("rivers", None) # remove old rivers if in staticgeoms
self.rivers # add new rivers to staticgeoms
# Add hydrologically conditioned elevation map for the river, if required
if river_routing == "local-inertial":
postfix = {"wflow_dem": "_avg", "dem_subgrid": "_subgrid"}.get(
elevtn_map, ""
)
name = f"hydrodem{postfix}"
ds_out = flw.dem_adjust(
da_flwdir=self.staticmaps[self._MAPS["flwdir"]],
da_elevtn=self.staticmaps[elevtn_map],
da_rivmsk=self.staticmaps[self._MAPS["rivmsk"]],
flwdir=self.flwdir,
connectivity=connectivity,
river_d8=True,
logger=self.logger,
).rename(name)
self.set_staticmaps(ds_out)
# update toml model.river_routing
self.logger.debug(
f'Update wflow config model.river_routing="{river_routing}"'
)
self.set_config("model.river_routing", river_routing)
self.set_config("input.lateral.river.bankfull_depth", self._MAPS["rivdph"])
self.set_config("input.lateral.river.bankfull_elevation", name)
else:
self.set_config("model.river_routing", river_routing)
def setup_floodplains(
self,
hydrography_fn,
floodplain_type: str,
### Options for 1D floodplains
river_upa: Optional[float] = None,
flood_depths: List = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0],
### Options for 2D floodplains
elevtn_map: str = "wflow_dem",
connectivity: int = 4,
):
"""
This components adds floodplain information to the model schematistation. The user can
define what type of floodplains are required (1D or 2D), through the ``floodplain_type``
argument.
If ``floodplain_type`` is set to "1d", a floodplain profile is derived for every river
cell. It adds a map with floodplain volume per flood depth, which is used in the wflow
1D floodplain schematisation.
Note, it is important to use the same river uparea value as used in the :py:meth:`setup_rivers`
method.
If ``floodplain_type`` is set to "2d", this component adds a hydrologically conditioned
elevation (hydrodem) map for land routing (local-inertial). For this options, landcells
need to be conditioned to D4 flow directions otherwise pits may remain in the land
cells.
The conditioned elevation can be based on the average cell elevation ("wflow_dem") or
subgrid outlet pixel elevation ("dem_subgrid"). Note that the subgrid elevation will
likely overestimate the floodplain storage capacity.
Additionally, note that the same input elevation map should be used for river bankfull
elevation and land elevation when using local-inertial land routing.
Requires :py:meth:`setup_rivers` to be executed beforehand (with ``river_routing`` set to
"local-inertial").
Adds model layers:
* **floodplain_volume** map: map with floodplain volumes, has flood depth as third
dimension [m3] (for 1D floodplains)
* **hydrodem** map: hydrologically conditioned elevation [m+REF] (for 2D floodplains)
Parameters
----------
floodplain_type: {"1d", "2d"}
Option defining the type of floodplains, see below what arguments are related to
the different floodplain types
hydrography_fn : str, Path
Name of data source for hydrography data. Must be same as setup_basemaps for
consistent results.
* Required variables: ['flwdir', 'uparea', 'elevtn']
river_upa : float, optional
(1D floodplains) minimum upstream area threshold for drain in the HAND. Optional
value, as it is inferred from the staticmaps metadata, to be consistent with
setup_rivers.
flood_depths : tuple of float, optional
(1D floodplains) flood depths at which a volume is derived, by default
[0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0]
elevtn_map: {"wflow_dem", "dem_subgrid"}
(2D floodplains) Name of staticmap to hydrologically condition, by default
"wflow_dem"
See Also
--------
workflows.river_floodplain_volume
hydromt.flw.dem_adjust
pyflwdir.FlwdirRaster.dem_adjust
setup_rivers
"""
if self.get_config("model.river_routing") != "local-inertial":
raise ValueError(
f"Floodplains (1d or 2d) are currently only supported with local intertial river routing"
)
r_list = ["1d", "2d"]
if floodplain_type not in r_list:
raise ValueError(
f'river_routing="{floodplain_type}" unknown. Select from {r_list}.'
)
# Adjust settings based on floodplain_type selection
if floodplain_type == "1d":
floodplain_1d = True
land_routing = "kinematic-wave"
if not hasattr(pyflwdir.FlwdirRaster, "ucat_volume"):
self.logger.warning("This method requires pyflwdir >= 0.5.6")
return
self.logger.info("Preparing 1D river floodplain_volume map.")
# read data
ds_hydro = self.data_catalog.get_rasterdataset(
hydrography_fn, geom=self.region, buffer=10
)
ds_hydro.coords["mask"] = ds_hydro.raster.geometry_mask(self.region)
# try to get river uparea from staticmaps, throw error if not specified or when found but different from specified value
new_river_upa = self.staticmaps[self._MAPS["rivmsk"]].attrs.get(
"river_upa", river_upa
)
if new_river_upa is None:
raise ValueError(
"No value for `river_upa` specified, and the value cannot be inferred from the staticmaps attributes"
)
elif new_river_upa != river_upa and river_upa is not None:
raise ValueError(
f"Value specified for river_upa ({river_upa}) is different from the value found in the staticmaps ({new_river_upa})"
)
self.logger.debug(f"Using river_upa value value of: {new_river_upa}")
# get river floodplain volume
inv_rename = {v: k for k, v in self._MAPS.items() if v in self.staticmaps}
da_fldpln = workflows.river_floodplain_volume(
ds=ds_hydro,
ds_model=self.staticmaps.rename(inv_rename),
river_upa=new_river_upa,
flood_depths=flood_depths,
logger=self.logger,
)
# check if the layer already exists, since overwriting with different flood_depth values is not working properly if this is the case
if "floodplain_volume" in self.staticmaps:
self.logger.warning(
"Layer `floodplain_volume` already in staticmaps, removing layer and `flood_depth` dimension to ensure correctly setting new flood_depth dimensions"
)
self._staticmaps = self._staticmaps.drop_dims("flood_depth")
self.set_staticmaps(da_fldpln, "floodplain_volume")
elif floodplain_type == "2d":
floodplain_1d = False
land_routing = "local-inertial"
if not elevtn_map in self.staticmaps:
raise ValueError(f'"{elevtn_map}" not found in staticmaps')
postfix = {"wflow_dem": "_avg", "dem_subgrid": "_subgrid"}.get(
elevtn_map, ""
)
name = f"hydrodem{postfix}"
self.logger.info(f"Preparing {name} map for land routing.")
name = f"hydrodem{postfix}_D{connectivity}"
self.logger.info(f"Preparing {name} map for land routing.")
ds_out = flw.dem_adjust(
da_flwdir=self.staticmaps[self._MAPS["flwdir"]],
da_elevtn=self.staticmaps[elevtn_map],
da_rivmsk=self.staticmaps[self._MAPS["rivmsk"]],
flwdir=self.flwdir,
connectivity=connectivity,
river_d8=True,
logger=self.logger,
).rename(name)
self.set_staticmaps(ds_out)
# Update config
self.logger.debug(f'Update wflow config model.floodplain_1d="{floodplain_1d}"')
self.set_config("model.floodplain_1d", floodplain_1d)
self.logger.debug(f'Update wflow config model.land_routing="{land_routing}"')
self.set_config("model.land_routing", land_routing)
if floodplain_type == "1d":
# include new input data
self.set_config(
"input.lateral.river.floodplain.volume", "floodplain_volume"
)
# Add states
self.set_config("state.lateral.river.floodplain.q", "q_floodplain")
self.set_config("state.lateral.river.floodplain.h", "h_floodplain")
self.set_config("state.lateral.land.q", "q_land")
# Remove local-inertial land states
if self.get_config("state.lateral.land.qx") is not None:
self.config["state"]["lateral"]["land"].pop("qx", None)
if self.get_config("state.lateral.land.qy") is not None:
self.config["state"]["lateral"]["land"].pop("qy", None)
if self.get_config("output.lateral.land.qx") is not None:
self.config["output"]["lateral"]["land"].pop("qx", None)
if self.get_config("output.lateral.land.qy") is not None:
self.config["output"]["lateral"]["land"].pop("qy", None)
else:
# include new input data
self.set_config("input.lateral.river.bankfull_elevation", name)
self.set_config("input.lateral.land.elevation", name)
# Add local-inertial land routing states
self.set_config("state.lateral.land.qx", "qx_land")
self.set_config("state.lateral.land.qy", "qy_land")
# Remove kinematic-wave and 1d floodplain states
if self.get_config("state.lateral.land.q") is not None:
self.config["state"]["lateral"]["land"].pop("q", None)
if self.get_config("state.lateral.river.floodplain.q") is not None:
self.config["state"]["lateral"]["river"]["floodplain"].pop("q", None)
if self.get_config("state.lateral.river.floodplain.h") is not None:
self.config["state"]["lateral"]["river"]["floodplain"].pop("h", None)
if self.get_config("output.lateral.land.q") is not None:
self.config["output"]["lateral"]["land"].pop("q", None)
def setup_riverwidth(
self,
predictor="discharge",
fill=False,
fit=False,
min_wth=1.0,
precip_fn="chelsa",
climate_fn="koppen_geiger",
**kwargs,
):
"""
This component sets the river width parameter based on a power-lay relationship with a predictor.
By default the riverwidth is estimated based on discharge as ``predictor`` and used to
set the riverwidth globally based on pre-defined power-law parameters per climate class.
With ``fit`` set to True, the power-law relationsship paramters are set on-the-fly.
With ``fill`` set to True, the estimated river widths are only used fill gaps in the
observed data. Alternative ``predictor`` values are precip (accumulated precipitation)
and uparea (upstream area). For these predictors values ``fit`` default to True.
By default the predictor is based on discharge which is estimated through multiple linear
regression with precipitation and upstream area per climate zone.
* **wflow_riverwidth** map: river width [m]
Parameters
----------
predictor : {"discharge", "precip", "uparea"}
Predictor used in the power-law equation: width = a * predictor ^ b.
Discharge is based on multiple linear regression per climate zone.
Precip is based on the 10x the daily average accumulated precipitation [m3/s].
Uparea is based on the upstream area grid [km2].
Other variables, e.g. bankfull discharge, can also be provided if present in the staticmaps
fill : bool, optional
If True (default), use estimate to fill gaps, outliers and lake/res areas in observed width data (if present);
if False, set all riverwidths based on predictor (automatic choice if no observations found)
fit : bool, optional kwarg
If True, the power-law parameters are fitted on the fly
By default True for all but "discharge" predictor.
A-priori derived parameters will be overwritten if True.
a, b : float, optional kwarg
Manual power-law parameters
min_wth : float
minimum river width
precip_fn : {'chelsa'}
Source of long term precipitation grid if the predictor is set to 'discharge' or 'precip'.
climate_fn: {'koppen_geiger'}
Source of long-term climate grid if the predictor is set to 'discharge'.
"""
self.logger.warning(
'The "setup_riverwidth" method has been deprecated and will soon be removed. '
'You can now use the "setup_river" method for all river parameters.'
)
if not self._MAPS["rivmsk"] in self.staticmaps:
raise ValueError(
'The "setup_riverwidth" method requires to run setup_river method first.'
)
# derive river width
data = {}
if predictor in ["discharge", "precip"]:
da_precip = self.data_catalog.get_rasterdataset(
precip_fn, geom=self.region, buffer=2
)
da_precip.name = precip_fn
data["da_precip"] = da_precip
if predictor == "discharge":
da_climate = self.data_catalog.get_rasterdataset(
climate_fn, geom=self.region, buffer=2
)
da_climate.name = climate_fn
data["da_climate"] = da_climate
inv_rename = {v: k for k, v in self._MAPS.items() if v in self.staticmaps}
da_rivwth = workflows.river_width(
ds_like=self.staticmaps.rename(inv_rename),
flwdir=self.flwdir,
data=data,
fill=fill,
fill_outliers=kwargs.pop("fill_outliers", fill),
min_wth=min_wth,
mask_names=["lakeareas", "resareas", "glacareas"],
predictor=predictor,
a=kwargs.get("a", None),
b=kwargs.get("b", None),
logger=self.logger,
fit=fit,
**kwargs,
)
self.set_staticmaps(da_rivwth, name=self._MAPS["rivwth"])
def setup_lulcmaps(
self,
lulc_fn="globcover",
lulc_mapping_fn=None,
lulc_vars=[
"landuse",
"Kext",
"N",
"PathFrac",
"RootingDepth",
"Sl",
"Swood",
"WaterFrac",
],
):
"""
This component derives several wflow maps are derived based on landuse-
landcover (LULC) data.
Currently, ``lulc_fn`` can be set to the "vito", "globcover", "esa_worldcover"
or "corine", of which lookup tables are constructed to convert lulc classses to
model parameters based on literature. The data is remapped at its original
resolution and then resampled to the model resolution using the average
value, unless noted differently.
Adds model layers:
* **landuse** map: Landuse class [-]
* **Kext** map: Extinction coefficient in the canopy gap fraction equation [-]
* **Sl** map: Specific leaf storage [mm]
* **Swood** map: Fraction of wood in the vegetation/plant [-]
* **RootingDepth** map: Length of vegetation roots [mm]
* **PathFrac** map: The fraction of compacted or urban area per grid cell [-]
* **WaterFrac** map: The fraction of open water per grid cell [-]
* **N** map: Manning Roughness [-]
Parameters
----------
lulc_fn : {"globcover", "vito", "corine"}
Name of data source in data_sources.yml file.
lulc_mapping_fn : str
Path to a mapping csv file from landuse in source name to parameter values in lulc_vars.
lulc_vars : list
List of landuse parameters to keep.\
By default ["landuse","Kext","N","PathFrac","RootingDepth","Sl","Swood","WaterFrac"]
"""
self.logger.info("Preparing LULC parameter maps.")
if lulc_mapping_fn is None:
fn_map = f"{lulc_fn}_mapping_default"
else:
fn_map = lulc_mapping_fn
if not isfile(fn_map) and fn_map not in self.data_catalog:
raise ValueError(f"LULC mapping file not found: {fn_map}")
# read landuse map to DataArray
da = self.data_catalog.get_rasterdataset(
lulc_fn, geom=self.region, buffer=2, variables=["landuse"]
)
df_map = self.data_catalog.get_dataframe(
fn_map,
driver_kwargs={"index_col": 0}, # only used if fn_map is a file path
)
# process landuse
ds_lulc_maps = workflows.landuse(
da=da,
ds_like=self.staticmaps,
df=df_map,
params=lulc_vars,
logger=self.logger,
)
rmdict = {k: v for k, v in self._MAPS.items() if k in ds_lulc_maps.data_vars}
self.set_staticmaps(ds_lulc_maps.rename(rmdict))
def setup_laimaps(self, lai_fn="modis_lai"):
"""
This component sets leaf area index (LAI) climatology maps per month.
The values are resampled to the model resolution using the average value.
The only ``lai_fn`` currently supported is "modis_lai" based on MODIS data.
Adds model layers:
* **LAI** map: Leaf Area Index climatology [-]
Resampled from source data using average. Assuming that missing values
correspond to bare soil, these are set to zero before resampling.
Parameters
----------
lai_fn : {'modis_lai'}
Name of data source for LAI parameters, see data/data_sources.yml.
* Required variables: ['LAI']
"""
if lai_fn not in self.data_catalog:
self.logger.warning(f"Invalid source '{lai_fn}', skipping setup_laimaps.")
return
# retrieve data for region
self.logger.info(f"Preparing LAI maps.")
da = self.data_catalog.get_rasterdataset(lai_fn, geom=self.region, buffer=2)
da_lai = workflows.lai(
da=da,
ds_like=self.staticmaps,
logger=self.logger,
)
# Rename the first dimension to time
rmdict = {da_lai.dims[0]: "time"}
self.set_staticmaps(da_lai.rename(rmdict), name="LAI")
def setup_config_output_timeseries(
self,
mapname: str,
toml_output: Optional[str] = "csv",
header: Optional[List[str]] = ["Q"],
param: Optional[List[str]] = ["lateral.river.q_av"],
reducer: Optional[List[str]] = None,
):
"""This components sets the default gauge map based on basin outlets.
Adds model layers:
* **csv.column** config: csv timeseries to save based on mapname locations
* **netcdf.variable** config: netcdf timeseries to save based on mapname locations
Parameters
----------
mapname : str
Name of the gauge map (in staticmaps.nc) to use for scalar output.
toml_output : str, optional
One of ['csv', 'netcdf', None] to update [csv] or [netcdf] section of wflow toml file or do nothing. By default, 'csv'.
header : list, optional
Save specific model parameters in csv section. This option defines the header of the csv file./
By default saves Q (for lateral.river.q_av).
param: list, optional
Save specific model parameters in csv section. This option defines the wflow variable corresponding to the/
names in gauge_toml_header. By default saves lateral.river.q_av (for Q).
reducer: list, optional
If map is an area rather than a point location, provides the reducer for the parameters to save. By default None.
"""
# # Add new outputcsv section in the config
if toml_output == "csv" or toml_output == "netcdf":
self.logger.info(f"Adding {param} to {toml_output} section of toml.")
# Add map to the input section of config
basename = (
mapname
if not mapname.startswith("wflow")
else mapname.replace("wflow_", "")
)
self.set_config(f"input.{basename}", mapname)
# Settings and add csv or netcdf sections if not already in config
# csv
if toml_output == "csv":
header_name = "header"
var_name = "column"
if self.get_config("csv") is None:
self.set_config("csv.path", "output.csv")
# netcdf
if toml_output == "netcdf":
header_name = "name"
var_name = "variable"
if self.get_config("netcdf") is None:
self.set_config("netcdf.path", "output_scalar.nc")
# initialise column / varibale section
if self.get_config(f"{toml_output}.{var_name}") is None:
self.set_config(f"{toml_output}.{var_name}", [])
# Add new output column/variable to config
for o in range(len(param)):
gauge_toml_dict = {
header_name: header[o],
"map": basename,
"parameter": param[o],
}
if reducer is not None:
gauge_toml_dict["reducer"]: reducer[o]
# If the gauge column/variable already exists skip writting twice
if gauge_toml_dict not in self.config[toml_output][var_name]:
self.config[toml_output][var_name].append(gauge_toml_dict)
else:
self.logger.info(
f"toml_output set to {toml_output}, skipping adding gauge specific outputs to the toml."
)
def setup_outlets(
self,
river_only=True,
toml_output="csv",
gauge_toml_header=["Q"],
gauge_toml_param=["lateral.river.q_av"],
):
"""This components sets the default gauge map based on basin outlets.
Adds model layers:
* **wflow_gauges** map: gauge IDs map from catchment outlets [-]
* **gauges** geom: polygon of catchment outlets
Parameters
----------
river_only : bool, optional
Only derive outlet locations if they are located on a river instead of locations for all catchments, by default True.
toml_output : str, optional
One of ['csv', 'netcdf', None] to update [csv] or [netcdf] section of wflow toml file or do nothing. By default, 'csv'.
gauge_toml_header : list, optional
Save specific model parameters in csv section. This option defines the header of the csv file./
By default saves Q (for lateral.river.q_av).
gauge_toml_param: list, optional
Save specific model parameters in csv section. This option defines the wflow variable corresponding to the/
names in gauge_toml_header. By default saves lateral.river.q_av (for Q).
"""
# read existing staticgeoms; important to get the right basin when updating
# fix in set_staticgeoms / set_geoms method
self.staticgeoms
self.logger.info(f"Gauges locations set based on river outlets.")
idxs_out = self.flwdir.idxs_pit
# Only keep river outlets for gauges
if river_only:
idxs_out = idxs_out[
(self.staticmaps[self._MAPS["rivmsk"]] > 0).values.flat[idxs_out]
]