From 0072b0915aa89b44862622662e56d318d480c68e Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Mon, 18 Mar 2024 19:55:31 +0100
Subject: [PATCH] No need for MeshKernel

---
 coupling/delwaq/delwaq_util.py | 43 +++-------------------------------
 coupling/delwaq/gen_delwaq.py  | 13 +++-------
 2 files changed, 6 insertions(+), 50 deletions(-)

diff --git a/coupling/delwaq/delwaq_util.py b/coupling/delwaq/delwaq_util.py
index 180253953..5c0e323d8 100644
--- a/coupling/delwaq/delwaq_util.py
+++ b/coupling/delwaq/delwaq_util.py
@@ -90,6 +90,8 @@ def ugridify(model: ribasim.Model):
         model.filepath.parent / "database.gpkg", layer="Node", fid_as_index=True
     )
     edge_df = model.edge.df[model.edge.df.edge_type == "flow"]
+    edge_df.set_crs(epsg=28992, inplace=True, allow_override=True)
+    node_df.set_crs(epsg=28992, inplace=True, allow_override=True)
 
     # from node_id to the node_dim index
     node_lookup = pd.Series(
@@ -97,12 +99,6 @@ def ugridify(model: ribasim.Model):
         data=node_df.index.argsort(),
         name="node_index",
     )
-    # from edge_id to the edge_dim index
-    edge_lookup = pd.Series(
-        index=edge_df.index.rename("edge_id"),
-        data=edge_df.index.argsort(),
-        name="edge_index",
-    )
 
     grid = xu.Ugrid1d(
         node_x=node_df.geometry.x,
@@ -127,38 +123,5 @@ def ugridify(model: ribasim.Model):
     uds = uds.assign_coords(from_node_id=(edge_dim, edge_df.from_node_id))
     uds = uds.assign_coords(to_node_id=(edge_dim, edge_df.to_node_id))
     uds = uds.assign_coords(edge_id=(edge_dim, edge_df.index))
-    # MDAL doesn't like string coordinates
-    # uds = uds.assign_coords(node_name=(node_dim, node_df.name))
-    # uds = uds.assign_coords(node_type=(node_dim, node_df["type"]))
-    # uds = uds.assign_coords(edge_name=(edge_dim, edge_df.name))
-
-    results_dir = model.filepath.parent / model.results_dir
-
-    # Split out the boundary condition flows, since are on nodes, not edges.
-    # Use pyarrow backend since it doesn't convert the edge_id to float to handle missings.
-    all_flow_df = pd.read_feather(results_dir / "flow.arrow")
-
-    # https://github.com/pydata/xarray/issues/6318 datetime64[ms] gives trouble
-    all_flow_df.time = all_flow_df.time.astype("datetime64[ns]")
-
-    flow_df = all_flow_df[all_flow_df.edge_id.notna()].copy()
-    # The numpy_nullable backend converts to float to handle missing,
-    # now we can convert it back since there are no missings left here.
-    # The pyarrow backend fixes this, but then we get object dtypes after to_xarray()
-    flow_df.edge_id = flow_df.edge_id.astype(np.int64)
-
-    flow_df[edge_dim] = edge_lookup[flow_df.edge_id].to_numpy()
-    flow_da = flow_df.set_index(["time", edge_dim]).flow_rate.to_xarray()
-    uds["flow"] = flow_da
-
-    bc_flow_df = all_flow_df[all_flow_df.edge_id.isna()].copy()
-    bc_flow_df = bc_flow_df.rename(
-        columns={"flow_rate": "boundary_flow", "from_node_id": "node_id"}
-    ).drop(columns=["edge_id", "to_node_id"])
-    bc_flow_df[node_dim] = node_lookup[bc_flow_df.node_id].to_numpy()
-    bc_flow_df
-    bc_flow_da = bc_flow_df.set_index(["time", node_dim]).boundary_flow.to_xarray()
-    # perhaps not the best name?
-    # this does not visualize properly, and is selected on load, maybe due to alphabetical order
-    uds["boundary_flow"] = bc_flow_da
+
     return uds
diff --git a/coupling/delwaq/gen_delwaq.py b/coupling/delwaq/gen_delwaq.py
index f508156d9..ad045e57c 100644
--- a/coupling/delwaq/gen_delwaq.py
+++ b/coupling/delwaq/gen_delwaq.py
@@ -6,11 +6,9 @@
 from pathlib import Path
 
 import geopandas as gpd
-import meshkernel as mk
 import numpy as np
 import pandas as pd
 import ribasim
-import xugrid as xu
 from delwaq_util import (
     strfdelta,
     ugridify,
@@ -50,6 +48,7 @@
 edge = model.edge.df
 assert edge is not None
 edge = edge[edge.edge_type == "flow"]  # no control or allocation stuff please
+edge.set_crs(epsg=28992, inplace=True, allow_override=True)
 assert edge is not None
 
 # Flows on non-existing edges indicate where the boundaries are
@@ -84,14 +83,8 @@
     )
 
 # Generate mesh and write to NetCDF
-edges = np.array(list(zip(edge.from_node_id, edge.to_node_id))).flatten()
-mesh1d = mk.Mesh1d(node.geometry.x, node.geometry.y, edges)
-ugrid = xu.Ugrid1d.from_meshkernel(mesh1d)
-ugrid.set_crs(epsg=28992)
-ds = ugrid.to_dataset()
-# ds.to_netcdf("ribasim.nc")
-ds = ugridify(model)
-ds.ugrid.to_netcdf(output_folder / "ribasim.nc")
+uds = ugridify(model)
+uds.ugrid.to_netcdf(output_folder / "ribasim.nc")
 
 # Generate area and flows
 # File format is int32, float32 based