5. Update mask with water level and outflow boundary cells - including use o
NOTE: - As you can see now, also msk=2 values (red line) have been added along the coastal boundary - As you can see now, also msk=3 values (purple line) have been added along the lateral inland boundaries within the gdf_include shapefile - reset_bounds=True means that you start without initial boundary cells (of the specified type), if reset_bounds=False (default) you build on the existing boundary cells (if available)
# --> this is used in making subgrid tables
-
-# to ensure a correct lower roughness of the rivers in the domain, we 'burn in' the river with a lower manning roughness value.
-# the position of the river is based on global data, that is loaded below
+
# derive river from hydrography data based on a minimum river length (river_len)
+# and minimum upstream area (river_upa)
-# read river shapefile and add manning value to the attributes
-gdf=sf.data_catalog.get_geodataframe("rivers_lin2019_v1",geom=sf.region).to_crs(
- sf.crs
+sf.setup_river_inflow(
+ hydrography='merit_hydro',
+ river_len=1000,
+ river_upa=50,
+ keep_rivers_geom=True)
-gdf["geometry"]=gdf.buffer(50)
-gdf["manning"]=0.03
-# rasterize the manning value of gdf to the model grid
-da_manning=sf.grid.raster.rasterize(gdf,"manning",nodata=np.nan)
-
-# uncomment to plot either the raster or the vector data:
-# da_manning.plot(vmin=0, x='xc', y='yc', cmap='viridis')
-# gdf.plot()
+# Make a plot of model
+# note the src points and derived river network
+fig,ax=sf.plot_basemap(variable="dep",plot_bounds=False,bmap="sat",zoomlevel=12)
-
+
+
+
+
+
+
+
+
[10]:
-
# use the river manning raster in combination with vito land to derive the manning roughness file
+
# --> this is used when the making subgrid tables (step 8)
+
+# use the derived rivers to burn these into the model subgrid
+# note that the width and depth are arbitrary here
+
+gdf_riv=sf.geoms['rivers_inflow'].copy()
+gdf_riv['rivwth']=[100,50,50]# width [m]
+gdf_riv['rivdph']=1.5# depth [m]
+gdf_riv["manning"]=0.03# manning coefficient [s.m-1/3]
+gdf_riv[['geometry','rivwth','rivdph','manning']]
+
+# prepare the river dataset for the subgrid
+# instead of using the derived rivers, you could also use an aribitrary shapefile with centerlines
+# Other options include to add a river mask shapefile rather than a river width attribute
+# Note that the width is either specified on the river centerline or with a river mask
+# Also the river bed level (rivbed) can be specified instead of the river depth (rivdph).
+
+datasets_riv=[{'centerlines':gdf_riv}]
+
To ensure a correct lower roughness of the rivers in the domain, we ‘burn in’ the river with a lower manning roughness value. There is two different ways to do this:
+
Rasterize the manning value of gdf to the model grid and use this as a manning raster in datasets_rgh
+
Provide a manning roughness for datasets_riv while burning in rivers into the subgrid
+
+
[11]:
+
+
+
# 1. rasterize the manning value of gdf to the model grid and use this as a manning raster
+# gdf_riv_buf = gdf_riv.assign(geometry=gdf_riv.buffer(gdf_riv['rivwth']/2))
+# da_manning = sf.grid.raster.rasterize(gdf_riv_buf, "manning", nodata=np.nan)
+
+# use the river manning raster in combination with vito land to derive the manning roughness file# NOTE that we can combine in-memory data with data from the data catalog
-datasets_rgh=[{"manning":da_manning},{"lulc":"vito"}]
+# datasets_rgh = [{"manning": da_manning}, {"lulc": "vito"}]
-sf.setup_manning_roughness(
- datasets_rgh=datasets_rgh,
- manning_land=0.04,
- manning_sea=0.02,
- rgh_lev_land=0,# the minimum elevation of the land
-)
-_=sf.plot_basemap(variable="manning",plot_bounds=False,bmap="sat",zoomlevel=12)
+# uncomment to plot either the raster or the vector data:
+# da_manning.plot(vmin=0, x='xc', y='yc', cmap='viridis')
+# gdf.plot()
-
-
+
+
[12]:
+
-
-
-The nodata value None is not in the reclass table.None will be used for the params.
-No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
-
+
# 2. provide the manning value to the datasets_riv and burn in the river with a lower manning value
+# a global manning value can be provided to the datasets_riv, or a manning value per river segment can be provided
+# NOTE when using this method, we don't need to provide the river manning values to the datasets_rgh
+datasets_rgh=[{"lulc":"vito"}]
+
-
-
-
-
+
+
[13]:
+
+
+
# NOTE This step can be skipped because we prepare a subgrid model,
+# the manning roughness is than included in the subgrid tables
+# For models without subgrid uncomment the following lines.
+
+# sf.setup_manning_roughness(
+# datasets_rgh=datasets_rgh,
+# manning_land=0.04,
+# manning_sea=0.02,
+# rgh_lev_land=0, # the minimum elevation of the land
+# )
+# _ = sf.plot_basemap(variable="manning", plot_bounds=False, bmap="sat", zoomlevel=12)
+
Subgrid derived tables are used to better capture the elevation and roughness of your domain, to either improve your results, or to allow you to run on a courser grid resolution (means faster simulation). For more info about subgrid tables, click here.
You as user can specify multiple settings about how the subgrid derived tables should be made.
Every single grid cell of the flux grid of the size inp.dx by inp.dy is defined into subgrid pixels (default nr_subgrid_pixels = 20). For every subgrid pixel the topobathy data is loaded, ideally this consists of high-resolution DEM datasets that you specify as user.
In this example with dx=dy=50m, having nr_subgrid_pixels = 20 means we are loading data onto a 2.5 m subpixel grid However, the input data of Gebco and Merit_hydro is way coarser, therefore let’s set the ratio to 5 for now.
-
[11]:
+
[14]:
# Every single grid cell of the flux grid of the size inp.dx by inp.dy is defined into subgrid pixels (default is 20, nr_subgrid_pixels = 20).
@@ -758,6 +812,7 @@
11. Add an upstream discharge time-series as forcing:#
a. specify discharge input locations: srcfile
For more info about what the srcfile is, click here
b. specify discharge time-series: disfile
For more info about what the disfile is, click here
-
[18]:
+
[21]:
-
# We follow exactly the same steps as for the water level forcing, but now specify 1 location where we specify discharges in m3/s
-x=[321732]
-y=[5047136]
+
# We now use the previously created src discharge points (step 6)
+# Alternativly you can also create a new geodataframe with points, similar to the water level forcing points
-# add to Geopandas dataframe as needed by HydroMT
-pnts=gpd.points_from_xy(x,y)
-index=[1]# NOTE that the index should start at one
-src=gpd.GeoDataFrame(index=index,geometry=pnts,crs=sf.crs)
-
-time=pd.date_range(
- start=utils.parse_datetime(sf.config["tstart"]),
- end=utils.parse_datetime(sf.config["tstop"]),
- periods=3,
+# make up some discharge data
+index=sf.forcing['dis'].index
+dis=np.array(
+ [[2.0,1.0],
+ [5.0,2.0],
+ [2.0,1.0]])
-
-# and the actual water levels, in this case for input location 1 a water level rising from 0 to 2 meters and back to 0:
-dis=[[2.0],[5.0],[2.0]]
-
dispd=pd.DataFrame(index=time,columns=index,data=dis)
-# now we call the function setup_discharge_forcing
-sf.setup_discharge_forcing(timeseries=dispd,locations=src)
+# now we call the function setup_discharge_forcing, which adds the discharge forcing to the src points
+sf.setup_discharge_forcing(timeseries=dispd)
-# NOTE: the discharge forcing data is now stored in the sf.forcing dictionary
+# # NOTE: the discharge forcing data is now stored in the sf.forcing dictionarysf.forcing.keys()
In case you want to add other types of forcing, see the example notebook example_forcing.ipynb for other types. Or read more about this in the SFINCS manual
# hourly rainfall rates of ECMWF' ERA5 data for the specific area and period have been made available for this period in the artefact data
@@ -1083,20 +1121,20 @@
In SFINCS, a weirfile is often used to explicity account for line-element features such as dikes, dunes or floodwalls. Read more about structures in the SFINCS manual
-
[21]:
+
[24]:
# In this example specify a 'line' style shapefile for the location of the weir to be added
@@ -1129,20 +1167,20 @@
Your basemap and forcing figures are saved in the folder ‘figs’, GIS files (tiff & geojson) of your model setup in ‘gis’ and merged elevation and manning roughness on subgrid resolution in ‘subgrid’.
List of dictionaries with topobathy data.
Each should minimally contain a data catalog source name, data file path,
-or xarray raster object (‘elevtn’)
-
Optional merge arguments include ‘zmin’, ‘zmax’, ‘mask’, ‘offset’,
-‘reproj_method’, and ‘merge_method’. e.g.:
-[
List of dictionaries with Manning’s n datasets. Each dictionary should at
+least contain one of the following:
* (1) manning: filename (or Path) of gridded data with manning values
-* (2) lulc (and reclass_table) :a combination of a filename of gridded landuse/landcover and a mapping table.
-In additon, optional merge arguments can be provided e.g.: merge_method, gdf_valid_fn
-
buffer_cells (int, optional) – Number of cells between datasets to ensure smooth transition of bed levels, by default 0
+* (2) lulc (and reclass_table) :a combination of a filename of gridded
+landuse/landcover and a mapping table. In additon, optional merge arguments
+can be provided, e.g.: [
+
List of dictionaries with river datasets. Each dictionary should at least
contain a river centerline data and optionally a river mask:
+* centerlines: filename (or Path) of river centerline with attributes
+
+
rivwth (river width [m]; required if not river mask provided),
+rivdph or rivbed (river depth [m]; river bedlevel [m+REF]),
+manning (Manning’s n [s/m^(1/3)]; optional)
+
+
+
mask (optional): filename (or Path) of river mask
+
+
river attributes (optional): “rivdph”, “rivbed”, “rivwth”, “manning”
to fill missing values
+
+
+
+
+
arguments to the river burn method (optional):
segment_length [m] (default 500m) and riv_bank_q [0-1] (default 0.5)
+which used to estimate the river bank height in case river depth is provided.
+
+
+
For more info see :py:function:~hydromt.workflows.bathymetry.burn_river_rect
buffer_cells (int, optional) – Number of cells between datasets to ensure smooth transition of bed levels,
+by default 0
nbins (int, optional) – Number of bins in which hypsometry is subdivided, by default 10
nr_subgrid_pixels (int, optional) – Number of subgrid pixels per computational cell, by default 20
nrmax (int, optional) – Maximum number of cells per subgrid-block, by default 2000
These blocks are used to prevent memory issues while working with large datasets
-
max_gradient (float, optional) – If slope in hypsometry exceeds this value, then smoothing is applied, to prevent numerical stability problems, by default 5.0
+
max_gradient (float, optional) – If slope in hypsometry exceeds this value, then smoothing is applied,
+to prevent numerical stability problems, by default 5.0
z_minimum (float, optional) – Minimum depth in the subgrid tables, by default -99999.0
manning_land (float, optional) – Constant manning roughness values for land and sea, by default 0.04 and 0.02 s.m-1/3
-Note that these values are only used when no Manning’s n datasets are provided, or to fill the nodata values
+Note that these values are only used when no Manning’s n datasets are provided,
+or to fill the nodata values
manning_sea (float, optional) – Constant manning roughness values for land and sea, by default 0.04 and 0.02 s.m-1/3
-Note that these values are only used when no Manning’s n datasets are provided, or to fill the nodata values
-
rgh_lev_land (float, optional) – Elevation level to distinguish land and sea roughness (when using manning_land and manning_sea), by default 0.0
-
write_dep_tif (bool, optional) – Create geotiff of the merged topobathy on the subgrid resolution, by default False
-
write_man_tif (bool, optional) – Create geotiff of the merged roughness on the subgrid resolution, by default False
+Note that these values are only used when no Manning’s n datasets are provided,
+or to fill the nodata values
+
rgh_lev_land (float, optional) – Elevation level to distinguish land and sea roughness
+(when using manning_land and manning_sea), by default 0.0
+
write_dep_tif (bool, optional) – Write geotiff of the merged topobathy / roughness on the subgrid resolution.
+These files are not used by SFINCS, but can be used for visualisation and
+downscaling of the floodmaps. Unlinke the SFINCS files it is written
+to disk at execution of this method. By default False
+
write_man_tif (bool, optional) – Write geotiff of the merged topobathy / roughness on the subgrid resolution.
+These files are not used by SFINCS, but can be used for visualisation and
+downscaling of the floodmaps. Unlinke the SFINCS files it is written
+to disk at execution of this method. By default False
Return river bank height estimated as from height above nearest drainage
-(HAND) values adjecent to river cells. For each feature in gdf_riv the nearest
-river bank cells are identified and the bank heigth is estimated based on a quantile
-value q.
nmin (int, optional) – Minimum threshold for valid river bank cells, by default 20
-
q (float, optional) – quantile [0-100] for river bank estimate, by default 25.0
+
da_elv (xr.DataArray) – DEM and manning raster to burn river depth and manning values into
+
da_man (xr.DataArray) – DEM and manning raster to burn river depth and manning values into
+
gdf_riv (gpd.GeoDataFrame) – River center lines.
+
gdf_zb (gpd.GeoDataFrame, optional) – Point locations with a ‘rivbed’ river bed level [m+REF] column, by defualt None
+
gdf_riv_mask (gpd.GeoDataFrame, optional) – Mask in which to interpolate z values, by default None.
+If provided, ‘rivwth’ column is not required in gdf_riv.
+
segment_length (float, optional) – Approximate river segment length [m], by default 500
+
riv_bank_q (float, optional) – quantile [0-1] for river bank estimation, by default 0.25
+
rivwth_name (str, optional) – river width [m], depth [m], bed level [m+REF], & manning [s.m-1/3] column names
+in gdf_riv, by default “rivwth”, “rivdph”, “rivbed”, and “manning”
+
rivdph_name (str, optional) – river width [m], depth [m], bed level [m+REF], & manning [s.m-1/3] column names
+in gdf_riv, by default “rivwth”, “rivdph”, “rivbed”, and “manning”
+
rivbed_name (str, optional) – river width [m], depth [m], bed level [m+REF], & manning [s.m-1/3] column names
+in gdf_riv, by default “rivwth”, “rivdph”, “rivbed”, and “manning”
+
manning_name (str, optional) – river width [m], depth [m], bed level [m+REF], & manning [s.m-1/3] column names
+in gdf_riv, by default “rivwth”, “rivdph”, “rivbed”, and “manning”
-
Returns:
-
-
rivbank_dz (np.ndarray) – riverbank elevations for each segment in gdf_riv
-
da_riv_mask, da_bnk_mask (xr.DataArray:) – River and river-bank masks
Estimate river bedlevel zb using gradually varying flow (gvf), manning’s equation
-(manning) or a power-law relation (powlaw) rivdph_method. The river is based on flow
-directions with and minimum upstream area threshold.
-
-
Parameters:
-
-
ds (xr.Dataset) – Model map layers containing elevnt_name, uparea_name and rivmsk_name (optional)
-variables.
-
flwdir (pyflwdir.FlwdirRaster) – Flow direction object
-
gdf_riv (gpd.GeoDataFrame, optional) – River attribute data with “qbankfull” and “rivwth” data, by default None
-
gdf_qbf (gpd.GeoDataFrame, optional) – Bankfull river discharge data with “qbankfull” column, by default None
-
rivdph_method ({'gvf', 'manning', 'powlaw', 'geom'}) – River depth estimate method, by default ‘gvf’
-
rivwth_method ({'geom', 'mask'}) – River width estimate method, by default ‘mask’
-
river_upa (float, optional) – Minimum upstream area threshold for rivers [km2], by default 100.0
-
river_len (float, optional) – Mimimum river length [m] within the model domain to define river cells, by default 1000.
-
min_rivwth (float, optional) – Minimum river width [m] (by default 50.0 m) and depth [m] (by default 1.0 m)
-
min_rivdph (float, optional) – Minimum river width [m] (by default 50.0 m) and depth [m] (by default 1.0 m)
-
segment_length (float, optional) – Approximate river segment length [m], by default 5e3
min_convergence (float, optional) – Minimum width convergence threshold to define estuaries [m/m], by default 0.01
-
max_dist (float, optional) – Maximum distance threshold to spatially merge gdf_riv and gdf_qbf, by default 100.0
-
rivbank (bool, optional) – If True (default), approximate the reference elevation for the river depth based
-on the river bankfull elevation at cells neighboring river cells. Otherwise
-use the elevation of the local river cell as reference level.
-
rivbankq (float, optional) – quantile [1-100] for river bank estimation, by default 25
-
constrain_estuary (bool, optional) – If True (default) fix the river depth in estuaries based on the upstream river depth.
-
constrain_rivbed (bool, optional) – If True (default) correct the river bed level to be hydrologically correct
-
-
-
Returns:
-
-
gdf_riv (gpd.GeoDataFrame) – River segments with bed level (zb) estimates
self.write_forcing()self.write_states()# config last; might be udpated when writing maps, states or forcing
- self.write_config()
+ self.write_config()# write data catalog with used data sources
- # self.write_data_catalog() # new in hydromt v0.4.4
+ self.write_data_catalog()# new in hydromt v0.4.4
[docs]defread_grid(self,data_vars:Union[List,str]=None)->None:"""Read SFINCS binary grid files and save to `grid` attribute.
@@ -3834,6 +3863,79 @@
Source code for hydromt_sfincs.sfincs
self.logger.warning(f"Unknown key {key} in datasets_rgh. Ignoring.")datasets_out.append(dd)
+ returndatasets_out
+
+ def_parse_datasets_riv(self,datasets_riv):
+"""Parse filenames or paths of Datasets in list of dictionaries
+ datasets_riv into xr.DataArrays and gdf.GeoDataFrames:
+
+ see SfincsModel.setup_subgrid for details
+ """
+ # option 1: rectangular river cross-sections based on river centerline
+ # depth/bedlevel, manning attributes are specified on the river centerline
+ # TODO: make this work with LineStringZ geometries for bedlevel
+ # the width is either specified on the river centerline or river mask
+ # option 2: (TODO): irregular river cross-sections
+ # cross-sections are specified as a series of points (river_crosssections)
+ parse_keys=[
+ "centerlines",
+ "mask",
+ "gdf_riv",
+ "gdf_riv_mask",
+ ]
+ copy_keys=[]
+ attrs=["rivwth","rivdph","rivbed","manning"]
+
+ datasets_out=[]
+ fordatasetindatasets_riv:
+ dd={}
+
+ # parse rivers
+ if"centerlines"indataset:
+ rivers=dataset.get("centerlines")
+ ifisinstance(rivers,str)andriversinself.geoms:
+ gdf_riv=self.geoms[rivers].copy()
+ else:
+ gdf_riv=self.data_catalog.get_geodataframe(
+ rivers,
+ geom=self.mask.raster.box,
+ buffer=1e3,# 1km
+ ).to_crs(self.crs)
+ # update missing attributes based on global values
+ forkeyinattrs:
+ ifkeyindataset:
+ value=dataset.pop(key)
+ ifkeynotingdf_riv.columns:# update all
+ gdf_riv[key]=value
+ elifnp.any(np.isnan(gdf_riv[key])):# fill na
+ gdf_riv[key]=gdf_riv[key].fillna(value)
+ ifnotgdf_riv.columns.isin(["rivbed","rivdph"]).any():
+ raiseValueError("No 'rivbed' or 'rivdph' attribute found.")
+ else:
+ raiseValueError("No 'centerlines' dataset provided.")
+ dd.update({"gdf_riv":gdf_riv})
+
+ # parse mask
+ if"mask"indataset:
+ gdf_riv_mask=self.data_catalog.get_geodataframe(
+ dataset.get("mask"),
+ geom=self.mask.raster.box,
+ )
+ dd.update({"gdf_riv_mask":gdf_riv_mask})
+ elif"rivwth"notingdf_riv:
+ raiseValueError(
+ "Either mask must be provided or centerlines "
+ "should contain a 'rivwth' attribute."
+ )
+
+ # copy remaining keys
+ forkey,valueindataset.items():
+ ifkeyincopy_keysandkeynotindd:
+ dd.update({key:value})
+ elifkeynotincopy_keys+parse_keys:
+ self.logger.warning(f"Unknown key {key} in datasets_riv. Ignoring.")
+ datasets_out.append(dd)
+
returndatasets_out
Source code for hydromt_sfincs.workflows.bathymetry
"""Workflows, to estimate river bathymetry and burn these in a DEM."""importlogging
-fromtypingimportTupleimportgeopandasasgpdimportnumpyasnp
-importpyflwdirimportxarrayasxr
-fromhydromt.gis_utilsimportnearest,nearest_merge,spread2d
-fromhydromt.workflowsimportrivers
+fromhydromt.gis_utilsimportnearest,parse_crsfromscipyimportndimage
+fromscipy.interpolateimportinterp1d
+fromshapely.geometryimportLineString,MultiLineString,MultiPoint,Point
+fromshapely.opsimportlinemerge,snap,split,unary_unionlogger=logging.getLogger(__name__)__all__=[
- "get_rivbank_dz",
- "get_river_bathymetry",
- "burn_river_zb",
+ "burn_river_rect",]
-
[docs]defget_rivbank_dz(
- gdf_riv:gpd.GeoDataFrame,
- da_msk:xr.DataArray,
- da_hnd:xr.DataArray,
- nmin:int=20,
- q:float=25.0,
-)->np.ndarray:
-"""Return river bank height estimated as from height above nearest drainage
- (HAND) values adjecent to river cells. For each feature in `gdf_riv` the nearest
- river bank cells are identified and the bank heigth is estimated based on a quantile
- value `q`.
+def_split_line_by_point(
+ line:LineString,point:Point,tolerance:float=1.0e-5
+)->MultiLineString:
+"""Split a single line with a point. Parameters ----------
- gdf_riv : gpd.GeoDataFrame
- River segments
- da_msk : xr.DataArray of bool
- River mask
- da_hnd : xr.DataArray of float
- Height above nearest drain (HAND) map
- nmin : int, optional
- Minimum threshold for valid river bank cells, by default 20
- q : float, optional
- quantile [0-100] for river bank estimate, by default 25.0
+ line : LineString
+ line
+ point : Point
+ point
+ tolerance : float, optional
+ tolerance to snap the point to the line, by default 1.0e-5 Returns -------
- rivbank_dz: np.ndarray
- riverbank elevations for each segment in `gdf_riv`
- da_riv_mask, da_bnk_mask: xr.DataArray:
- River and river-bank masks
+ MultiLineString
+ splitted line """
- # rasterize streams
- gdf_riv=gdf_riv.copy()
- gdf_riv["segid"]=np.arange(1,gdf_riv.index.size+1,dtype=np.int32)
- valid=gdf_riv.length>0# drop pits or any zero length segments
- segid=da_hnd.raster.rasterize(gdf_riv[valid],"segid").astype(np.int32)
- segid.raster.set_nodata(0)
- segid.name="segid"
- # NOTE: the assumption is that banks are found in cells adjacent to any da_msk cell
- da_msk=da_msk.raster.reproject_like(da_hnd,method="nearest")
- _mask=ndimage.binary_fill_holes(da_msk)# remove islands
- mask=ndimage.binary_dilation(_mask,np.ones((3,3)))
- da_mask=xr.DataArray(
- coords=da_hnd.raster.coords,dims=da_hnd.raster.dims,data=mask
- )
- da_mask.raster.set_crs(da_hnd.raster.crs)
- # find nearest stream segment for all river bank cells
- segid_spread=spread2d(da_obs=segid,da_mask=da_mask)
- # get edge of riv mask -> riv banks
- da_bnk_mask=np.logical_and(da_hnd>0,np.logical_xor(da_mask,_mask))
- da_riv_mask=np.logical_and(
- np.logical_and(da_hnd>=0,da_msk),np.logical_xor(da_bnk_mask,da_mask)
- )
- # get median HAND for each stream -> riv bank dz
- rivbank_dz=ndimage.labeled_comprehension(
- da_hnd.values,
- labels=np.where(da_bnk_mask,segid_spread["segid"].values,np.int32(0)),
- index=gdf_riv["segid"].values,
- func=lambdax:0ifx.size<nminelsenp.percentile(x,q),
- out_dtype=da_hnd.dtype,
- default=-9999,
- )
- returnrivbank_dz,da_riv_mask,da_bnk_mask
-
-
-
[docs]defget_river_bathymetry(
- ds:xr.Dataset,
- flwdir:pyflwdir.FlwdirRaster,
- gdf_riv:gpd.GeoDataFrame=None,
- gdf_qbf:gpd.GeoDataFrame=None,
- rivdph_method:str="gvf",
- rivwth_method:str="mask",
- river_upa:float=100.0,
- river_len:float=1e3,
- min_rivdph:float=1.0,
- min_rivwth:float=50.0,
- segment_length:float=5e3,
- smooth_length:float=10e3,
- min_convergence:float=0.01,
- max_dist:float=100.0,
- rivbank:bool=True,
- rivbankq:float=25,
- constrain_estuary:bool=True,
- constrain_rivbed:bool=True,
- elevtn_name:str="elevtn",
- uparea_name:str="uparea",
- rivmsk_name:str="rivmsk",
- logger=logger,
- **kwargs,
-)->Tuple[gpd.GeoDataFrame,xr.DataArray]:
-"""Estimate river bedlevel zb using gradually varying flow (gvf), manning's equation
- (manning) or a power-law relation (powlaw) rivdph_method. The river is based on flow
- directions with and minimum upstream area threshold.
+ returnsplit(snap(line,point,tolerance),point)
+
+
+def_line_to_points(line:LineString,dist:float=None,n:int=None)->MultiPoint:
+"""Get points along line based on a distance `dist` or number of points `n`. Parameters ----------
- ds : xr.Dataset
- Model map layers containing `elevnt_name`, `uparea_name` and `rivmsk_name` (optional)
- variables.
- flwdir : pyflwdir.FlwdirRaster
- Flow direction object
- gdf_riv : gpd.GeoDataFrame, optional
- River attribute data with "qbankfull" and "rivwth" data, by default None
- gdf_qbf : gpd.GeoDataFrame, optional
- Bankfull river discharge data with "qbankfull" column, by default None
- rivdph_method : {'gvf', 'manning', 'powlaw', 'geom'}
- River depth estimate method, by default 'gvf'
- rivwth_method : {'geom', 'mask'}
- River width estimate method, by default 'mask'
- river_upa : float, optional
- Minimum upstream area threshold for rivers [km2], by default 100.0
- river_len: float, optional
- Mimimum river length [m] within the model domain to define river cells, by default 1000.
- min_rivwth, min_rivdph: float, optional
- Minimum river width [m] (by default 50.0 m) and depth [m] (by default 1.0 m)
- segment_length : float, optional
- Approximate river segment length [m], by default 5e3
- smooth_length : float, optional
- Approximate smoothing length [m], by default 10e3
- min_convergence : float, optional
- Minimum width convergence threshold to define estuaries [m/m], by default 0.01
- max_dist : float, optional
- Maximum distance threshold to spatially merge `gdf_riv` and `gdf_qbf`, by default 100.0
- rivbank: bool, optional
- If True (default), approximate the reference elevation for the river depth based
- on the river bankfull elevation at cells neighboring river cells. Otherwise
- use the elevation of the local river cell as reference level.
- rivbankq : float, optional
- quantile [1-100] for river bank estimation, by default 25
- constrain_estuary : bool, optional
- If True (default) fix the river depth in estuaries based on the upstream river depth.
- constrain_rivbed : bool, optional
- If True (default) correct the river bed level to be hydrologically correct
+ line : LineString
+ line
+ dist : float
+ distance between points
+ n: integer
+ numer of points Returns -------
- gdf_riv: gpd.GeoDataFrame
- River segments with bed level (zb) estimates
- da_msk: xr.DataArray:
- River mask
+ MultiPoint
+ points """
- raster_kwargs=dict(coords=ds.raster.coords,dims=ds.raster.dims)
- da_elv=ds[elevtn_name]
-
- # get vector of stream segments
- da_upa=ds[uparea_name]
- rivd8=da_upa>river_upa
- ifriver_len>0:
- dx_headwater=np.where(flwdir.upstream_sum(rivd8)==0,flwdir.distnc,0)
- rivlen=flwdir.fillnodata(dx_headwater,nodata=0,direction="down",how="max")
- rivd8=np.logical_and(rivd8,rivlen>=river_len)
- feats=flwdir.streams(
- max_len=int(round(segment_length/ds.raster.res[0])),
- uparea=da_upa.values,
- elevtn=da_elv.values,
- rivdst=flwdir.distnc,
- strord=flwdir.stream_order(mask=rivd8),
- mask=rivd8,
- )
- gdf_stream=gpd.GeoDataFrame.from_features(feats,crs=ds.raster.crs)
- # gdf_stream = gdf_stream[gdf_stream['rivdst']>0] # remove pits
- flw=pyflwdir.from_dataframe(gdf_stream.set_index("idx"))
- _=flw.main_upstream(uparea=gdf_stream["uparea"].values)
-
- # merge gdf_riv with gdf_stream
- ifgdf_rivisnotNone:
- cols=[cforcin["rivwth","qbankfull"]ifcingdf_riv]
- gdf_riv=nearest_merge(gdf_stream,gdf_riv,columns=cols,max_dist=max_dist)
- gdf_riv["rivlen"]=gdf_riv["rivdst"]-flw.downstream(gdf_riv["rivdst"])
- else:
- gdf_riv=gdf_stream
- # merge gdf_qbf (qbankfull) with gdf_riv
- ifgdf_qbfisnotNoneand"qbankfull"ingdf_qbf.columns:
- if"qbankfull"ingdf_riv:
- gdf_riv=gdf_riv.drop(colums="qbankfull")
- idx_nn,dists=nearest(gdf_qbf,gdf_riv)
- valid=dists<max_dist
- gdf_riv.loc[idx_nn[valid],"qbankfull"]=gdf_qbf["qbankfull"].values[valid]
- logger.info(f"{sum(valid)}/{len(idx_nn)} qbankfull boundary points set.")
- check_q=rivdph_method=="geom"
- assertcheck_qor"qbankfull"ingdf_riv.columns,'gdf_riv has no "qbankfull" data'
- # propagate qbankfull and rivwth values
- ifrivwth_method=="geom"and"rivwth"notingdf_riv.columns:
- logger.info(f'Setting missing "rivwth" to min_rivwth {min_rivwth:.2f}m')
- gdf_riv["rivwth"]=min_rivwth
- ifrivdph_method=="geom"and"rivdph"notingdf_riv.columns:
- logger.info(f'Setting missing "rivdph" to min_rivdph {min_rivdph:.2f}m')
- gdf_riv["rivdph"]=min_rivdph
- forcolin["qbankfull","rivwth","rivdph"]:
- ifcolnotingdf_riv.columns:
- continue
- data=gdf_riv[col].fillna(-9999)
- data=flw.fillnodata(data,-9999,direction="down",how="max")
- gdf_riv[col]=np.maximum(0,data)
-
- # create river mask with river polygon
- ifrivmsk_namenotindsand"rivwth":
- assertgdf_riv.crs.is_projected
- gdf_riv_buf=gdf_riv.copy()
- buf=np.maximum(gdf_riv_buf["rivwth"]/2,1)
- gdf_riv_buf["geometry"]=gdf_riv_buf.buffer(buf)
- da_msk=np.logical_and(
- ds.raster.geometry_mask(gdf_riv_buf),da_elv!=da_elv.raster.nodata
- )
- elifrivmsk_nameinds:
- # merge river mask with river line
- da_msk=ds.raster.geometry_mask(gdf_riv,all_touched=True)
- da_msk=np.logical_or(da_msk,ds[rivmsk_name]>0)
+ ifdistisnotNone:
+ distances=np.arange(0,line.length,dist)
+ elifnisnotNone:
+ distances=np.linspace(0,line.length,n)else:
- raiseValueError("No river width or river mask provided.")
-
- ## get zs
- da_hnd=xr.DataArray(flwdir.hand(rivd8.values,da_elv.values),**raster_kwargs)
- da_hnd=da_hnd.where(da_elv>=0,0)
- da_hnd.raster.set_crs(ds.raster.crs)
- ifrivbank:
- logger.info("Deriving bankfull river surface elevation from its banks.")
- dz=get_rivbank_dz(gdf_riv,da_msk=da_msk,da_hnd=da_hnd,q=rivbankq)[0]
- dz=flw.fillnodata(dz,-9999,direction="down",how="max")
+ ValueError('Either "dist" or "n" should be provided')
+ points=unary_union(line.interpolate(distances))
+ returnpoints
+
+
+def_split_line_equal(
+ line:LineString,approx_length:float,tolerance:float=1.0e-5
+)->MultiLineString:
+"""Split line into segments with equal length.
+
+ Parameters
+ ----------
+ line : LineString
+ line to split
+ approx_length : float
+ Based in this approximate length the number of line segments is determined.
+ tolerance : float, optional
+ tolerance to snap the point to the line, by default 1.0e-5
+
+ Returns
+ -------
+ MultiLineString
+ line splitted in segments of equal length
+ """
+ n=int(np.floor(line.length/approx_length))
+ ifn<=1:
+ returnlineelse:
- logger.info("Deriving bankfull river surface elevation at its center line.")
- dz=0
- gdf_riv["zs"]=np.maximum(
- gdf_riv["elevtn"],flw.dem_adjust(gdf_riv["elevtn"]+dz)
- )
- gdf_riv["rivbank_dz"]=gdf_riv["zs"]-gdf_riv["elevtn"]
-
- # estimate stream segment average width from river mask
- ifrivwth_method=="mask":
- logger.info("Deriving river segment average width from permanent water mask.")
- rivwth=rivers.river_width(gdf_riv,da_rivmask=da_msk)
- gdf_riv["rivwth"]=flw.fillnodata(rivwth,-9999,direction="down",how="max")
- smooth_n=int(np.round(smooth_length/segment_length/2))
- ifsmooth_n>0:
- logger.info(f"Smoothing river width (n={smooth_n}).")
- gdf_riv["rivwth"]=flw.moving_average(
- gdf_riv["rivwth"],n=smooth_n,restrict_strord=True
- )
- gdf_riv["rivwth"]=np.maximum(gdf_riv["rivwth"],min_rivwth)
-
- # estimate river depth, smooth and correct
- if"rivdph"notingdf_riv.columns:
- gdf_riv["rivdph"]=rivers.river_depth(
- data=gdf_riv,
- flwdir=flw,
- method=rivdph_method,
- rivzs_name="zs",
- min_rivdph=min_rivdph,
- **kwargs,
- )
- ifconstrain_estuary:
- # set width from mask and depth constant in estuaries
- # estuaries based on convergence of width from river mask
- gdf_riv["estuary"]=flw.classify_estuaries(
- elevtn=gdf_riv["elevtn"],
- rivwth=gdf_riv["rivwth"],
- rivdst=gdf_riv["rivdst"],
- min_convergence=min_convergence,
- )
- gdf_riv["rivdph"]=np.where(
- gdf_riv["estuary"]==1,-9999,gdf_riv["rivdph"].values
- )
- gdf_riv["rivdph"]=flw.fillnodata(gdf_riv["rivdph"],-9999,"down")
- gdf_riv["rivdph"]=np.maximum(gdf_riv["rivdph"],min_rivdph)
-
- # calculate bed level from river depth
- gdf_riv["zb"]=gdf_riv["zs"]-gdf_riv["rivdph"]
- ifconstrain_rivbed:
- gdf_riv["zb"]=flw.dem_adjust(gdf_riv["zb"])
- gdf_riv["zb"]=np.minimum(gdf_riv["zb"],gdf_riv["elevtn"])
- gdf_riv["rivdph"]=gdf_riv["zs"]-gdf_riv["zb"]
-
- # calculate rivslp
- if"rivslp"notingdf_riv:
- dz=gdf_riv["zb"]-flw.downstream(gdf_riv["zb"])
- dx=gdf_riv["rivdst"]-flw.downstream(gdf_riv["rivdst"])
- # fill nodata with upstream neighbors and set lower bound of zero
- gdf_riv["rivslp"]=np.maximum(
- 0,flw.fillnodata(np.where(dx>0,dz/dx,-1),-1)
- )
+ split_points=_line_to_points(line,n=n)
+ return_split_line_by_point(line,split_points,tolerance=tolerance)
- returngdf_riv,da_msk
+defsplit_line_equal(gdf:gpd.GeoDataFrame,dist:float)->gpd.GeoDataFrame:
+"""Split lines in `gdf` into segments with equal length `dist`.
-
[docs]defburn_river_zb(
- gdf_riv:gpd.GeoDataFrame,
- da_elv:xr.DataArray,
- da_msk:xr.DataArray,
- flwdir:pyflwdir.FlwdirRaster=None,
- river_d4:bool=True,
+ Parameters
+ ----------
+ gdf : gpd.GeoDataFrame
+ lines
+ dist : float
+ approximate length of splitted line segments
+
+ Returns
+ -------
+ gpd.GeoDataFrame
+ splitted lines
+ """
+
+ def_split_geom(x:gpd.GeoSeries,dist:float=dist)->MultiLineString:
+ return_split_line_equal(x.geometry,dist)
+
+ gdf_splitted=gdf.assign(geometry=gdf.apply(_split_geom,axis=1)).explode()
+ returngdf_splitted
+
+
+definterp_along_line_to_grid(
+ da_mask:xr.DataArray,
+ gdf_lines:gpd.GeoDataFrame,
+ gdf_zb:gpd.GeoDataFrame,
+ column_names:list=["z"],logger=logger,
- **kwargs,
-):
-"""Burn bedlevels from `gdf_riv` (column zb) into the DEM `da_elv` at river cells
- indicated in `da_msk`.
+)->xr.DataArray:
+"""Interpolate values from `gdf_zb` along
+ lines from `gdf_lines` to grid cells from `da_mask`. Parameters ----------
- gdf_riv: gpd.GeoDataFrame
- River segments with bed level (zb) estimates
- da_elv : xr.DataArray of float
- Elevation raster
- da_msk: xr.DataArray of bool:
- River mask
- flwdir : pyflwdir.FlwdirRaster, optional
- Flow direction object
- river_d4 : bool, optional
- If True (default) ensure river cells have D4 connectivity
+ da_mask : xr.DataArray, optional
+ Boolean mask. Only cells with True values are interpolated.
+ gdf_lines : gpd.GeoDataFrame
+ center lines
+ gdf_zb : gpd.GeoDataFrame
+ point locations with alues to interpolate
+ column_names : list, optional
+ column names to interpolate, by default ["z"] Returns -------
- da_elv1: xr.DataArray
- DEM with bedlevels burned in.
+ xr.Dataset
+ interpolated values
+ """
+ ifnotall([cingdf_zb.columnsforcincolumn_names]):
+ missing=[cforcincolumn_namesifcnotingdf_zb.columns]
+ raiseValueError(f"Missing columns in gdf_zb: {missing}")
+ # get cell centers of cells to interpolate
+ xs,ys=da_mask.raster.xy(*np.where(da_mask.values))
+ cc=gpd.GeoDataFrame(geometry=gpd.points_from_xy(xs,ys),crs=da_mask.raster.crs)
+
+ # find nearest line and calculate relative distance along line for all z points
+ gdf_zb=gdf_zb[["geometry"]+column_names].copy()
+ gdf_zb["idx0"],gdf_zb["dist"]=nearest(gdf_zb,gdf_lines)
+ nearest_lines=gdf_lines.loc[gdf_zb["idx0"],"geometry"].values
+ gdf_zb["x"]=nearest_lines.project(gdf_zb["geometry"])
+ gdf_zb.set_index("idx0",inplace=True)
+ # keep only lines with associated points
+ gdf_lines=gdf_lines.loc[np.unique(gdf_zb.index.values)]
+
+ # find nearest line and calculate relative distance along line for all cell centers
+ cc["idx0"],cc["dist"]=nearest(cc,gdf_lines)
+ nearest_lines=gdf_lines.loc[cc["idx0"],"geometry"].values
+ cc["x"]=nearest_lines.project(cc["geometry"].to_crs(gdf_lines.crs))
+
+ # interpolate z values per line
+ def_interp(cc0,gdf_zb=gdf_zb,column_names=column_names):
+ kwargs=dict(kind="linear",fill_value="extrapolate")
+ idx0=cc0["idx0"].values[0]# iterpolate per line with ID idx0
+ fornameincolumn_names:
+ x0=np.atleast_1d(gdf_zb.loc[idx0,"x"])
+ z0=np.atleast_1d(gdf_zb.loc[idx0,name]).astype(np.float32)
+ valid=np.isfinite(z0)
+ x0,z0=x0[valid],z0[valid]
+ ifx0.size==0:
+ cc0[name]=np.nan
+ logger.warning(f"River segment {idx0} has no valid values for {name}.")
+ elifx0.size==1:
+ cc0[name]=z0[0]
+ else:
+ x1=cc0["x"].values
+ cc0[name]=interp1d(x0,z0,**kwargs)(x1)
+ returncc0
+
+ cc=cc.groupby("idx0").apply(_interp)[["geometry"]+column_names]
+
+ # rasterize interpolated z values
+ ds_out=xr.Dataset()
+ fornameincolumn_names:
+ da0=da_mask.raster.rasterize(cc,name,nodata=np.nan).where(da_mask)
+ ds_out=ds_out.assign(**{name:da0})
+
+ returnds_out
+
+
+
[docs]defburn_river_rect(
+ da_elv:xr.DataArray,
+ gdf_riv:gpd.GeoDataFrame,
+ da_man:xr.DataArray=None,
+ gdf_zb:gpd.GeoDataFrame=None,
+ gdf_riv_mask:gpd.GeoDataFrame=None,
+ segment_length:float=500,
+ riv_bank_q:float=0.5,# TODO add default value in docstring of setup_subgrid
+ rivwth_name:str="rivwth",
+ rivdph_name:str="rivdph",
+ rivbed_name:str="rivbed",
+ manning_name:str="manning",
+ logger=logger,
+):
+"""Burn rivers with a rectangular cross profile into a DEM.
+
+ Parameters
+ ----------
+ da_elv, da_man : xr.DataArray
+ DEM and manning raster to burn river depth and manning values into
+ gdf_riv : gpd.GeoDataFrame
+ River center lines.
+ gdf_zb : gpd.GeoDataFrame, optional
+ Point locations with a 'rivbed' river bed level [m+REF] column, by defualt None
+ gdf_riv_mask : gpd.GeoDataFrame, optional
+ Mask in which to interpolate z values, by default None.
+ If provided, 'rivwth' column is not required in `gdf_riv`.
+ segment_length : float, optional
+ Approximate river segment length [m], by default 500
+ riv_bank_q : float, optional
+ quantile [0-1] for river bank estimation, by default 0.25
+ rivwth_name, rivdph_name, rivbed_name, manning_name: str, optional
+ river width [m], depth [m], bed level [m+REF], & manning [s.m-1/3] column names
+ in gdf_riv, by default "rivwth", "rivdph", "rivbed", and "manning"
+
+
"""
- assertda_elv.raster.identical_grid(da_msk)
- logger.debug("Burn bedlevel values into DEM.")
- nodata=da_elv.raster.nodata
- # drop invalid segments, e.g. at pits
- gdf_riv=gdf_riv[np.logical_and(gdf_riv.length>0,np.isfinite(gdf_riv["zb"]))]
- zb=da_elv.raster.rasterize(gdf_riv,col_name="zb",nodata=-9999,**kwargs)
- # interpolate values if rivslp and rivdst is given
- ifnp.all(np.isin(["rivslp","rivdst"],gdf_riv.columns))andflwdirisnotNone:
- logger.debug("Interpolate bedlevel values")
- gdf_riv1=gdf_riv[gdf_riv["rivdst"]>0]
- slp=da_elv.raster.rasterize(gdf_riv1,col_name="rivslp",nodata=0)
- dst0=da_elv.raster.rasterize(gdf_riv1,col_name="rivdst",nodata=0)
- dst0=np.where(dst0>0,flwdir.distnc-dst0,0)
- zb=(zb+dst0*slp).where(zb!=-9999,-9999)
- zb.raster.set_nodata(-9999)
- zb.name="elevtn"
- # spread values inside river mask and replace da_elv values
- da_elv1=spread2d(zb,da_msk)["elevtn"]
- da_msk=np.logical_and(da_msk,da_elv1!=-9999)
- da_elv1=np.minimum(da_elv,da_elv1.where(da_msk,da_elv))
- da_elv1=da_elv1.where(da_elv!=nodata,nodata)
-
- ifriver_d4andflwdirisnotNone:
- logger.debug("Correct for D4 connectivity bed level")
- elevtn=flwdir.dem_dig_d4(da_elv1.values,rivmsk=da_msk.values,nodata=nodata)
- da_elv1=xr.DataArray(
- data=elevtn,
- coords=da_elv.raster.coords,
- dims=da_elv.raster.dims,
+ # clip
+ gdf_riv=gdf_riv.clip(da_elv.raster.box.to_crs(gdf_riv.crs))
+ ifgdf_riv.index.size==0:# no rivers in domain
+ returnda_elv,da_man
+
+ # reproject to utm if geographic
+ dst_crs=gdf_riv.crs
+ ifdst_crs.is_geographic:
+ dst_crs=parse_crs("utm",gdf_riv.to_crs(4326).total_bounds)
+ gdf_riv=gdf_riv.to_crs(dst_crs)
+
+ # check river mask
+ # create gdf_riv_mask based on buffered river center line only
+ # make sure the river is at least one cell wide
+ res=abs(da_elv.raster.res[0])
+ ifrivwth_nameingdf_riv.columns:
+ gdf_riv["buf"]=np.maximum(gdf_riv[rivwth_name].fillna(2),res)/2
+ else:
+ gdf_riv["buf"]=res/2
+ ifgdf_riv_maskisNone:
+ gdf_riv_mask=gdf_riv.assign(geometry=gdf_riv.buffer(gdf_riv["buf"]))
+ else:
+ # get gdf_riv outside of mask and buffer these lines
+ # then merge with gdf_riv_mask to get the full river mask
+ gdf_mask=gpd.GeoDataFrame(
+ geometry=[gdf_riv_mask.buffer(0).unary_union],
+ crs=gdf_riv_mask.crs,
+ )# create single polygon to clip
+ gdf_riv_clip=gdf_riv.overlay(gdf_mask,how="difference")
+ gdf_riv_mask1=gdf_riv_clip.assign(
+ geometry=gdf_riv_clip.buffer(gdf_riv_clip["buf"])
+ )
+ gdf_riv_mask=gpd.overlay(gdf_riv_mask,gdf_riv_mask1,how="union")
+ da_riv_mask=da_elv.raster.geometry_mask(gdf_riv_mask)
+
+ ifgdf_zbisNoneandrivbed_namenotingdf_riv.columns:
+ # calculate river bedlevel based on river depth per segment
+ gdf_riv_seg=split_line_equal(gdf_riv,segment_length).reset_index(drop=True)
+ # create mask or river bank cells adjacent to river mask -> numpy array
+ riv_mask_closed=ndimage.binary_closing(da_riv_mask)# remove islands
+ riv_bank=np.logical_xor(
+ riv_mask_closed,ndimage.binary_dilation(riv_mask_closed)
+ )
+ # make sure river bank cells are not nodata
+ riv_bank=np.logical_and(
+ riv_bank,np.isfinite(da_elv.raster.mask_nodata().values))
- da_msk=np.logical_or(da_msk,da_elv1<da_elv)
+ # sample elevation at river bank cells
+ rows,cols=np.where(riv_bank)
+ riv_bank_cc=gpd.GeoDataFrame(
+ geometry=gpd.points_from_xy(*da_elv.raster.xy(rows,cols)),
+ data={"z":da_elv.values[rows,cols]},
+ crs=da_elv.raster.crs,
+ )
+ # find nearest river center line for each river bank cell
+ riv_bank_cc["idx0"],_=nearest(riv_bank_cc,gdf_riv_seg)
+ # calculate segment river bank elevation as percentile of river bank cells
+ gdf_riv_seg["z"]=(
+ riv_bank_cc[["idx0","z"]].groupby("idx0").quantile(q=riv_bank_q)
+ )
+ # calculate river bed elevation per segment
+ gdf_riv_seg[rivbed_name]=gdf_riv_seg["z"]-gdf_riv_seg[rivdph_name]
+ # get zb points at center of line segments
+ points=gdf_riv_seg.geometry.interpolate(0.5,normalized=True)
+ gdf_zb=gdf_riv_seg.assign(geometry=points)
+ elifgdf_zbisNone:
+ # get zb points at center of line segments
+ points=gdf_riv.geometry.interpolate(0.5,normalized=True)
+ gdf_zb=gdf_riv.assign(geometry=points)
+ elifgdf_zbisnotNone:
+ ifrivbed_namenotingdf_zb.columns:
+ raiseValueError(f"Missing {rivbed_name} attribute in gdf_zb")
+ # fill missing manning values based on centerlines
+ # TODO manning always defined on centerline?
+ ifmanning_namenotingdf_zb.columnsandmanning_nameingdf_riv.columns:
+ gdf_zb[manning_name]=np.nan
+ ifnp.any(np.isnan(gdf_zb[manning_name])):
+ gdf_zb["idx0"],_=nearest(gdf_zb,gdf_riv)
+ man_nearest=gdf_riv.loc[gdf_zb["idx0"],manning_name]
+ gdf_zb[manning_name]=gdf_zb[manning_name].fillna(man_nearest)
+ elifrivbed_namenotingdf_zb.columns:
+ raiseValueError(f"Missing {rivbed_name} or {rivdph_name} attributes")
+
+ # merge river lines > z points are interpolated along merged line
+ ifgdf_riv.index.size>1:
+ gdf_riv_merged=gpd.GeoDataFrame(
+ geometry=[linemerge(gdf_riv.unary_union)],crs=gdf_riv.crs
+ )
+ gdf_riv_merged=gdf_riv_merged.explode().reset_index(drop=True)
+ else:
+ gdf_riv_merged=gdf_riv
+
+ # interpolate river depth and manning along river center line
+ # TODO nearest interpolation for manning?
+ column_names=[rivbed_name]
+ ifmanning_nameingdf_zb.columns:
+ column_names+=[manning_name]
+ fornameincolumn_names:
+ gdf_zb=gdf_zb[np.isfinite(gdf_zb[name])]
+ ds=interp_along_line_to_grid(
+ da_mask=da_riv_mask,
+ gdf_lines=gdf_riv_merged,
+ gdf_zb=gdf_zb,
+ column_names=column_names,
+ logger=logger,
+ )
+
+ # update elevation with river bottom elevations
+ # river bed elevation must be lower than original elevation
+ da_elv1=da_elv.where(
+ np.logical_or(np.isnan(ds[rivbed_name]),da_elv<ds[rivbed_name]),
+ ds[rivbed_name],
+ )
+
+ # update manning:
+ da_man1=da_man
+ ifmanning_nameindsandda_manisnotNone:
+ da_man1=da_man1.where(np.isnan(ds[manning_name]),ds[manning_name])
- # set attrs and return
- da_elv1.raster.set_nodata(nodata)
- da_elv1.raster.set_crs(da_elv.raster.crs)
- returnda_elv1,da_msk
[docs]defmerge_multi_dataarrays(da_list:List[dict],
+ gdf_list:List[dict]=[],da_like:xr.DataArray=None,reproj_kwargs:Dict={},buffer_cells:int=0,# not in list
@@ -555,6 +558,22 @@
Source code for hydromt_sfincs.workflows.merge
interp_method=interp_method,)
+ # burn in rivers
+ foriinrange(len(gdf_list)):
+ cs_type=gdf_list[i].get("type","rectangular")
+ ifcs_type=="rectangular":
+ # width or gdf_riv_mask is required
+ # zb is used if provided, otherwise depth is used
+ da1=burn_river_rect(
+ da_elv=da1,
+ gdf_riv=gdf_list[i].get("gdf"),
+ gdf_riv_mask=gdf_list[i].get("gdf_mask",None),
+ rivdph_name=gdf_list[i].get("depth","rivdph"),
+ rivwth_name=gdf_list[i].get("width","rivwth"),
+ rivbed_name=gdf_list[i].get("zb","rivbed"),
+ )
+ else:
+ raiseNotImplementedError(f"Cross section type {cs_type} not implemented.")returnda1