From 80fdcc46f6d20f5359d3aa3801d8bc03f6bebaaa Mon Sep 17 00:00:00 2001 From: hboisgon Date: Tue, 28 Mar 2023 21:45:14 +0800 Subject: [PATCH 01/17] split setup gauges and outlets + fixes #76 --- examples/wflow_build.ini | 3 + examples/wflow_sediment_build.ini | 3 + hydromt_wflow/data/wflow/wflow_sbm.toml | 24 +- .../data/wflow_sediment/wflow_sediment.toml | 8 +- hydromt_wflow/wflow.py | 320 ++++++++++-------- hydromt_wflow/wflow_sediment.py | 48 ++- tests/data/wflow_piave_build_subbasin.ini | 3 + .../wflow_sediment_piave_build_subbasin.ini | 3 + 8 files changed, 230 insertions(+), 182 deletions(-) diff --git a/examples/wflow_build.ini b/examples/wflow_build.ini index 427df2a7..c230d5b9 100644 --- a/examples/wflow_build.ini +++ b/examples/wflow_build.ini @@ -48,6 +48,9 @@ lai_fn = modis_lai # source for LAI: {modis_lai} soil_fn = soilgrids # source for soilmaps: {soilgrids} ptf_ksatver = brakensiek # pedotransfer function to calculate hydraulic conductivity: {brakensiek, cosby} +[setup_outlets] +river_only = True + [setup_gauges] gauges_fn = grdc # if not None add gaugemap. Either a path or known gauges_fn: {grdc} snap_to_river = True # if True snaps gauges from source to river diff --git a/examples/wflow_sediment_build.ini b/examples/wflow_sediment_build.ini index 8a1d5700..f47ef91f 100644 --- a/examples/wflow_sediment_build.ini +++ b/examples/wflow_sediment_build.ini @@ -47,6 +47,9 @@ canopy_fn = simard # source for vegetation canopy height: {sima soil_fn = soilgrids # source for soilmaps: {soilgrids} usleK_method = renard # method to compute the USLE K factor: {renard, epic} +[setup_outlets] +river_only = True + [setup_gauges] gauges_fn = grdc # If not None add gaugemap. Either a path or known gauges_fn: {grdc} snap_to_river = True # If True snaps gauges from source to river diff --git a/hydromt_wflow/data/wflow/wflow_sbm.toml b/hydromt_wflow/data/wflow/wflow_sbm.toml index 25e55955..e3155680 100644 --- a/hydromt_wflow/data/wflow/wflow_sbm.toml +++ b/hydromt_wflow/data/wflow/wflow_sbm.toml @@ -46,7 +46,6 @@ path_forcing = "inmaps.nc" path_static = "staticmaps.nc" # these are not directly part of the model -gauges = "wflow_gauges" ldd = "wflow_ldd" river_location = "wflow_river" subcatchment = "wflow_subcatch" @@ -127,29 +126,8 @@ thicknesslayers = [100, 300, 800] [output] path = "output.nc" -#[output.vertical] -#satwaterdepth = "satwaterdepth" -#snow = "snow" -#tsoil = "tsoil" -#ustorelayerdepth = "ustorelayerdepth" -#snowwater = "snowwater" -#canopystorage = "canopystorage" - [output.lateral.river] q_av = "q_river" -#h = "h_river" - -#[output.lateral.subsurface] -#ssf = "ssf" - -[output.lateral.land] -q = "q_land" -h = "h_land" [csv] -path = "output.csv" - -[[csv.column]] -header = "Q" -map = "gauges" -parameter = "lateral.river.q_av" +path = "output.csv" \ No newline at end of file diff --git a/hydromt_wflow/data/wflow_sediment/wflow_sediment.toml b/hydromt_wflow/data/wflow_sediment/wflow_sediment.toml index 0d230e58..9aed806a 100644 --- a/hydromt_wflow/data/wflow_sediment/wflow_sediment.toml +++ b/hydromt_wflow/data/wflow_sediment/wflow_sediment.toml @@ -45,7 +45,6 @@ path_forcing = "inmaps.nc" path_static = "staticmaps.nc" # these are not directly part of the model -gauges = "wflow_gauges" ldd = "wflow_ldd" river_location = "wflow_river" subcatchment = "wflow_subcatch" @@ -136,9 +135,4 @@ soilloss = "soilloss" SSconc = "SSconc" [csv] -path = "output.csv" - -[[csv.column]] -header = "TSS" -map = "gauges" -parameter = "lateral.river.SSconc" \ No newline at end of file +path = "output.csv" \ No newline at end of file diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 915d1295..dfa06758 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -467,20 +467,26 @@ def setup_hydrodem( self.set_config("model.river_routing", river_routing) self.logger.debug(f'Update wflow config model.land_routing="{land_routing}"') self.set_config("model.land_routing", land_routing) + remove_landq = False + if "output" in self.config["output"]: + if "land" in self.config["output"]["lateral"]: + remove_landq = True if river_routing == "local-inertial": self.set_config("input.lateral.river.bankfull_depth", self._MAPS["rivdph"]) self.set_config("input.lateral.river.bankfull_elevation", name) if land_routing == "local-inertial": self.set_config("input.lateral.land.elevation", name) self.config["state"]["lateral"]["land"].pop("q", None) - self.config["output"]["lateral"]["land"].pop("q", None) + if remove_landq: + self.config["output"]["lateral"]["land"].pop("q", None) self.set_config("state.lateral.land.qx", "qx_land") self.set_config("state.lateral.land.qy", "qy_land") else: self.config["state"]["lateral"]["land"].pop("qx", None) self.config["state"]["lateral"]["land"].pop("qy", None) - self.config["output"]["lateral"]["land"].pop("qx", None) - self.config["output"]["lateral"]["land"].pop("qy", None) + if remove_landq: + self.config["output"]["lateral"]["land"].pop("qx", None) + self.config["output"]["lateral"]["land"].pop("qy", None) self.set_config("state.lateral.land.q", "q_land") def setup_riverwidth( @@ -678,6 +684,70 @@ def setup_laimaps(self, lai_fn="modis_lai"): rmdict = {da_lai.dims[0]: "time"} self.set_staticmaps(da_lai.rename(rmdict), name="LAI") + def setup_outlets( + self, + river_only=True, + update_toml=True, + 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. + update_toml : boolean, optional + Update [outputcsv] section of wflow toml file. + 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 + 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] + ] + da_out, idxs_out, ids_out = flw.gauge_map( + self.staticmaps, idxs=self.flwdir.idxs_out + ) + self.set_staticmaps(da_out, name=self._MAPS["gauges"]) + points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs_out)) + gdf = gpd.GeoDataFrame( + index=ids_out.astype(np.int32), geometry=points, crs=self.crs + ) + gdf["fid"] = ids_out.astype(np.int32) + self.set_staticgeoms(gdf, name="gauges") + self.logger.info(f"Gauges map based on catchment river outlets added.") + + # # Add new outputcsv section in the config + if update_toml: + self.logger.info(f"Adding {gauge_toml_param} to csv output.") + self.set_config("input.gauges", "wflow_gauges") + if self.get_config("csv") is not None: + for o in range(len(gauge_toml_param)): + gauge_toml_dict = { + "header": gauge_toml_header[o], + "map": "gauges", + "parameter": gauge_toml_param[o], + } + # If the gauge outcsv column already exists skip writting twice + if gauge_toml_dict not in self.config["csv"]["column"]: + self.config["csv"]["column"].append(gauge_toml_dict) + def setup_gauges( self, gauges_fn="grdc", @@ -686,15 +756,13 @@ def setup_gauges( snap_to_river=True, mask=None, derive_subcatch=False, - derive_outlet=True, basename=None, update_toml=True, - gauge_toml_header=None, - gauge_toml_param=None, + gauge_toml_header=["Q", "P"], + gauge_toml_param=["lateral.river.q_av", "vertical.precipitation"], **kwargs, ): - """This components sets the default gauge map based on basin outlets and additional - gauge maps based on ``gauges_fn`` data. + """This components sets a gauge map based on ``gauges_fn`` data. Supported gauge datasets include "grdc" or "" for user supplied csv or geometry files with gauge locations. @@ -706,10 +774,8 @@ def setup_gauges( Adds model layers: - * **wflow_gauges** map: gauge IDs map from catchment outlets [-] * **wflow_gauges_source** map: gauge IDs map from source [-] (if gauges_fn) * **wflow_subcatch_source** map: subcatchment based on gauge locations [-] (if derive_subcatch) - * **gauges** geom: polygon of catchment outlets * **gauges_source** geom: polygon of gauges from source * **subcatch_source** geom: polygon of subcatchment based on gauge locations [-] (if derive_subcatch) @@ -727,8 +793,6 @@ def setup_gauges( If provided snaps to the mask, else snaps to the river (default). derive_subcatch : bool, optional Derive subcatch map for gauges, by default False - derive_outlet : bool, optional - Derive gaugemap based on catchment outlets, by default True basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. update_toml : boolean, optional @@ -743,139 +807,113 @@ def setup_gauges( # read existing staticgeoms; important to get the right basin when updating self.staticgeoms - if derive_outlet: - self.logger.info(f"Gauges locations set based on river outlets.") - da, idxs, ids = flw.gauge_map(self.staticmaps, idxs=self.flwdir.idxs_pit) - # Only keep river outlets for gauges - da = da.where(self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata) - ids_da = np.unique(da.values[da.values > 0]) - idxs_da = idxs[np.isin(ids, ids_da)] - self.set_staticmaps(da, name=self._MAPS["gauges"]) - points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs_da)) - gdf = gpd.GeoDataFrame( - index=ids_da.astype(np.int32), geometry=points, crs=self.crs + # TODO check snapping locations based on upstream area attribute of the gauge data + if gauges_fn is not None: + kwargs = {} + if isfile(gauges_fn): + # try to get epsg number directly, important when writting back data_catalog + if hasattr(self.crs, "to_epsg"): + code = self.crs.to_epsg() + else: + code = self.crs + kwargs.update(crs=code) + gdf = self.data_catalog.get_geodataframe( + gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs + ) + gdf = gdf.to_crs(self.crs) + elif source_gdf is not None and basename is None: + raise ValueError( + "Basename is required when setting gauges based on source_gdf" + ) + elif source_gdf is not None: + self.logger.info(f"Gauges locations read from source_gdf") + gdf = source_gdf.to_crs(self.crs) + else: + self.logger.warning( + "Either gauges_fn or source_gdf should be provided, skipping setup_gauges." ) - gdf["fid"] = ids_da.astype(np.int32) - self.set_staticgeoms(gdf, name="gauges") - self.logger.info(f"Gauges map based on catchment river outlets added.") - - if gauges_fn is not None or source_gdf is not None: - # append location from geometry - # TODO check snapping locations based on upstream area attribute of the gauge data - if gauges_fn is not None: - kwargs = {} - if isfile(gauges_fn): - # try to get epsg number directly, important when writting back data_catalog - if hasattr(self.crs, "to_epsg"): - code = self.crs.to_epsg() - else: - code = self.crs - kwargs.update(crs=code) - gdf = self.data_catalog.get_geodataframe( - gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs - ) - gdf = gdf.to_crs(self.crs) - elif source_gdf is not None and basename is None: - raise ValueError( - "Basename is required when setting gauges based on source_gdf" - ) - elif source_gdf is not None: - self.logger.info(f"Gauges locations read from source_gdf") - gdf = source_gdf.to_crs(self.crs) - if gdf.index.size == 0: - self.logger.warning( - f"No {gauges_fn} gauge locations found within domain" + if gdf.index.size == 0: + self.logger.warning(f"No {gauges_fn} gauge locations found within domain") + else: + if basename is None: + basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") + self.logger.info( + f"{gdf.index.size} {basename} gauge locations found within domain" + ) + # Set index to index_col + if index_col is not None and index_col in gdf: + gdf = gdf.set_index(index_col) + xs, ys = np.vectorize(lambda p: (p.xy[0][0], p.xy[1][0]))(gdf["geometry"]) + idxs = self.staticmaps.raster.xy_to_idx(xs, ys) + ids = gdf.index.values + + if snap_to_river and mask is None: + mask = self.staticmaps[self._MAPS["rivmsk"]].values + da, idxs, ids = flw.gauge_map( + self.staticmaps, + idxs=idxs, + ids=ids, + stream=mask, + flwdir=self.flwdir, + logger=self.logger, + ) + # Filter gauges that could not be snapped to rivers + if snap_to_river: + ids_old = ids.copy() + da = da.where( + self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata ) + ids_new = np.unique(da.values[da.values > 0]) + idxs = idxs[np.isin(ids_old, ids_new)] + ids = da.values.flat[idxs] + # Add to staticmaps + mapname = f'{str(self._MAPS["gauges"])}_{basename}' + self.set_staticmaps(da, name=mapname) + + # geoms + points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs)) + # if csv contains additional columns, these are also written in the staticgeoms + gdf_snapped = gpd.GeoDataFrame( + index=ids.astype(np.int32), geometry=points, crs=self.crs + ) + # Set the index name of gdf snapped based on original gdf + if gdf.index.name is not None: + gdf_snapped.index.name = gdf.index.name else: - if basename is None: - basename = ( - os.path.basename(gauges_fn).split(".")[0].replace("_", "-") - ) - self.logger.info( - f"{gdf.index.size} {basename} gauge locations found within domain" - ) - # Set index to index_col - if index_col is not None and index_col in gdf: - gdf = gdf.set_index(index_col) - xs, ys = np.vectorize(lambda p: (p.xy[0][0], p.xy[1][0]))( - gdf["geometry"] - ) - idxs = self.staticmaps.raster.xy_to_idx(xs, ys) - ids = gdf.index.values - - if snap_to_river and mask is None: - mask = self.staticmaps[self._MAPS["rivmsk"]].values - da, idxs, ids = flw.gauge_map( - self.staticmaps, - idxs=idxs, - ids=ids, - stream=mask, - flwdir=self.flwdir, - logger=self.logger, - ) - # Filter gauges that could not be snapped to rivers - if snap_to_river: - ids_old = ids.copy() - da = da.where( - self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata - ) - ids_new = np.unique(da.values[da.values > 0]) - idxs = idxs[np.isin(ids_old, ids_new)] - ids = da.values.flat[idxs] - # Add to staticmaps - mapname = f'{str(self._MAPS["gauges"])}_{basename}' - self.set_staticmaps(da, name=mapname) - - # geoms - points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs)) - # if csv contains additional columns, these are also written in the staticgeoms - gdf_snapped = gpd.GeoDataFrame( - index=ids.astype(np.int32), geometry=points, crs=self.crs - ) - # Set the index name of gdf snapped based on original gdf - if gdf.index.name is not None: - gdf_snapped.index.name = gdf.index.name - else: - gdf_snapped.index.name = "fid" - gdf.index.name = "fid" - # Add gdf attributes to gdf_snapped (filter on snapped index before merging) - df_attrs = pd.DataFrame(gdf.drop(columns="geometry")) - df_attrs = df_attrs[np.isin(df_attrs.index, gdf_snapped.index)] - gdf_snapped = gdf_snapped.merge( - df_attrs, how="inner", on=gdf.index.name - ) - # Add gdf_snapped to staticgeoms - self.set_staticgeoms(gdf_snapped, name=mapname.replace("wflow_", "")) - - # # Add new outputcsv section in the config - if gauge_toml_param is None and update_toml: - gauge_toml_header = ["Q", "P"] - gauge_toml_param = ["lateral.river.q_av", "vertical.precipitation"] - - if update_toml: - self.set_config(f"input.gauges_{basename}", f"{mapname}") - if self.get_config("csv") is not None: - for o in range(len(gauge_toml_param)): - gauge_toml_dict = { - "header": gauge_toml_header[o], - "map": f"gauges_{basename}", - "parameter": gauge_toml_param[o], - } - # If the gauge outcsv column already exists skip writting twice - if gauge_toml_dict not in self.config["csv"]["column"]: - self.config["csv"]["column"].append(gauge_toml_dict) - self.logger.info(f"Gauges map from {basename} added.") - - # add subcatch - if derive_subcatch: - da_basins = flw.basin_map( - self.staticmaps, self.flwdir, idxs=idxs, ids=ids - )[0] - mapname = self._MAPS["basins"] + "_" + basename - self.set_staticmaps(da_basins, name=mapname) - gdf_basins = self.staticmaps[mapname].raster.vectorize() - self.set_staticgeoms(gdf_basins, name=mapname.replace("wflow_", "")) + gdf_snapped.index.name = "fid" + gdf.index.name = "fid" + # Add gdf attributes to gdf_snapped (filter on snapped index before merging) + df_attrs = pd.DataFrame(gdf.drop(columns="geometry")) + df_attrs = df_attrs[np.isin(df_attrs.index, gdf_snapped.index)] + gdf_snapped = gdf_snapped.merge(df_attrs, how="inner", on=gdf.index.name) + # Add gdf_snapped to staticgeoms + self.set_staticgeoms(gdf_snapped, name=mapname.replace("wflow_", "")) + + # # Add new outputcsv section in the config + if update_toml: + self.set_config(f"input.gauges_{basename}", f"{mapname}") + if self.get_config("csv") is not None: + for o in range(len(gauge_toml_param)): + gauge_toml_dict = { + "header": gauge_toml_header[o], + "map": f"gauges_{basename}", + "parameter": gauge_toml_param[o], + } + # If the gauge outcsv column already exists skip writting twice + if gauge_toml_dict not in self.config["csv"]["column"]: + self.config["csv"]["column"].append(gauge_toml_dict) + self.logger.info(f"Gauges map from {basename} added.") + + # add subcatch + if derive_subcatch: + da_basins = flw.basin_map( + self.staticmaps, self.flwdir, idxs=idxs, ids=ids + )[0] + mapname = self._MAPS["basins"] + "_" + basename + self.set_staticmaps(da_basins, name=mapname) + gdf_basins = self.staticmaps[mapname].raster.vectorize() + self.set_staticgeoms(gdf_basins, name=mapname.replace("wflow_", "")) def setup_areamap( self, diff --git a/hydromt_wflow/wflow_sediment.py b/hydromt_wflow/wflow_sediment.py index 1522b06f..cb1d2d03 100644 --- a/hydromt_wflow/wflow_sediment.py +++ b/hydromt_wflow/wflow_sediment.py @@ -170,6 +170,40 @@ def setup_reservoirs( if self.get_config("input.lateral.river.reservoir") is not None: del self.config["input"]["lateral"]["river"]["reservoir"] + def setup_outlets( + self, + river_only=True, + update_toml=True, + gauge_toml_header=["TSS"], + gauge_toml_param=["lateral.river.SSconc"], + ): + """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. + update_toml : boolean, optional + Update [outputcsv] section of wflow toml file. + gauge_toml_header : list, optional + Save specific model parameters in csv section. This option defines the header of the csv file./ + By default saves TSS (for lateral.river.SSconc). + 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.SSconc (for TSS). + """ + super().setup_outlets( + river_only=river_only, + update_toml=update_toml, + gauge_toml_header=gauge_toml_header, + gauge_toml_param=gauge_toml_param, + ) + def setup_gauges( self, gauges_fn=None, @@ -177,15 +211,13 @@ def setup_gauges( snap_to_river=True, mask=None, derive_subcatch=False, - derive_outlet=True, basename=None, update_toml=True, - gauge_toml_header=None, - gauge_toml_param=None, + gauge_toml_header=["Q", "TSS"], + gauge_toml_param=["lateral.river.q_riv", "lateral.river.SSconc"], **kwargs, ): - """This components sets the default gauge map based on basin outlets and additional - gauge maps based on ``gauges_fn`` data. + """This components sets a gauge map based on ``gauges_fn`` data. Supported gauge datasets include "grdc" or "" for user supplied csv or geometry files with gauge locations. @@ -197,10 +229,8 @@ def setup_gauges( Adds model layers: - * **wflow_gauges** map: gauge IDs map from catchment outlets [-] * **wflow_gauges_source** map: gauge IDs map from source [-] (if gauges_fn) * **wflow_subcatch_source** map: subcatchment based on gauge locations [-] (if derive_subcatch) - * **gauges** geom: polygon of catchment outlets * **gauges_source** geom: polygon of gauges from source * **subcatch_source** geom: polygon of subcatchment based on gauge locations [-] (if derive_subcatch) @@ -230,16 +260,12 @@ def setup_gauges( names in gauge_toml_header. By default saves lateral.river.q_riv (for Q) and lateral.river.SSconc (for TSS). """ # # Add new outputcsv section in the config - if gauge_toml_param is None and update_toml: - gauge_toml_header = ["Q", "TSS"] - gauge_toml_param = ["lateral.river.q_riv", "lateral.river.SSconc"] super().setup_gauges( gauges_fn=gauges_fn, source_gdf=source_gdf, snap_to_river=snap_to_river, mask=mask, derive_subcatch=derive_subcatch, - derive_outlet=derive_outlet, basename=basename, update_toml=update_toml, gauge_toml_header=gauge_toml_header, diff --git a/tests/data/wflow_piave_build_subbasin.ini b/tests/data/wflow_piave_build_subbasin.ini index aea5e46d..9c105fbf 100644 --- a/tests/data/wflow_piave_build_subbasin.ini +++ b/tests/data/wflow_piave_build_subbasin.ini @@ -43,6 +43,9 @@ lai_fn = modis_lai soil_fn = soilgrids ptf_ksatver = brakensiek +[setup_outlets] +river_only = True + [setup_gauges] gauges_fn = grdc snap_to_river = True diff --git a/tests/data/wflow_sediment_piave_build_subbasin.ini b/tests/data/wflow_sediment_piave_build_subbasin.ini index 31e63c68..080cc1f2 100644 --- a/tests/data/wflow_sediment_piave_build_subbasin.ini +++ b/tests/data/wflow_sediment_piave_build_subbasin.ini @@ -43,6 +43,9 @@ canopy_fn = simard soil_fn = soilgrids usleK_method = renard +[setup_outlets] +river_only = True + [setup_gauges] gauges_fn = grdc snap_to_river = True From ebbeb5e43d277863a8ab179cbc0b043420870ac0 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Tue, 28 Mar 2023 22:22:35 +0800 Subject: [PATCH 02/17] support saving netcdf scalar output for gauges --- hydromt_wflow/wflow.py | 117 ++++++++++++++++++++++---------- hydromt_wflow/wflow_sediment.py | 46 ++++++------- 2 files changed, 106 insertions(+), 57 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index dfa06758..85892786 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -684,10 +684,75 @@ def setup_laimaps(self, lai_fn="modis_lai"): rmdict = {da_lai.dims[0]: "time"} self.set_staticmaps(da_lai.rename(rmdict), name="LAI") + def _setup_config_timeseries( + self, + mapname: str, + 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 + ---------- + mapname : str + Name of the gauge map to use for scalar output (without wflow_). + 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). + """ + + # # Add new outputcsv section in the config + if toml_output == "csv" or toml_output == "netcdf": + self.logger.info( + f"Adding {gauge_toml_param} to {toml_output} section of toml." + ) + self.set_config(f"input.{mapname}", f"wflow_{mapname}") + # csv + if toml_output == "csv": + if self.get_config("csv") is None: + self.set_config("csv.path", "output.csv") + for o in range(len(gauge_toml_param)): + gauge_toml_dict = { + "header": gauge_toml_header[o], + "map": mapname, + "parameter": gauge_toml_param[o], + } + # If the gauge outcsv column already exists skip writting twice + if gauge_toml_dict not in self.config["csv"]["column"]: + self.config["csv"]["column"].append(gauge_toml_dict) + # netcdf + if toml_output == "netcdf": + if self.get_config("netcdf") is None: + self.set_config("netcdf.path", "output_scalar.nc") + for o in range(len(gauge_toml_param)): + gauge_toml_dict = { + "name": gauge_toml_header[o], + "map": mapname, + "parameter": gauge_toml_param[o], + } + # If the gauge outcsv column already exists skip writting twice + if gauge_toml_dict not in self.config["netcdf"]["variable"]: + self.config["netcdf"]["variable"].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, - update_toml=True, + toml_output="csv", gauge_toml_header=["Q"], gauge_toml_param=["lateral.river.q_av"], ): @@ -702,8 +767,8 @@ def setup_outlets( ---------- river_only : bool, optional Only derive outlet locations if they are located on a river instead of locations for all catchments, by default True. - update_toml : boolean, optional - Update [outputcsv] section of wflow toml file. + 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). @@ -733,20 +798,12 @@ def setup_outlets( self.set_staticgeoms(gdf, name="gauges") self.logger.info(f"Gauges map based on catchment river outlets added.") - # # Add new outputcsv section in the config - if update_toml: - self.logger.info(f"Adding {gauge_toml_param} to csv output.") - self.set_config("input.gauges", "wflow_gauges") - if self.get_config("csv") is not None: - for o in range(len(gauge_toml_param)): - gauge_toml_dict = { - "header": gauge_toml_header[o], - "map": "gauges", - "parameter": gauge_toml_param[o], - } - # If the gauge outcsv column already exists skip writting twice - if gauge_toml_dict not in self.config["csv"]["column"]: - self.config["csv"]["column"].append(gauge_toml_dict) + self._setup_config_timeseries( + mapname="gauges", + toml_output=toml_output, + gauge_toml_header=gauge_toml_header, + gauge_toml_param=gauge_toml_param, + ) def setup_gauges( self, @@ -757,7 +814,7 @@ def setup_gauges( mask=None, derive_subcatch=False, basename=None, - update_toml=True, + toml_output="csv", gauge_toml_header=["Q", "P"], gauge_toml_param=["lateral.river.q_av", "vertical.precipitation"], **kwargs, @@ -795,8 +852,8 @@ def setup_gauges( Derive subcatch map for gauges, by default False basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. - update_toml : boolean, optional - Update [outputcsv] section of wflow toml file. + 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) and P (for vertical.precipitation). @@ -890,20 +947,12 @@ def setup_gauges( # Add gdf_snapped to staticgeoms self.set_staticgeoms(gdf_snapped, name=mapname.replace("wflow_", "")) - # # Add new outputcsv section in the config - if update_toml: - self.set_config(f"input.gauges_{basename}", f"{mapname}") - if self.get_config("csv") is not None: - for o in range(len(gauge_toml_param)): - gauge_toml_dict = { - "header": gauge_toml_header[o], - "map": f"gauges_{basename}", - "parameter": gauge_toml_param[o], - } - # If the gauge outcsv column already exists skip writting twice - if gauge_toml_dict not in self.config["csv"]["column"]: - self.config["csv"]["column"].append(gauge_toml_dict) - self.logger.info(f"Gauges map from {basename} added.") + self._setup_config_timeseries( + mapname=mapname, + toml_output=toml_output, + gauge_toml_header=gauge_toml_header, + gauge_toml_param=gauge_toml_param, + ) # add subcatch if derive_subcatch: diff --git a/hydromt_wflow/wflow_sediment.py b/hydromt_wflow/wflow_sediment.py index cb1d2d03..8c6cca9f 100644 --- a/hydromt_wflow/wflow_sediment.py +++ b/hydromt_wflow/wflow_sediment.py @@ -173,33 +173,33 @@ def setup_reservoirs( def setup_outlets( self, river_only=True, - update_toml=True, + toml_output="csv", gauge_toml_header=["TSS"], gauge_toml_param=["lateral.river.SSconc"], ): """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. - update_toml : boolean, optional - Update [outputcsv] section of wflow toml file. - gauge_toml_header : list, optional - Save specific model parameters in csv section. This option defines the header of the csv file./ - By default saves TSS (for lateral.river.SSconc). - 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.SSconc (for TSS). + 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 TSS (for lateral.river.SSconc). + 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.SSconc (for TSS). """ super().setup_outlets( river_only=river_only, - update_toml=update_toml, + toml_output=toml_output, gauge_toml_header=gauge_toml_header, gauge_toml_param=gauge_toml_param, ) @@ -212,7 +212,7 @@ def setup_gauges( mask=None, derive_subcatch=False, basename=None, - update_toml=True, + toml_output="csv", gauge_toml_header=["Q", "TSS"], gauge_toml_param=["lateral.river.q_riv", "lateral.river.SSconc"], **kwargs, @@ -250,8 +250,8 @@ def setup_gauges( Derive gaugemap based on catchment outlets, by default True basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. - update_toml : boolean, optional - Update [outputcsv] section of wflow toml file. + 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_riv) and TSS (for lateral.river.SSconc). @@ -267,7 +267,7 @@ def setup_gauges( mask=mask, derive_subcatch=derive_subcatch, basename=basename, - update_toml=update_toml, + toml_output=toml_output, gauge_toml_header=gauge_toml_header, gauge_toml_param=gauge_toml_param, ) From a5b068522b62be6ff04d4ed4a1b2b6ca2e055204 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 30 Mar 2023 00:05:03 +0800 Subject: [PATCH 03/17] new setup_config_output_timeseries + fixes #68 --- hydromt_wflow/wflow.py | 342 +++++++++++++++++++++++++---------------- 1 file changed, 210 insertions(+), 132 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 85892786..4f19e013 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -684,66 +684,75 @@ def setup_laimaps(self, lai_fn="modis_lai"): rmdict = {da_lai.dims[0]: "time"} self.set_staticmaps(da_lai.rename(rmdict), name="LAI") - def _setup_config_timeseries( + def setup_config_output_timeseries( self, mapname: str, - toml_output="csv", - gauge_toml_header=["Q"], - gauge_toml_param=["lateral.river.q_av"], + 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: - * **wflow_gauges** map: gauge IDs map from catchment outlets [-] - * **gauges** geom: polygon of catchment outlets + * **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 to use for scalar output (without wflow_). + 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'. - gauge_toml_header : list, optional + 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 + 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 {gauge_toml_param} to {toml_output} section of toml." + 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.{mapname}", f"wflow_{mapname}") + 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") - for o in range(len(gauge_toml_param)): - gauge_toml_dict = { - "header": gauge_toml_header[o], - "map": mapname, - "parameter": gauge_toml_param[o], - } - # If the gauge outcsv column already exists skip writting twice - if gauge_toml_dict not in self.config["csv"]["column"]: - self.config["csv"]["column"].append(gauge_toml_dict) # 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") - for o in range(len(gauge_toml_param)): - gauge_toml_dict = { - "name": gauge_toml_header[o], - "map": mapname, - "parameter": gauge_toml_param[o], - } - # If the gauge outcsv column already exists skip writting twice - if gauge_toml_dict not in self.config["netcdf"]["variable"]: - self.config["netcdf"]["variable"].append(gauge_toml_dict) + 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." @@ -777,6 +786,7 @@ def setup_outlets( 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.") @@ -798,8 +808,8 @@ def setup_outlets( self.set_staticgeoms(gdf, name="gauges") self.logger.info(f"Gauges map based on catchment river outlets added.") - self._setup_config_timeseries( - mapname="gauges", + self.setup_config_output_timeseries( + mapname="wflow_gauges", toml_output=toml_output, gauge_toml_header=gauge_toml_header, gauge_toml_param=gauge_toml_param, @@ -807,16 +817,19 @@ def setup_outlets( def setup_gauges( self, - gauges_fn="grdc", - source_gdf=None, - index_col=None, - snap_to_river=True, - mask=None, - derive_subcatch=False, - basename=None, - toml_output="csv", - gauge_toml_header=["Q", "P"], - gauge_toml_param=["lateral.river.q_av", "vertical.precipitation"], + gauges_fn: str, + index_col: Optional[str] = None, + snap_to_river: Optional[bool] = True, + mask: Optional[np.ndarray] = None, + max_dist: Optional[float] = 10e3, + derive_subcatch: Optional[bool] = False, + basename: Optional[str] = None, + toml_output: Optional[str] = "csv", + gauge_toml_header: Optional[List[str]] = ["Q", "P"], + gauge_toml_param: Optional[List[str]] = [ + "lateral.river.q_av", + "vertical.precipitation", + ], **kwargs, ): """This components sets a gauge map based on ``gauges_fn`` data. @@ -848,6 +861,9 @@ def setup_gauges( Snap point locations to the closest downstream river cell, by default True mask : np.boolean, optional If provided snaps to the mask, else snaps to the river (default). + max_dist : float, optional + Maximum distance between original and snapped point location. + A warning is logged if exceeded. By default 10 km. derive_subcatch : bool, optional Derive subcatch map for gauges, by default False basename : str, optional @@ -861,108 +877,170 @@ def setup_gauges( 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) and vertical.precipitation (for P). """ - # read existing staticgeoms; important to get the right basin when updating - self.staticgeoms - # TODO check snapping locations based on upstream area attribute of the gauge data - if gauges_fn is not None: - kwargs = {} - if isfile(gauges_fn): - # try to get epsg number directly, important when writting back data_catalog - if hasattr(self.crs, "to_epsg"): - code = self.crs.to_epsg() - else: - code = self.crs - kwargs.update(crs=code) - gdf = self.data_catalog.get_geodataframe( - gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs - ) - gdf = gdf.to_crs(self.crs) - elif source_gdf is not None and basename is None: - raise ValueError( - "Basename is required when setting gauges based on source_gdf" - ) - elif source_gdf is not None: - self.logger.info(f"Gauges locations read from source_gdf") - gdf = source_gdf.to_crs(self.crs) - else: - self.logger.warning( - "Either gauges_fn or source_gdf should be provided, skipping setup_gauges." - ) - + # Read data + kwargs = {} + if isfile(gauges_fn): + # try to get epsg number directly, important when writting back data_catalog + if hasattr(self.crs, "to_epsg"): + code = self.crs.to_epsg() + else: + code = self.crs + kwargs.update(crs=code) + gdf = self.data_catalog.get_geodataframe( + gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs + ) + # Create basename + if basename is None: + basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") + # Set index to index_col + if index_col is not None and index_col in gdf: + gdf = gdf.set_index(index_col) + + # Create gauge map, subcatch map and update toml if gdf.index.size == 0: self.logger.warning(f"No {gauges_fn} gauge locations found within domain") else: - if basename is None: - basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") self.logger.info( f"{gdf.index.size} {basename} gauge locations found within domain" ) - # Set index to index_col - if index_col is not None and index_col in gdf: - gdf = gdf.set_index(index_col) - xs, ys = np.vectorize(lambda p: (p.xy[0][0], p.xy[1][0]))(gdf["geometry"]) - idxs = self.staticmaps.raster.xy_to_idx(xs, ys) - ids = gdf.index.values - - if snap_to_river and mask is None: - mask = self.staticmaps[self._MAPS["rivmsk"]].values - da, idxs, ids = flw.gauge_map( - self.staticmaps, - idxs=idxs, - ids=ids, - stream=mask, - flwdir=self.flwdir, - logger=self.logger, - ) - # Filter gauges that could not be snapped to rivers - if snap_to_river: - ids_old = ids.copy() - da = da.where( - self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata - ) - ids_new = np.unique(da.values[da.values > 0]) - idxs = idxs[np.isin(ids_old, ids_new)] - ids = da.values.flat[idxs] - # Add to staticmaps - mapname = f'{str(self._MAPS["gauges"])}_{basename}' - self.set_staticmaps(da, name=mapname) - - # geoms - points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs)) - # if csv contains additional columns, these are also written in the staticgeoms - gdf_snapped = gpd.GeoDataFrame( - index=ids.astype(np.int32), geometry=points, crs=self.crs - ) - # Set the index name of gdf snapped based on original gdf - if gdf.index.name is not None: - gdf_snapped.index.name = gdf.index.name - else: - gdf_snapped.index.name = "fid" - gdf.index.name = "fid" - # Add gdf attributes to gdf_snapped (filter on snapped index before merging) - df_attrs = pd.DataFrame(gdf.drop(columns="geometry")) - df_attrs = df_attrs[np.isin(df_attrs.index, gdf_snapped.index)] - gdf_snapped = gdf_snapped.merge(df_attrs, how="inner", on=gdf.index.name) - # Add gdf_snapped to staticgeoms - self.set_staticgeoms(gdf_snapped, name=mapname.replace("wflow_", "")) - - self._setup_config_timeseries( - mapname=mapname, + self.create_gauges( + gdf_gauges=gdf, + basename=basename, + snap_to_river=snap_to_river, + mask=mask, + max_dist=max_dist, + derive_subcatch=derive_subcatch, toml_output=toml_output, gauge_toml_header=gauge_toml_header, gauge_toml_param=gauge_toml_param, ) - # add subcatch - if derive_subcatch: - da_basins = flw.basin_map( - self.staticmaps, self.flwdir, idxs=idxs, ids=ids - )[0] - mapname = self._MAPS["basins"] + "_" + basename - self.set_staticmaps(da_basins, name=mapname) - gdf_basins = self.staticmaps[mapname].raster.vectorize() - self.set_staticgeoms(gdf_basins, name=mapname.replace("wflow_", "")) + def create_gauges( + self, + gdf_gauges: gpd.GeoDataFrame, + basename: str, + snap_to_river: Optional[bool] = True, + mask: Optional[np.ndarray] = None, + max_dist: Optional[float] = 10e3, + derive_subcatch: Optional[bool] = False, + toml_output: Optional[str] = "csv", + gauge_toml_header: Optional[List[str]] = ["Q", "P"], + gauge_toml_param: Optional[List[str]] = [ + "lateral.river.q_av", + "vertical.precipitation", + ], + ): + """This components sets a gauge map based on ``gdf_gauges`` geodataframe. + + The gdf index will be used as IDs in the map. If ``snap_to_river`` is set + to True, the gauge location will be snapped to the boolean river mask. If + ``derive_subcatch`` is set to True, an additional subcatch map is derived from + the gauge locations. + + Adds model layers: + + * **wflow_gauges_source** map: gauge IDs map from source [-] (if gauges_fn) + * **wflow_subcatch_source** map: subcatchment based on gauge locations [-] (if derive_subcatch) + * **gauges_source** geom: polygon of gauges from source + * **subcatch_source** geom: polygon of subcatchment based on gauge locations [-] (if derive_subcatch) + + Parameters + ---------- + source_gdf : geopandas.GeoDataFame, optional + Direct gauges file geometry, by default None. + basename : str, optional + Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. + snap_to_river : bool, optional + Snap point locations to the closest downstream river cell, by default True + mask : np.boolean, optional + If provided snaps to the mask, else snaps to the river (default). + max_dist : float, optional + Maximum distance between original and snapped point location. + A warning is logged if exceeded. By default 10 km. + derive_subcatch : bool, optional + Derive subcatch map for gauges, by default False + 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) and P (for vertical.precipitation). + 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) and vertical.precipitation (for P). + """ + # read existing staticgeoms; important to get the right basin when updating + self.staticgeoms + # Reproject to model crs + gdf_gauges = gdf_gauges.to_crs(self.crs) + + # Get coords, index and ID + xs, ys = np.vectorize(lambda p: (p.xy[0][0], p.xy[1][0]))( + gdf_gauges["geometry"] + ) + idxs = self.staticmaps.raster.xy_to_idx(xs, ys) + ids = gdf_gauges.index.values + + # if snap_to_river use river map as the mask + if snap_to_river and mask is None: + mask = self.staticmaps[self._MAPS["rivmsk"]].values + # Derive gauge map + da, idxs, ids = flw.gauge_map( + self.staticmaps, + idxs=idxs, + ids=ids, + stream=mask, + flwdir=self.flwdir, + max_dist=max_dist, + logger=self.logger, + ) + # Filter gauges that could not be snapped to rivers + if snap_to_river: + ids_old = ids.copy() + da = da.where(self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata) + ids_new = np.unique(da.values[da.values > 0]) + idxs = idxs[np.isin(ids_old, ids_new)] + ids = da.values.flat[idxs] + # Add to staticmaps + mapname = f'{str(self._MAPS["gauges"])}_{basename}' + self.set_staticmaps(da, name=mapname) + + # geoms + points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs)) + # if csv contains additional columns, these are also written in the staticgeoms + gdf_snapped = gpd.GeoDataFrame( + index=ids.astype(np.int32), geometry=points, crs=self.crs + ) + # Set the index name of gdf snapped based on original gdf + if gdf_gauges.index.name is not None: + gdf_snapped.index.name = gdf_gauges.index.name + else: + gdf_snapped.index.name = "fid" + gdf_gauges.index.name = "fid" + # Add gdf attributes to gdf_snapped (filter on snapped index before merging) + df_attrs = pd.DataFrame(gdf_gauges.drop(columns="geometry")) + df_attrs = df_attrs[np.isin(df_attrs.index, gdf_snapped.index)] + gdf_snapped = gdf_snapped.merge(df_attrs, how="inner", on=gdf_gauges.index.name) + # Add gdf_snapped to staticgeoms + self.set_staticgeoms(gdf_snapped, name=mapname.replace("wflow_", "")) + + # Add output timeseries for gauges in the toml + self.setup_config_output_timeseries( + mapname=mapname, + toml_output=toml_output, + gauge_toml_header=gauge_toml_header, + gauge_toml_param=gauge_toml_param, + ) + + # add subcatch + if derive_subcatch: + da_basins = flw.basin_map(self.staticmaps, self.flwdir, idxs=idxs, ids=ids)[ + 0 + ] + mapname = self._MAPS["basins"] + "_" + basename + self.set_staticmaps(da_basins, name=mapname) + gdf_basins = self.staticmaps[mapname].raster.vectorize() + self.set_staticgeoms(gdf_basins, name=mapname.replace("wflow_", "")) def setup_areamap( self, From c5e9ad993830074b2bc88d850f9b2d51316ef6ce Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 30 Mar 2023 00:19:04 +0800 Subject: [PATCH 04/17] support gauges geodataset fixes #85 --- hydromt_wflow/wflow.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 4f19e013..d85665d7 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -887,9 +887,23 @@ def setup_gauges( else: code = self.crs kwargs.update(crs=code) - gdf = self.data_catalog.get_geodataframe( - gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs - ) + gdf = self.data_catalog.get_geodataframe( + gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs + ) + elif self.data_catalog[gauges_fn]["data_type"] == "GeoDataFrame": + gdf = self.data_catalog.get_geodataframe( + gauges_fn, geom=self.basins, **kwargs + ) + elif self.data_catalog[gauges_fn]["data_type"] == "GeoDataset": + da = self.data_catalog.get_geodataset(gauges_fn, geom=self.basins, **kwargs) + gdf = da.vector.to_gdf() + else: + raise ValueError( + f"{gauges_fn} data source not found or incorrect data_type (GeoDataFrame or GeoDataset)." + ) + # Check for point geometry + if not np.all(np.isin(gdf.geometry.type, "Point")): + raise ValueError(f"{gauges_fn} contains other geometries than Point") # Create basename if basename is None: basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") @@ -947,8 +961,8 @@ def create_gauges( Parameters ---------- - source_gdf : geopandas.GeoDataFame, optional - Direct gauges file geometry, by default None. + source_gdf : geopandas.GeoDataFame + Gauges GeoDataFrame geometry basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. snap_to_river : bool, optional From 1259aae382d33f69be15a565de6f2ca6980d6731 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 30 Mar 2023 17:43:08 +0800 Subject: [PATCH 05/17] fix tests after gauge changes --- hydromt_wflow/wflow.py | 23 +++++++++++++---------- tests/data/wflow_piave_build_subbasin.ini | 1 + 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index d85665d7..7a8a7f6e 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -750,9 +750,9 @@ def setup_config_output_timeseries( } 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) + # 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." @@ -797,7 +797,10 @@ def setup_outlets( (self.staticmaps[self._MAPS["rivmsk"]] > 0).values.flat[idxs_out] ] da_out, idxs_out, ids_out = flw.gauge_map( - self.staticmaps, idxs=self.flwdir.idxs_out + self.staticmaps, + idxs=idxs_out, + flwdir=self.flwdir, + logger=self.logger, ) self.set_staticmaps(da_out, name=self._MAPS["gauges"]) points = gpd.points_from_xy(*self.staticmaps.raster.idx_to_xy(idxs_out)) @@ -811,8 +814,8 @@ def setup_outlets( self.setup_config_output_timeseries( mapname="wflow_gauges", toml_output=toml_output, - gauge_toml_header=gauge_toml_header, - gauge_toml_param=gauge_toml_param, + header=gauge_toml_header, + param=gauge_toml_param, ) def setup_gauges( @@ -890,11 +893,11 @@ def setup_gauges( gdf = self.data_catalog.get_geodataframe( gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs ) - elif self.data_catalog[gauges_fn]["data_type"] == "GeoDataFrame": + elif self.data_catalog[gauges_fn].data_type == "GeoDataFrame": gdf = self.data_catalog.get_geodataframe( gauges_fn, geom=self.basins, **kwargs ) - elif self.data_catalog[gauges_fn]["data_type"] == "GeoDataset": + elif self.data_catalog[gauges_fn].data_type == "GeoDataset": da = self.data_catalog.get_geodataset(gauges_fn, geom=self.basins, **kwargs) gdf = da.vector.to_gdf() else: @@ -1042,8 +1045,8 @@ def create_gauges( self.setup_config_output_timeseries( mapname=mapname, toml_output=toml_output, - gauge_toml_header=gauge_toml_header, - gauge_toml_param=gauge_toml_param, + header=gauge_toml_header, + param=gauge_toml_param, ) # add subcatch diff --git a/tests/data/wflow_piave_build_subbasin.ini b/tests/data/wflow_piave_build_subbasin.ini index 9c105fbf..335875fe 100644 --- a/tests/data/wflow_piave_build_subbasin.ini +++ b/tests/data/wflow_piave_build_subbasin.ini @@ -1,4 +1,5 @@ [setup_config] +output.lateral.land.h = "h_land" [setup_basemaps] res = 0.016666666666666666 From db829b5352df7e7b690b966726700bedad56c928 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 30 Mar 2023 17:43:36 +0800 Subject: [PATCH 06/17] update docs and changelog --- docs/api.rst | 4 ++++ docs/changelog.rst | 10 +++++++--- docs/user_guide/sediment_model_setup.rst | 8 ++++++-- docs/user_guide/wflow_model_setup.rst | 6 +++++- 4 files changed, 22 insertions(+), 6 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 6a300711..76cdfccb 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -37,8 +37,10 @@ Setup components WflowModel.setup_laimaps WflowModel.setup_soilmaps WflowModel.setup_hydrodem + WflowModel.setup_outlets WflowModel.setup_gauges WflowModel.setup_areamap + WflowModel.setup_config_output_timeseries WflowModel.setup_precip_forcing WflowModel.setup_temp_pet_forcing WflowModel.setup_constant_pars @@ -146,8 +148,10 @@ Setup components WflowSedimentModel.setup_soilmaps WflowSedimentModel.setup_riverwidth WflowSedimentModel.setup_riverbedsed + WflowSedimentModel.setup_outlets WflowSedimentModel.setup_gauges WflowSedimentModel.setup_areamap + WflowSedimentModel.setup_config_output_timeseries WflowSedimentModel.setup_constant_pars WflowSedimentModel.setup_staticmaps_from_raster diff --git a/docs/changelog.rst b/docs/changelog.rst index b13f6522..fee195f0 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -11,15 +11,22 @@ Unreleased Added ----- +- Add options to calculate daily Penman-Monteith potential evaporation using the pyet package. Depending on the available variables, two options are defined ``penman-monteith_tdew`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind_u', 'wind_v', 'temp_dew', 'kin', 'press_msl']) and ``penman-monteith_rh_simple`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind', 'rh', 'kin']). - Support in toml for dir_input and dir_output options. `PR #140 `_ - In **setup_reservoirs**: Global Water Watch compatibility for determining reservoir parameters. - In **setup_reservoirs**: All dowloaded reservoir timeseries are saved to root in 1 csv file. Column headers indicate reservoir id. +- **setup_oulets**: Add map/geom of basin outlets (on river or all) and optionally updates outputs in toml file. +- **setup_config_output_timeseries**: add new variable/column to the netcf/csv output section of the toml based on a selected gauge/area map. +- **setup_gauges**: support for snapping based on a user defined max distance and snapping based on upstream area attribute. +- **setup_gauges**: gauges_fn can be both GeoDataFrame or GeoDataset (new) data_type. +- **create_gauges**: submethod to add gauges directly from a geopandas.geodataframe object rather than a data source from catalog. Changed ------- - Default tomls are now using the dir_output option to specify *run_default* folder. - in **setup_reservoirs**: options 'usehe' and 'priorityjrc' are removed and replaced with 'timeseries_fn'. Options are ['jrc', 'gww']. By default None to use reservoir_fn data directly. - in **setup_areamap**: name of the added map is based on column name of the vector data (col2raster) instead of name of the vector data file (area_fn). Allows to add several maps from one vector data file. +- in **setup_gauges**: basins outlets are not generated by this function anymore but by the new **setup_outlets** method. Fixed ----- @@ -33,9 +40,6 @@ Fixed Deprecated ---------- -Added -^^^^^ -- add options to calculate daily Penman-Monteith potential evaporation using the pyet package. Depending on the available variables, two options are defined ``penman-monteith_tdew`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind_u', 'wind_v', 'temp_dew', 'kin', 'press_msl']) and ``penman-monteith_rh_simple`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind', 'rh', 'kin']). v0.2.1 (22 November 2022) ========================= diff --git a/docs/user_guide/sediment_model_setup.rst b/docs/user_guide/sediment_model_setup.rst index 8e3ffc7b..eb7cac84 100644 --- a/docs/user_guide/sediment_model_setup.rst +++ b/docs/user_guide/sediment_model_setup.rst @@ -52,10 +52,14 @@ a specific method see its documentation. - This component sets the river width parameter based on a power-lay relationship with a predictor. * - :py:func:`~WflowSedimentModel.setup_riverbedsed` - Setup sediments based river bed characteristics maps. - * - :py:func:`~WflowSedimentModel.setup_gauges` - - This components sets the default gauge map based on basin outlets and additional gauge maps based on gauges_fn data. + * - :py:func:`~WflowModel.setup_outlets` + - This method sets the default gauge map based on basin outlets. + * - :py:func:`~WflowModel.setup_gauges` + - This method sets the default gauge map based on a gauges_fn data. * - :py:func:`~WflowModel.setup_areamap` - Setup area map from vector data to save wflow outputs for specific area. + * - :py:func:`~WflowModel.setup_config_output_timeseries` + - This method add a new variable/column to the netcf/csv output section of the toml based on a selected gauge/area map. * - :py:func:`~WflowSedimentModel.setup_constant_pars` - Setup constant parameter maps. * - :py:func:`~WflowModel.setup_staticmaps_from_raster` diff --git a/docs/user_guide/wflow_model_setup.rst b/docs/user_guide/wflow_model_setup.rst index 100de4b6..8cb43f31 100644 --- a/docs/user_guide/wflow_model_setup.rst +++ b/docs/user_guide/wflow_model_setup.rst @@ -52,10 +52,14 @@ a specific method see its documentation. - This component derives several (layered) soil parameters based on a database with physical soil properties using available point-scale (pedo)transfer functions (PTFs) from literature with upscaling rulesto ensure flux matching across scales. * - :py:func:`~WflowModel.setup_hydrodem` - This component adds a hydrologically conditioned elevation (hydrodem) map for river and/or land local-inertial routing. + * - :py:func:`~WflowModel.setup_outlets` + - This method sets the default gauge map based on basin outlets. * - :py:func:`~WflowModel.setup_gauges` - - This method sets the default gauge map based on basin outlets and additional gauge maps based on gauges_fn data. + - This method sets the default gauge map based on a gauges_fn data. * - :py:func:`~WflowModel.setup_areamap` - Setup area map from vector data to save wflow outputs for specific area. + * - :py:func:`~WflowModel.setup_config_output_timeseries` + - This method add new variable/column to the netcf/csv output section of the toml based on a selected gauge/area map. * - :py:func:`~WflowModel.setup_precip_forcing` - Setup gridded precipitation forcing at model resolution. * - :py:func:`~WflowModel.setup_temp_pet_forcing` From 79eb1369c84c4e3e09e69da6f045f7220c97be3e Mon Sep 17 00:00:00 2001 From: hboisgon Date: Fri, 31 Mar 2023 23:45:53 +0800 Subject: [PATCH 07/17] add gauge snapping based on upstream area fixes #170 --- hydromt_wflow/wflow.py | 97 +++++++++++++++++++++-------- hydromt_wflow/workflows/__init__.py | 1 + hydromt_wflow/workflows/gauges.py | 61 ++++++++++++++++++ tests/test_model_methods.py | 24 +++++++ 4 files changed, 156 insertions(+), 27 deletions(-) create mode 100644 hydromt_wflow/workflows/gauges.py diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 7a8a7f6e..5fb1729d 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -824,7 +824,10 @@ def setup_gauges( index_col: Optional[str] = None, snap_to_river: Optional[bool] = True, mask: Optional[np.ndarray] = None, + snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, + wdw: Optional[int] = 5, + rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, basename: Optional[str] = None, toml_output: Optional[str] = "csv", @@ -856,17 +859,26 @@ def setup_gauges( ---------- gauges_fn : str, {"grdc"}, optional Known source name or path to gauges file geometry file, by default None. - source_gdf : geopandas.GeoDataFame, optional - Direct gauges file geometry, by default None. + + * Required variables if snap_uparea is True: ["uparea"] index_col : str, optional Column in gauges_fn to use for ID values, by default None (use the default index column) - snap_to_river : bool, optional - Snap point locations to the closest downstream river cell, by default True mask : np.boolean, optional If provided snaps to the mask, else snaps to the river (default). + snap_to_river : bool, optional + Snap point locations to the closest downstream river cell, by default True + snap_uparea: bool, optional + Snap gauges based on upstream area. Gauges_fn should have "uparea" in its attributes. max_dist : float, optional Maximum distance between original and snapped point location. A warning is logged if exceeded. By default 10 km. + wdw: int, optional + Window size in number of cells around discharge boundary locations + to snap to, only used if ``snap_uparea`` is True. By default 5. + rel_error: float, optional + Maximum relative error (default 0.05) + between the gauge location upstream area and the upstream area of + the best fit grid cell, only used if snap_area is True. derive_subcatch : bool, optional Derive subcatch map for gauges, by default False basename : str, optional @@ -895,18 +907,19 @@ def setup_gauges( ) elif self.data_catalog[gauges_fn].data_type == "GeoDataFrame": gdf = self.data_catalog.get_geodataframe( - gauges_fn, geom=self.basins, **kwargs + gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs ) elif self.data_catalog[gauges_fn].data_type == "GeoDataset": da = self.data_catalog.get_geodataset(gauges_fn, geom=self.basins, **kwargs) gdf = da.vector.to_gdf() + # Check for point geometry + if not np.all(np.isin(gdf.geometry.type, "Point")): + raise ValueError(f"{gauges_fn} contains other geometries than Point") else: raise ValueError( f"{gauges_fn} data source not found or incorrect data_type (GeoDataFrame or GeoDataset)." ) - # Check for point geometry - if not np.all(np.isin(gdf.geometry.type, "Point")): - raise ValueError(f"{gauges_fn} contains other geometries than Point") + # Create basename if basename is None: basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") @@ -926,7 +939,10 @@ def setup_gauges( basename=basename, snap_to_river=snap_to_river, mask=mask, + snap_uparea=snap_uparea, max_dist=max_dist, + wdw=wdw, + rel_error=rel_error, derive_subcatch=derive_subcatch, toml_output=toml_output, gauge_toml_header=gauge_toml_header, @@ -938,8 +954,11 @@ def create_gauges( gdf_gauges: gpd.GeoDataFrame, basename: str, snap_to_river: Optional[bool] = True, - mask: Optional[np.ndarray] = None, + mask: Optional[str] = None, + snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, + wdw: Optional[int] = 5, + rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, toml_output: Optional[str] = "csv", gauge_toml_header: Optional[List[str]] = ["Q", "P"], @@ -972,9 +991,18 @@ def create_gauges( Snap point locations to the closest downstream river cell, by default True mask : np.boolean, optional If provided snaps to the mask, else snaps to the river (default). + snap_uparea: bool, optional + Snap gauges based on upstream area. Gauges_fn should have "uparea" in its attributes. max_dist : float, optional Maximum distance between original and snapped point location. A warning is logged if exceeded. By default 10 km. + wdw: int, optional + Window size in number of cells around discharge boundary locations + to snap to, only used if ``snap_uparea`` is True. By default 5. + rel_error: float, optional + Maximum relative error (default 0.05) + between the gauge location upstream area and the upstream area of + the best fit grid cell, only used if snap_area is True. derive_subcatch : bool, optional Derive subcatch map for gauges, by default False toml_output : str, optional @@ -1000,24 +1028,39 @@ def create_gauges( # if snap_to_river use river map as the mask if snap_to_river and mask is None: - mask = self.staticmaps[self._MAPS["rivmsk"]].values - # Derive gauge map - da, idxs, ids = flw.gauge_map( - self.staticmaps, - idxs=idxs, - ids=ids, - stream=mask, - flwdir=self.flwdir, - max_dist=max_dist, - logger=self.logger, - ) - # Filter gauges that could not be snapped to rivers - if snap_to_river: - ids_old = ids.copy() - da = da.where(self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata) - ids_new = np.unique(da.values[da.values > 0]) - idxs = idxs[np.isin(ids_old, ids_new)] - ids = da.values.flat[idxs] + mask = self._MAPS["rivmsk"] + if mask is not None: + mask = self.staticmaps[mask].values + if mask is not None: + # Derive gauge map + da, idxs, ids = flw.gauge_map( + self.staticmaps, + idxs=idxs, + ids=ids, + stream=mask, + flwdir=self.flwdir, + max_dist=max_dist, + logger=self.logger, + ) + # Filter gauges that could not be snapped to rivers + if snap_to_river: + ids_old = ids.copy() + da = da.where( + self.staticmaps[self._MAPS["rivmsk"]] != 0, da.raster.nodata + ) + ids_new = np.unique(da.values[da.values > 0]) + idxs = idxs[np.isin(ids_old, ids_new)] + ids = da.values.flat[idxs] + elif snap_uparea and "uparea" in gdf_gauges.columns: + # Derive gauge map based on upstream area snapping + da, idxs, ids = workflows.gauge_map_uparea( + self.staticmaps, + gdf_gauges, + uparea_name="wflow_uparea", + wdw=wdw, + rel_error=rel_error, + logger=self.logger, + ) # Add to staticmaps mapname = f'{str(self._MAPS["gauges"])}_{basename}' self.set_staticmaps(da, name=mapname) diff --git a/hydromt_wflow/workflows/__init__.py b/hydromt_wflow/workflows/__init__.py index a7068912..ef67e4dd 100644 --- a/hydromt_wflow/workflows/__init__.py +++ b/hydromt_wflow/workflows/__init__.py @@ -6,3 +6,4 @@ from .landuse import * from .soilgrids import * from .glaciers import * +from .gauges import * diff --git a/hydromt_wflow/workflows/gauges.py b/hydromt_wflow/workflows/gauges.py new file mode 100644 index 00000000..a6fce47a --- /dev/null +++ b/hydromt_wflow/workflows/gauges.py @@ -0,0 +1,61 @@ +import numpy as np +import xarray as xr +import geopandas as gpd +from hydromt import flw +from typing import Optional +import logging + + +logger = logging.getLogger(__name__) + + +__all__ = ["gauge_map_uparea"] + + +def gauge_map_uparea( + ds: xr.Dataset, + gdf: gpd.GeoDataFrame, + uparea_name: Optional[str] = "wflow_uparea", + wdw: Optional[int] = 1, + rel_error: float = 0.05, + abs_error: float = 50, + logger=logger, +): + """ """ + if uparea_name not in ds: + return + + ds_wdw = ds.raster.sample(gdf, wdw=wdw) + logger.debug( + f"Snapping gauges points to best matching uparea cell within wdw (size={wdw})." + ) + upa0 = xr.DataArray(gdf["uparea"], dims=("index")) + upa_dff = np.abs(ds_wdw[uparea_name].where(ds_wdw[uparea_name] > 0).load() - upa0) + upa_check = np.logical_or((upa_dff / upa0) <= rel_error, upa_dff <= abs_error) + # find best matching uparea cell in window + i_wdw = upa_dff.argmin("wdw").load() + + idx_valid = np.where(upa_check.isel(wdw=i_wdw).values)[0] + if idx_valid.size < gdf.index.size: + logger.warning( + f"{idx_valid.size}/{gdf.index.size} gauge points successfully snapped." + ) + i_wdw = i_wdw.isel(index=idx_valid) + ds_out = ds_wdw.isel(wdw=i_wdw.load(), index=idx_valid) + + idxs_out = ds.raster.xy_to_idx( + xs=ds_out[ds.raster.x_dim].values, ys=ds_out[ds.raster.y_dim].values + ) + ids_out = ds_out.index.values + + # Derive gauge map + da, idxs_out, ids_out = flw.gauge_map( + ds, + idxs=idxs_out, + ids=ids_out, + stream=None, + flwdir=None, + logger=logger, + ) + + return da, idxs_out, ids_out diff --git a/tests/test_model_methods.py b/tests/test_model_methods.py index 2a3d19e1..c964cfb5 100644 --- a/tests/test_model_methods.py +++ b/tests/test_model_methods.py @@ -110,3 +110,27 @@ def test_setup_reservoirs(source, tmpdir): ) == number_of_reservoirs ), f"Number of non-null values in {i} not equal to number of reservoirs in model area" + + +def test_setup_gauges(): + logger = logging.getLogger(__name__) + # read model from examples folder + root = join(EXAMPLEDIR, "wflow_piave_subbasin") + + # Initialize model and read results + mod = WflowModel(root=root, mode="r", data_libs="artifact_data", logger=logger) + # uparea rename not in the latest artifact_data version + mod.data_catalog["grdc"].rename = {"area": "uparea"} + mod.setup_gauges( + gauges_fn="grdc", + basename="grdc_uparea", + snap_to_river=False, + mask=None, + snap_uparea=True, + wdw=5, + rel_error=0.05, + ) + gdf = mod.staticgeoms["gauges_grdc_uparea"] + ds_samp = mod.staticmaps[["wflow_river", "wflow_uparea"]].raster.sample(gdf, wdw=0) + assert np.all(ds_samp["wflow_river"].values == 1) + assert np.allclose(ds_samp["wflow_uparea"].values, gdf["uparea"].values, rtol=0.05) From 953b1e12c88ccf7814eee02a82d73ef26e6bf7f1 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 21 Apr 2023 17:25:47 +0200 Subject: [PATCH 08/17] fix gauges with index 0 --- hydromt_wflow/wflow.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 5fb1729d..1c05644b 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -826,9 +826,10 @@ def setup_gauges( mask: Optional[np.ndarray] = None, snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, - wdw: Optional[int] = 5, + wdw: Optional[int] = 2, rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, + index_col: Optional[str] = None, basename: Optional[str] = None, toml_output: Optional[str] = "csv", gauge_toml_header: Optional[List[str]] = ["Q", "P"], @@ -874,13 +875,15 @@ def setup_gauges( A warning is logged if exceeded. By default 10 km. wdw: int, optional Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 5. + to snap to, only used if ``snap_uparea`` is True. By default 2. rel_error: float, optional Maximum relative error (default 0.05) between the gauge location upstream area and the upstream area of the best fit grid cell, only used if snap_area is True. derive_subcatch : bool, optional Derive subcatch map for gauges, by default False + index_col : str, optional + Column in gauges_fn to use for ID values, by default None (use the default index column) basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. toml_output : str, optional @@ -937,6 +940,7 @@ def setup_gauges( self.create_gauges( gdf_gauges=gdf, basename=basename, + index_col=index_col, snap_to_river=snap_to_river, mask=mask, snap_uparea=snap_uparea, @@ -953,11 +957,12 @@ def create_gauges( self, gdf_gauges: gpd.GeoDataFrame, basename: str, + index_col: Optional[str] = None, snap_to_river: Optional[bool] = True, mask: Optional[str] = None, snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, - wdw: Optional[int] = 5, + wdw: Optional[int] = 2, rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, toml_output: Optional[str] = "csv", @@ -987,6 +992,8 @@ def create_gauges( Gauges GeoDataFrame geometry basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. + index_col : str, optional + Column name to use as index, if None use the source_gdf index. snap_to_river : bool, optional Snap point locations to the closest downstream river cell, by default True mask : np.boolean, optional @@ -998,7 +1005,7 @@ def create_gauges( A warning is logged if exceeded. By default 10 km. wdw: int, optional Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 5. + to snap to, only used if ``snap_uparea`` is True. By default 2. rel_error: float, optional Maximum relative error (default 0.05) between the gauge location upstream area and the upstream area of @@ -1024,7 +1031,13 @@ def create_gauges( gdf_gauges["geometry"] ) idxs = self.staticmaps.raster.xy_to_idx(xs, ys) - ids = gdf_gauges.index.values + if index_col is not None and index_col in gdf_gauges.columns: + ids = gdf_gauges[index_col].values + else: + ids = gdf_gauges.index.values + if np.any(ids == 0): + ids = ids + 1 + self.logger.warning("Gauge ID 0 is not allowed, setting to 1") # if snap_to_river use river map as the mask if snap_to_river and mask is None: From 497b72f8fd84afc33444d396738ae4e07032dcaa Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Fri, 21 Apr 2023 17:42:00 +0200 Subject: [PATCH 09/17] fix gauges with index 0 --- hydromt_wflow/wflow.py | 25 +++++++++---------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 1c05644b..4e62bd1d 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -826,10 +826,9 @@ def setup_gauges( mask: Optional[np.ndarray] = None, snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, - wdw: Optional[int] = 2, + wdw: Optional[int] = 3, rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, - index_col: Optional[str] = None, basename: Optional[str] = None, toml_output: Optional[str] = "csv", gauge_toml_header: Optional[List[str]] = ["Q", "P"], @@ -875,15 +874,13 @@ def setup_gauges( A warning is logged if exceeded. By default 10 km. wdw: int, optional Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 2. + to snap to, only used if ``snap_uparea`` is True. By default 3. rel_error: float, optional Maximum relative error (default 0.05) between the gauge location upstream area and the upstream area of the best fit grid cell, only used if snap_area is True. derive_subcatch : bool, optional Derive subcatch map for gauges, by default False - index_col : str, optional - Column in gauges_fn to use for ID values, by default None (use the default index column) basename : str, optional Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. toml_output : str, optional @@ -926,9 +923,6 @@ def setup_gauges( # Create basename if basename is None: basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") - # Set index to index_col - if index_col is not None and index_col in gdf: - gdf = gdf.set_index(index_col) # Create gauge map, subcatch map and update toml if gdf.index.size == 0: @@ -962,7 +956,7 @@ def create_gauges( mask: Optional[str] = None, snap_uparea: Optional[bool] = False, max_dist: Optional[float] = 10e3, - wdw: Optional[int] = 2, + wdw: Optional[int] = 3, rel_error: Optional[float] = 0.05, derive_subcatch: Optional[bool] = False, toml_output: Optional[str] = "csv", @@ -1005,7 +999,7 @@ def create_gauges( A warning is logged if exceeded. By default 10 km. wdw: int, optional Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 2. + to snap to, only used if ``snap_uparea`` is True. By default 3. rel_error: float, optional Maximum relative error (default 0.05) between the gauge location upstream area and the upstream area of @@ -1024,7 +1018,7 @@ def create_gauges( # read existing staticgeoms; important to get the right basin when updating self.staticgeoms # Reproject to model crs - gdf_gauges = gdf_gauges.to_crs(self.crs) + gdf_gauges = gdf_gauges.to_crs(self.crs).copy() # Get coords, index and ID xs, ys = np.vectorize(lambda p: (p.xy[0][0], p.xy[1][0]))( @@ -1032,12 +1026,11 @@ def create_gauges( ) idxs = self.staticmaps.raster.xy_to_idx(xs, ys) if index_col is not None and index_col in gdf_gauges.columns: - ids = gdf_gauges[index_col].values - else: - ids = gdf_gauges.index.values - if np.any(ids == 0): - ids = ids + 1 + gdf_gauges = gdf_gauges.set_index(index_col) + if np.any(gdf_gauges.index == 0): self.logger.warning("Gauge ID 0 is not allowed, setting to 1") + gdf_gauges.index = gdf_gauges.index.values + 1 + ids = gdf_gauges.index.values # if snap_to_river use river map as the mask if snap_to_river and mask is None: From f27a35695fa303dfb264fb0c5ac37736e756379c Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Mon, 24 Apr 2023 14:20:04 +0200 Subject: [PATCH 10/17] fix indexing with dask bool array error --- hydromt_wflow/workflows/waterbodies.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hydromt_wflow/workflows/waterbodies.py b/hydromt_wflow/workflows/waterbodies.py index a57c689b..e1704abd 100644 --- a/hydromt_wflow/workflows/waterbodies.py +++ b/hydromt_wflow/workflows/waterbodies.py @@ -79,8 +79,8 @@ def waterbodymaps( xdim = ds_like.raster.x_dim for i in res_id: res_acc = ds_like[uparea_name].where(ds_out["resareas"] == i) - # max_res_acc = np.amax(res_acc.values()) - max_res_acc = res_acc.where(res_acc == res_acc.max(), drop=True).squeeze() + max_res_acc = np.amax(res_acc.values()) + max_res_acc = res_acc.where(res_acc == max_res_acc, drop=True).squeeze() yacc = max_res_acc[ydim].values xacc = max_res_acc[xdim].values ds_out["reslocs"].loc[{f"{ydim}": yacc, f"{xdim}": xacc}] = i From 0b25e2473419c22045c0d5ffb9989b76ca16180e Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Mon, 24 Apr 2023 14:26:56 +0200 Subject: [PATCH 11/17] fix indexing with dask bool array error --- hydromt_wflow/workflows/waterbodies.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hydromt_wflow/workflows/waterbodies.py b/hydromt_wflow/workflows/waterbodies.py index e1704abd..428f1eb1 100644 --- a/hydromt_wflow/workflows/waterbodies.py +++ b/hydromt_wflow/workflows/waterbodies.py @@ -79,8 +79,9 @@ def waterbodymaps( xdim = ds_like.raster.x_dim for i in res_id: res_acc = ds_like[uparea_name].where(ds_out["resareas"] == i) - max_res_acc = np.amax(res_acc.values()) - max_res_acc = res_acc.where(res_acc == max_res_acc, drop=True).squeeze() + max_res_acc = res_acc.where( + res_acc == np.amax(res_acc.values), drop=True + ).squeeze() yacc = max_res_acc[ydim].values xacc = max_res_acc[xdim].values ds_out["reslocs"].loc[{f"{ydim}": yacc, f"{xdim}": xacc}] = i From 9175ea00ce7ab75a3f0096df0841f7de71109784 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Mon, 24 Apr 2023 14:29:42 +0200 Subject: [PATCH 12/17] fix indexing with dask bool array error --- hydromt_wflow/workflows/waterbodies.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/hydromt_wflow/workflows/waterbodies.py b/hydromt_wflow/workflows/waterbodies.py index 428f1eb1..614b41da 100644 --- a/hydromt_wflow/workflows/waterbodies.py +++ b/hydromt_wflow/workflows/waterbodies.py @@ -78,10 +78,8 @@ def waterbodymaps( ydim = ds_like.raster.y_dim xdim = ds_like.raster.x_dim for i in res_id: - res_acc = ds_like[uparea_name].where(ds_out["resareas"] == i) - max_res_acc = res_acc.where( - res_acc == np.amax(res_acc.values), drop=True - ).squeeze() + res_acc = ds_like[uparea_name].where(ds_out["resareas"] == i).load() + max_res_acc = res_acc.where(res_acc == res_acc.max(), drop=True).squeeze() yacc = max_res_acc[ydim].values xacc = max_res_acc[xdim].values ds_out["reslocs"].loc[{f"{ydim}": yacc, f"{xdim}": xacc}] = i From 9b82083da9dfaad374ecc92c584d71b15a429d48 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Tue, 6 Jun 2023 09:57:04 +0800 Subject: [PATCH 13/17] remove create_gauges method --- docs/changelog.rst | 1 - hydromt_wflow/wflow.py | 124 +++++++----------------------------- tests/test_model_methods.py | 25 ++++++++ 3 files changed, 47 insertions(+), 103 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index fee195f0..a9b5df8e 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -19,7 +19,6 @@ Added - **setup_config_output_timeseries**: add new variable/column to the netcf/csv output section of the toml based on a selected gauge/area map. - **setup_gauges**: support for snapping based on a user defined max distance and snapping based on upstream area attribute. - **setup_gauges**: gauges_fn can be both GeoDataFrame or GeoDataset (new) data_type. -- **create_gauges**: submethod to add gauges directly from a geopandas.geodataframe object rather than a data source from catalog. Changed ------- diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 4e62bd1d..19ba0973 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -892,7 +892,6 @@ def setup_gauges( 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) and vertical.precipitation (for P). """ - # TODO check snapping locations based on upstream area attribute of the gauge data # Read data kwargs = {} if isfile(gauges_fn): @@ -902,18 +901,18 @@ def setup_gauges( else: code = self.crs kwargs.update(crs=code) - gdf = self.data_catalog.get_geodataframe( + gdf_gauges = self.data_catalog.get_geodataframe( gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs ) elif self.data_catalog[gauges_fn].data_type == "GeoDataFrame": - gdf = self.data_catalog.get_geodataframe( + gdf_gauges = self.data_catalog.get_geodataframe( gauges_fn, geom=self.basins, assert_gtype="Point", **kwargs ) elif self.data_catalog[gauges_fn].data_type == "GeoDataset": da = self.data_catalog.get_geodataset(gauges_fn, geom=self.basins, **kwargs) - gdf = da.vector.to_gdf() + gdf_gauges = da.vector.to_gdf() # Check for point geometry - if not np.all(np.isin(gdf.geometry.type, "Point")): + if not np.all(np.isin(gdf_gauges.geometry.type, "Point")): raise ValueError(f"{gauges_fn} contains other geometries than Point") else: raise ValueError( @@ -925,96 +924,16 @@ def setup_gauges( basename = os.path.basename(gauges_fn).split(".")[0].replace("_", "-") # Create gauge map, subcatch map and update toml - if gdf.index.size == 0: + if gdf_gauges.index.size == 0: self.logger.warning(f"No {gauges_fn} gauge locations found within domain") - else: - self.logger.info( - f"{gdf.index.size} {basename} gauge locations found within domain" - ) - self.create_gauges( - gdf_gauges=gdf, - basename=basename, - index_col=index_col, - snap_to_river=snap_to_river, - mask=mask, - snap_uparea=snap_uparea, - max_dist=max_dist, - wdw=wdw, - rel_error=rel_error, - derive_subcatch=derive_subcatch, - toml_output=toml_output, - gauge_toml_header=gauge_toml_header, - gauge_toml_param=gauge_toml_param, - ) - - def create_gauges( - self, - gdf_gauges: gpd.GeoDataFrame, - basename: str, - index_col: Optional[str] = None, - snap_to_river: Optional[bool] = True, - mask: Optional[str] = None, - snap_uparea: Optional[bool] = False, - max_dist: Optional[float] = 10e3, - wdw: Optional[int] = 3, - rel_error: Optional[float] = 0.05, - derive_subcatch: Optional[bool] = False, - toml_output: Optional[str] = "csv", - gauge_toml_header: Optional[List[str]] = ["Q", "P"], - gauge_toml_param: Optional[List[str]] = [ - "lateral.river.q_av", - "vertical.precipitation", - ], - ): - """This components sets a gauge map based on ``gdf_gauges`` geodataframe. - - The gdf index will be used as IDs in the map. If ``snap_to_river`` is set - to True, the gauge location will be snapped to the boolean river mask. If - ``derive_subcatch`` is set to True, an additional subcatch map is derived from - the gauge locations. - Adds model layers: + return - * **wflow_gauges_source** map: gauge IDs map from source [-] (if gauges_fn) - * **wflow_subcatch_source** map: subcatchment based on gauge locations [-] (if derive_subcatch) - * **gauges_source** geom: polygon of gauges from source - * **subcatch_source** geom: polygon of subcatchment based on gauge locations [-] (if derive_subcatch) + # Create the gauges map + self.logger.info( + f"{gdf_gauges.index.size} {basename} gauge locations found within domain" + ) - Parameters - ---------- - source_gdf : geopandas.GeoDataFame - Gauges GeoDataFrame geometry - basename : str, optional - Map name in staticmaps (wflow_gauges_basename), if None use the gauges_fn basename. - index_col : str, optional - Column name to use as index, if None use the source_gdf index. - snap_to_river : bool, optional - Snap point locations to the closest downstream river cell, by default True - mask : np.boolean, optional - If provided snaps to the mask, else snaps to the river (default). - snap_uparea: bool, optional - Snap gauges based on upstream area. Gauges_fn should have "uparea" in its attributes. - max_dist : float, optional - Maximum distance between original and snapped point location. - A warning is logged if exceeded. By default 10 km. - wdw: int, optional - Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 3. - rel_error: float, optional - Maximum relative error (default 0.05) - between the gauge location upstream area and the upstream area of - the best fit grid cell, only used if snap_area is True. - derive_subcatch : bool, optional - Derive subcatch map for gauges, by default False - 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) and P (for vertical.precipitation). - 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) and vertical.precipitation (for P). - """ # read existing staticgeoms; important to get the right basin when updating self.staticgeoms # Reproject to model crs @@ -1037,7 +956,17 @@ def create_gauges( mask = self._MAPS["rivmsk"] if mask is not None: mask = self.staticmaps[mask].values - if mask is not None: + if snap_uparea and "uparea" in gdf_gauges.columns: + # Derive gauge map based on upstream area snapping + da, idxs, ids = workflows.gauge_map_uparea( + self.staticmaps, + gdf_gauges, + uparea_name="wflow_uparea", + wdw=wdw, + rel_error=rel_error, + logger=self.logger, + ) + else: # Derive gauge map da, idxs, ids = flw.gauge_map( self.staticmaps, @@ -1057,16 +986,7 @@ def create_gauges( ids_new = np.unique(da.values[da.values > 0]) idxs = idxs[np.isin(ids_old, ids_new)] ids = da.values.flat[idxs] - elif snap_uparea and "uparea" in gdf_gauges.columns: - # Derive gauge map based on upstream area snapping - da, idxs, ids = workflows.gauge_map_uparea( - self.staticmaps, - gdf_gauges, - uparea_name="wflow_uparea", - wdw=wdw, - rel_error=rel_error, - logger=self.logger, - ) + # Add to staticmaps mapname = f'{str(self._MAPS["gauges"])}_{basename}' self.set_staticmaps(da, name=mapname) diff --git a/tests/test_model_methods.py b/tests/test_model_methods.py index c964cfb5..523abf45 100644 --- a/tests/test_model_methods.py +++ b/tests/test_model_methods.py @@ -134,3 +134,28 @@ def test_setup_gauges(): ds_samp = mod.staticmaps[["wflow_river", "wflow_uparea"]].raster.sample(gdf, wdw=0) assert np.all(ds_samp["wflow_river"].values == 1) assert np.allclose(ds_samp["wflow_uparea"].values, gdf["uparea"].values, rtol=0.05) + + # Test with/without snapping + stations_fn = join(EXAMPLEDIR, "test_stations.csv") + mod.setup_gauges( + gauges_fn=stations_fn, + basename="stations_snapping", + snap_to_river=True, + mask=None, + ) + gdf_snap = mod.staticgeoms["gauges_stations_snapping"] + + mod.setup_gauges( + gauges_fn=stations_fn, + basename="stations_no_snapping", + snap_to_river=False, + mask=None, + ) + gdf_no_snap = mod.staticgeoms["gauges_stations_no_snapping"] + + # Check that not all geometries of gdf_snap and gdf_no_snap are equal + assert not gdf_snap.equals(gdf_no_snap) + # Find wich row is identical + equal = gdf_snap[gdf_snap["geometry"] == gdf_no_snap["geometry"]] + assert len(equal) == 1 + assert equal.index.values[0] == 1003 From 1c7bc1d3df07f73d1ee966d1c45c7da53654171f Mon Sep 17 00:00:00 2001 From: hboisgon Date: Wed, 7 Jun 2023 09:17:29 +0800 Subject: [PATCH 14/17] small bugfix for geodataframe --- hydromt_wflow/wflow.py | 11 ++++++++--- tests/test_model_methods.py | 3 ++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 1fcfe5a1..2e98a08d 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -3,6 +3,7 @@ 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 @@ -984,7 +985,7 @@ def setup_outlets( def setup_gauges( self, - gauges_fn: str, + gauges_fn: Union[str, Path, gpd.GeoDataFrame], index_col: Optional[str] = None, snap_to_river: Optional[bool] = True, mask: Optional[np.ndarray] = None, @@ -1021,7 +1022,7 @@ def setup_gauges( Parameters ---------- - gauges_fn : str, {"grdc"}, optional + gauges_fn : str, Path, gpd.GeoDataFrame, optional Known source name or path to gauges file geometry file, by default None. * Required variables if snap_uparea is True: ["uparea"] @@ -1058,7 +1059,11 @@ def setup_gauges( """ # Read data kwargs = {} - if isfile(gauges_fn): + if isinstance(gauges_fn, gpd.GeoDataFrame): + gdf_gauges = gauges_fn + if not np.all(np.isin(gdf_gauges.geometry.type, "Point")): + raise ValueError(f"{gauges_fn} contains other geometries than Point") + elif isfile(gauges_fn): # try to get epsg number directly, important when writting back data_catalog if hasattr(self.crs, "to_epsg"): code = self.crs.to_epsg() diff --git a/tests/test_model_methods.py b/tests/test_model_methods.py index 98682d53..217db625 100644 --- a/tests/test_model_methods.py +++ b/tests/test_model_methods.py @@ -119,7 +119,7 @@ def test_setup_gauges(): root = join(EXAMPLEDIR, "wflow_piave_subbasin") # Initialize model and read results - mod = WflowModel(root=root, mode="r", data_libs="artifact_data", logger=logger) + mod = WflowModel(root=root, mode="r", data_libs="artifact_data", logger=logger) # uparea rename not in the latest artifact_data version mod.data_catalog["grdc"].rename = {"area": "uparea"} mod.setup_gauges( @@ -161,6 +161,7 @@ def test_setup_gauges(): assert len(equal) == 1 assert equal.index.values[0] == 1003 + @pytest.mark.parametrize("elevtn_map", ["wflow_dem", "dem_subgrid"]) def test_setup_rivers(elevtn_map): # load netcdf file with floodplain layers From 0ebfdf2303847e8cc5cd065fdb1498c091970c66 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Wed, 7 Jun 2023 09:49:25 +0800 Subject: [PATCH 15/17] bugfix tests --- tests/data/wflow_piave_build_subbasin.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/data/wflow_piave_build_subbasin.ini b/tests/data/wflow_piave_build_subbasin.ini index 32f6b346..665484a7 100644 --- a/tests/data/wflow_piave_build_subbasin.ini +++ b/tests/data/wflow_piave_build_subbasin.ini @@ -1,4 +1,5 @@ [setup_config] +output.lateral.land.q = "q_land" output.lateral.land.h = "h_land" [setup_basemaps] From e93f8d9d1fd4fe2a1d5b87e336586e52d7a99f65 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Mon, 12 Jun 2023 15:43:54 +0800 Subject: [PATCH 16/17] update docstring --- hydromt_wflow/wflow.py | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 2e98a08d..20c7b5ca 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -903,6 +903,7 @@ def setup_config_output_timeseries( 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}", []) @@ -1008,11 +1009,27 @@ def setup_gauges( Supported gauge datasets include "grdc" or "" for user supplied csv or geometry files with gauge locations. If a csv file is provided, a "x" or "lon" and "y" or "lat" column is required - and the first column will be used as IDs in the map. If ``snap_to_river`` is set - to True, the gauge location will be snapped to the boolean river mask. If - ``derive_subcatch`` is set to True, an additional subcatch map is derived from + and the first column will be used as IDs in the map. + + There are three available methods to prepare the gauge map: + + * no snapping: ``mask=None``, ``snap_to_river=False``, ``snap_uparea=False``. + The gauge locations are used as is. + * snapping to mask: the gauge locations are snapped to a boolean mask map: either provide + ``mask`` or set ``snap_to_river=True`` to snap to the river (default). ``max_dist`` can be + used to set the maximum distance to snap to the mask. + * snapping based on upstream area matching: : ``snap_uparea=True``. The gauge locations + are snapped to the closest matching upstream area value. Requires gauges_fn to have + an `uparea` [m] column. The closest value will be looked for in a cell window of size ``wdw`` + and the difference between the gauge and the closest value should be smaller than ``rel_error``. + + If ``derive_subcatch`` is set to True, an additional subcatch map is derived from the gauge locations. + Finally the output locations can be added to wflow TOML file sections [csv] or [netcdf] + using the ``toml_output`` option. The ``gauge_toml_header`` and ``gauge_toml_param`` options + can be used to define the header and corresponding wflow variable names in the TOML file. + Adds model layers: * **wflow_gauges_source** map: gauge IDs map from source [-] (if gauges_fn) @@ -1035,11 +1052,11 @@ def setup_gauges( snap_uparea: bool, optional Snap gauges based on upstream area. Gauges_fn should have "uparea" in its attributes. max_dist : float, optional - Maximum distance between original and snapped point location. - A warning is logged if exceeded. By default 10 km. + Maximum distance [m] between original and snapped point location. + A warning is logged if exceeded. By default 10 000m. wdw: int, optional - Window size in number of cells around discharge boundary locations - to snap to, only used if ``snap_uparea`` is True. By default 3. + Window size in number of cells around the gauge locations + to snap uparea to, only used if ``snap_uparea`` is True. By default 3. rel_error: float, optional Maximum relative error (default 0.05) between the gauge location upstream area and the upstream area of From fc3c5692d30dc453b6fa2dbb198b0d1bf3e2d4e3 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Tue, 13 Jun 2023 08:49:33 +0800 Subject: [PATCH 17/17] update docstrings --- hydromt_wflow/wflow.py | 2 +- hydromt_wflow/workflows/gauges.py | 37 +++++++++++++++++++++++++++++-- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 20c7b5ca..72d93d5d 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -1020,7 +1020,7 @@ def setup_gauges( used to set the maximum distance to snap to the mask. * snapping based on upstream area matching: : ``snap_uparea=True``. The gauge locations are snapped to the closest matching upstream area value. Requires gauges_fn to have - an `uparea` [m] column. The closest value will be looked for in a cell window of size ``wdw`` + an `uparea` [km2] column. The closest value will be looked for in a cell window of size ``wdw`` and the difference between the gauge and the closest value should be smaller than ``rel_error``. If ``derive_subcatch`` is set to True, an additional subcatch map is derived from diff --git a/hydromt_wflow/workflows/gauges.py b/hydromt_wflow/workflows/gauges.py index a6fce47a..2b4ecf8a 100644 --- a/hydromt_wflow/workflows/gauges.py +++ b/hydromt_wflow/workflows/gauges.py @@ -21,9 +21,42 @@ def gauge_map_uparea( abs_error: float = 50, logger=logger, ): - """ """ + """ + Snaps point locations to grid cell with smallest difference in upstream area + within `wdw` around the original location if the local cell does not meet the + error criteria. + + Both the upstream area variable named ``uparea_name`` in ``ds`` and ``gdf`` as + well as ``abs_error`` should have the same unit (typically km2). + + Parameters + ---------- + ds : xr.Dataset + Dataset with upstream area variable. + gdf : gpd.GeoDataFrame + GeoDataFrame with gauge points and uparea column. + uparea_name : str, optional + Name of the upstream area variable in ``ds``, by default "wflow_uparea". + wdw : int, optional + Window size around the original location to search for the best matching cell, + by default 1. + rel_error : float, optional + Relative error threshold to accept the best matching cell, by default 0.05. + abs_error : float, optional + Absolute error threshold to accept the best matching cell, by default (50 km2). + + Returns + ------- + da : xr.DataArray + Gauge map with gauge points snapped to the best matching cell. + idxs_out : np.ndarray + Array of indices of the best matching cell. + ids_out : np.ndarray + Array of gauge point ids. + + """ if uparea_name not in ds: - return + raise ValueError(f"uparea_name {uparea_name} not found in ds.") ds_wdw = ds.raster.sample(gdf, wdw=wdw) logger.debug(