@dprincz suggested to add an element for the "outlet" channels to make them having a non-zero values of NEXT the next variable.

          I'm not sure `NAA=NA` should be an assumption set by the code. It could be set by the fields in the drainage database (e.g., by having no cell with `NEXT==0`, which could be handled by the workflow.

Originally posted by @dprincz in MESH-Model/MESH-Dev#36 (comment)

The MESH_parameters.nc file can now be prepared and fed into MESH >r1860 versions.

Reported by: @fuadyassin

Hello everybody,
Here is an error that occurred while running mesh_flow.ipynb (see below). It seems the problem comes from different dimensions.
Please find attached the files used in my case (domain shapefiles, modified landcover file, 3 remapped .nc files).

Thank you for your help

Adrien
files_east.zip

Here is the error message :

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 exp1.run()

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/meshflow/core.py:320, in MESHWorkflow.run(self)
    317 _ordered_dims = {self.main_id: self.main_seg}
    319 # 1 generate drainage database
--> 320 self.ddb = utility.prepare_mesh_ddb(riv=self.reordered_riv,
    321                                     cat=self.cat,
    322                                     landcover=self.landcover,
    323                                     cat_dim=self.main_id,
    324                                     gru_dim=self.gru_dim,
    325                                     gru_var=self.gru_var,
    326                                     hru_dim=self.hru_dim,
    327                                     hru_var=self.hru_var,
    328                                     gru_names=self.landcover_classes,
    329                                     include_vars=self.ddb_vars,
    330                                     attr_local=self.ddb_local_attrs,
    331                                     attr_global=self.ddb_global_attrs,
    332                                     min_values=self.ddb_min_values,
    333                                     fill_na=None,
    334                                     ordered_dims=_ordered_dims,
    335                                     ddb_units=self.ddb_units,
    336                                     ddb_to_units=self.ddb_to_units)
    338 # forcing unit registry
    339 _ureg = pint.UnitRegistry(force_ndarray_like=True)

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/meshflow/utility/network.py:452, in prepare_mesh_ddb(riv, cat, landcover, cat_dim, gru_dim, gru_var, hru_dim, hru_var, gru_names, include_vars, attr_local, attr_global, min_values, fill_na, ordered_dims, ddb_units, ddb_to_units)
    436 # creating dictionary for final merging
    437 lc_names = {
    438     "coords": {
    439         gru_dim: {
   (...)
    450     },
    451 }
--> 452 lc_names_da = xr.Dataset.from_dict(lc_names)
    454 # making a list of all xarray.Dataset objects
    455 ds_list = [landcover_ds,
    456            lc_names_da,
    457            riv_ds,
    458            cat_ds]

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/xarray/core/dataset.py:7438, in Dataset.from_dict(cls, d)
   7433 except KeyError as e:
   7434     raise ValueError(
   7435         "cannot convert dict without the key "
   7436         "'{dims_data}'".format(dims_data=str(e.args[0]))
   7437     )
-> 7438 obj = cls(variable_dict)
   7440 # what if coords aren't dims?
   7441 coords = set(d.get("coords", {})) - set(d.get("dims", {}))

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/xarray/core/dataset.py:685, in Dataset.__init__(self, data_vars, coords, attrs)
    682 if isinstance(coords, Dataset):
    683     coords = coords._variables
--> 685 variables, coord_names, dims, indexes, _ = merge_data_and_coords(
    686     data_vars, coords
    687 )
    689 self._attrs = dict(attrs) if attrs is not None else None
    690 self._close = None

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/xarray/core/dataset.py:416, in merge_data_and_coords(data_vars, coords)
    412     coords = create_coords_with_default_indexes(coords, data_vars)
    414 # exclude coords from alignment (all variables in a Coordinates object should
    415 # already be aligned together) and use coordinates' indexes to align data_vars
--> 416 return merge_core(
    417     [data_vars, coords],
    418     compat="broadcast_equals",
    419     join="outer",
    420     explicit_coords=tuple(coords),
    421     indexes=coords.xindexes,
    422     priority_arg=1,
    423     skip_align_args=[1],
    424 )

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/xarray/core/merge.py:723, in merge_core(objects, compat, join, combine_attrs, priority_arg, explicit_coords, indexes, fill_value, skip_align_args)
    718 prioritized = _get_priority_vars_and_indexes(aligned, priority_arg, compat=compat)
    719 variables, out_indexes = merge_collected(
    720     collected, prioritized, compat=compat, combine_attrs=combine_attrs
    721 )
--> 723 dims = calculate_dimensions(variables)
    725 coord_names, noncoord_names = determine_coords(coerced)
    726 if explicit_coords is not None:

File ~/FHIMP/virtual_env/fhimp_venv/lib/python3.10/site-packages/xarray/core/variable.py:3257, in calculate_dimensions(variables)
   3255             last_used[dim] = k
   3256         elif dims[dim] != size:
-> 3257             raise ValueError(
   3258                 f"conflicting sizes for dimension {dim!r}: "
   3259                 f"length {size} on {k!r} and length {dims[dim]} on {last_used!r}"
   3260             )
   3261 return dims

ValueError: conflicting sizes for dimension 'gru': length 1 on 'landcover_names' and length 0 on {'gru': 'gru'}

There seems to be an issue with the re-ordering of Rank and Next variables. It seems to be because of the "extra" variables added to the ddb_vars dictionary of the MESHWorkflow class.

While this issue needs to be taken care of in the longer term, but for now, users must avoid adding extra variable to the ddb, until this is facilitated in the workflow.

Reported-by: Sujata Budhathoki @sujata91

It is better to change the places, as it will be inline with the relevant dictionary of ddb.

Reported by: @sujata91

Mountain MESH workflows are missing plus their documentations.

Requested by: @mee067

The specifications of acceptable ESRI Shapefiles fed to the workflow is missing. It must be stated that:

1. Both rivers and catchments ESRI Shapefiles must have equal sizes, and if not there are certain assumptions,
2. There is a common ID data shared between the two ESRI Shapefiles,
3. The units of each Shapefile is a bit ambiguous, and it causes confusion for the users.

This is reported by: @sujata91

Requested by @ShervanGharari

It would be great if MESH setups could be subsetted into smaller domains. This can be implemented into meshflow.

I encountered a memory issue in generating the forcing file required for MESH using the meshflow workflow. The workflow runs well for smaller timeframe and smaller number of subbasin. But there is issue with the memory with larger forcing dataset. Previously I had ~1600 subbasin and ran from 1980-2018 forcing files and had to provide like 128GB memory specs. With another geofabric (waterloo Merithydro geofabric) there are ~4000 subbasin and same time frame. Providing larger memory makes the job pending in cluster and even if it runs there is a memory issue. May be there could be some alternative solution to this in the workflow?

Suggested by @MIsmlAhmed

It would be interesting to consider a threshold for fractions of each GRU to be finally considered in the final drainage database of MESH setups.

The community is welcome to comment and suggest.

There is an issue with MESH where the ddb dimension name is not consistent between various official versions, therefore, some flexibility is needed from the user side to change the dimensions names of the ddb object by the configuration dictionary.

The explanation is detailed in MESH-Model/MESH-Dev#6.

Possible solutions for the problem itself is using GRU and NGRU as the variables names instead of their lowercase forms.

This should be fixed with the next version of MESHFlow.

Problem statement

A warning is given as following:

/home/kasra545/virtual-envs/scienv/lib/python3.10/site-packages/xarray/core/accessor_dt.py:72: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self.
  values_as_series = pd.Series(values.ravel(), copy=False)

Unsure where it is coming from, so it warrants further investigations.

After hitting .run() on the MESHWorkFlow class, one gets the following warning:

path/to/meshflow/src/meshflow/utility/forcing_prep.py:116: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  var = [i for i in ds.dims.keys() if i != 'time']

The current version of MESH r1860 does not route streamflow in outlet sub-basins/channels. Therefore, as a quick fix, a fictitious channel could be added to each outlet sub-basin/channel, to facilitate the routing. Since these fictitious channels are considered as "outlet" instead, no routing will be done on them either and modellers could trust the result of vector-based routing.

Copying @mee067

Mohamed's comments on Slack:

for gridded setups, we used to have one more grid to account for that. GK does it automatically. It has a Rank with Next = 0 Earlier versions of the code used not to even run CLASS on that grid but newer code (1860) does run CLASS for that grid. Moving to vector-based setups, CLASS would thus run on the last subbasin which probably is not an extra basin. runoff will be routed to the instream channel but I do not think there will additional channel routing for the flow - I am not 100% sure though. I do not think that's critical. It does produce output at the outlet.

The thing is, MESH needs a NEXT value to route to, and assumes the channel to go from a grid/sub to its NEXT one. It's not from the most u/s point within a grid/sub to its most d/s point. I don't think there is channel routing for the last sub and thus we need to add another fictitious one d/s to get the effect. How much difference would that make for the outlet flow, I can't say, but I don't think it's that large.

The vector workflow does not add that extra sub and therefore, I found that Qo = Qi for the last sub for the vector-based Yukon setup. Qo is the channel outflow, Qi is the inflow from u/s subs + local RFF. I made a quick comparison for the Yukon outlet and I do not feel the local runoff is even included. I see all loops go only to NAA and then this is a special treatment for the outlet on lines 958-960 of rte_module.f90:

https://github.com/MESH-Model/MESH-Dev/blob/ff1de25e67ef552844f1e789c0faffb4175a8dac/Routing_Model/RPN_watroute/sa_mesh_process/rte_module.f90#L957-L959

For large basins such as the MRB, YRB, this may not be important - runoff from the last sub/grid would not be that large compared to the huge flow collected from the basin

but for small basins - like those little basins draining to the sea or a lake, ignoring the runoff and routing within the last sub could be significant - depends on the size and relative contribution of that last sub to the flow of the basin. It could be raining a lot on the shore line of a marine basin or a lake shore to make significant contribution.

If available, MESH can understand the river classes defined for a setup. For example, in MERIT-Basins, the order value could be used as the IAK in the drainage database. This could be informative for new users using this workflow.

Reported by @ShervanGharari

It is recommended to separate the processes to generate the forcing and ddb files.

This secondary check on the landcover fractions is currently missing.

Reported by: @ShervanGharari

In some setups, if the forcing files are created on various dates (I think this is the issue), meshflow cannot understand different files and sticks with the old ones.

This is a very rough description, as I do not understand the problem in the first place.

Reported-by: Jongho Keum (McMaster University) and Adrien Pierre (Laval University)