Skip to content

Commit

Permalink
solved conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
veenstrajelmer committed Nov 1, 2023
2 parents 2744f2b + c5c2d03 commit 395fa18
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 50 deletions.
63 changes: 38 additions & 25 deletions hydrolib/core/dflowfm/net/models.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ def split_by(gl: mk.GeometryList, by: float) -> list:


# TODO: this is an inconvenient name since meshkernel also has a Mesh2d class
class Mesh2d(BaseModel):
class Mesh2d(BaseModel):
"""Mesh2d defines a single two dimensional grid.
Attributes:
Expand Down Expand Up @@ -120,48 +120,52 @@ class Mesh2d(BaseModel):
# mesh2d_face_nodes: np.ndarray = Field(
# default_factory=lambda: np.empty((0, 0), dtype=np.int32)
# )

@property
def mesh2d_node_x(self):
return self.meshkernel.mesh2d_get().node_x

@property
def mesh2d_node_y(self):
return self.meshkernel.mesh2d_get().node_y

@property
def mesh2d_edge_x(self):
return self.meshkernel.mesh2d_get().edge_x

@property
def mesh2d_edge_y(self):
return self.meshkernel.mesh2d_get().edge_y

@property
def mesh2d_edge_nodes(self):
mesh2d_output = self.meshkernel.mesh2d_get()
edge_nodes = mesh2d_output.edge_nodes.reshape((-1, 2))
return edge_nodes

@property
def mesh2d_face_x(self):
return self.meshkernel.mesh2d_get().face_x

@property
def mesh2d_face_y(self):
return self.meshkernel.mesh2d_get().face_y

@property
def mesh2d_face_nodes(self):
mesh2d_output = self.meshkernel.mesh2d_get()
npf = mesh2d_output.nodes_per_face
if self.is_empty():
return np.empty((0, 0), dtype=np.int32)
face_node_connectivity = np.full((len(mesh2d_output.face_x), max(npf)), np.iinfo(np.int32).min)
idx = (np.ones_like(face_node_connectivity) * np.arange(max(npf))[None, :]) < npf[:, None]
face_node_connectivity = np.full(
(len(mesh2d_output.face_x), max(npf)), np.iinfo(np.int32).min
)
idx = (
np.ones_like(face_node_connectivity) * np.arange(max(npf))[None, :]
) < npf[:, None]
face_node_connectivity[idx] = mesh2d_output.face_nodes
return face_node_connectivity

def is_empty(self) -> bool:
"""Determine whether this Mesh2d is empty.
Expand All @@ -185,7 +189,7 @@ def _set_mesh2d(self, node_x, node_y, edge_nodes) -> None:
node_y=node_y.astype(np.float64),
edge_nodes=edge_nodes.ravel().astype(np.int32),
)

self.meshkernel.mesh2d_set(mesh2d)

def get_mesh2d(self) -> mk.Mesh2d:
Expand Down Expand Up @@ -234,7 +238,7 @@ def create_triangular(self, geometry_list: mk.GeometryList) -> None:
"""
# Call meshkernel
self.meshkernel.mesh2d_make_triangular_mesh_from_polygon(geometry_list)

def clip(
self,
geometrylist: mk.GeometryList,
Expand All @@ -247,7 +251,7 @@ def clip(
geometrylist (GeometryList): Polygon stored as GeometryList
deletemeshoption (int, optional): [description]. Defaults to 1.
"""

# For clipping outside
if not inside:
# Check if a multipolygon was provided when clipping outside
Expand Down Expand Up @@ -309,7 +313,7 @@ def refine(self, polygon: mk.GeometryList, level: int, min_edge_size=10.0):
polygon (GeometryList): Polygon in which to refine
level (int): Number of refinement steps
"""

# Check if parts are closed
# if not (polygon.x_coordinates[0], polygon.y_coordinates[0]) == (
# polygon.x_coordinates[-1],
Expand Down Expand Up @@ -729,8 +733,12 @@ class Mesh1d(BaseModel):
network1d_edge_nodes: np.ndarray = Field(
default_factory=lambda: np.empty((0, 2), np.int32)
)
network1d_geom_x: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double)) #TODO: sync with meshkernel.node_x
network1d_geom_y: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double)) #TODO: sync with meshkernel.node_y
network1d_geom_x: np.ndarray = Field(
default_factory=lambda: np.empty(0, np.double)
) # TODO: sync with meshkernel.node_x
network1d_geom_y: np.ndarray = Field(
default_factory=lambda: np.empty(0, np.double)
) # TODO: sync with meshkernel.node_y
network1d_part_node_count: np.ndarray = Field(
default_factory=lambda: np.empty(0, np.int32)
)
Expand All @@ -748,8 +756,10 @@ class Mesh1d(BaseModel):
default_factory=lambda: np.empty(0, np.double)
)

mesh1d_edge_nodes: np.ndarray = Field( #TODO: sync with meshkernel.edge_nodes (ravelled)
default_factory=lambda: np.empty((0, 2), np.int32)
mesh1d_edge_nodes: np.ndarray = (
Field( # TODO: sync with meshkernel.edge_nodes (ravelled)
default_factory=lambda: np.empty((0, 2), np.int32)
)
)
mesh1d_edge_x: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double))
mesh1d_edge_y: np.ndarray = Field(default_factory=lambda: np.empty(0, np.double))
Expand Down Expand Up @@ -1210,16 +1220,19 @@ def mesh1d_add_branch(
)
self._mesh1d._set_mesh1d()
return name

def plot(self):
import matplotlib.pyplot as plt
fig,ax = plt.subplots()

fig, ax = plt.subplots()
mesh2d_output = self._mesh2d.get_mesh2d()
mesh1d_output = self._mesh1d._get_mesh1d()
links_output = self._link1d2d.meshkernel.contacts_get()
mesh2d_output.plot_edges(ax=ax, color='r')
mesh1d_output.plot_edges(ax=ax, color='g')
links_output.plot_edges(ax=ax, mesh1d=mesh1d_output, mesh2d=mesh2d_output, color='k')
mesh2d_output.plot_edges(ax=ax, color="r")
mesh1d_output.plot_edges(ax=ax, color="g")
links_output.plot_edges(
ax=ax, mesh1d=mesh1d_output, mesh2d=mesh2d_output, color="k"
)


class NetworkModel(ParsableFileModel):
Expand Down
10 changes: 5 additions & 5 deletions hydrolib/core/dflowfm/net/reader.py
Original file line number Diff line number Diff line change
Expand Up @@ -86,17 +86,17 @@ def read_mesh2d(self, mesh2d: Mesh2d) -> None:
# for meshkey, nckey in self._explorer.mesh2d_var_name_mapping.items():
# setattr(mesh2d, meshkey, self._read_nc_attribute(ds[nckey]))
# TODO: replace with xugrid reader?

# set mesh2d on meshkernel instance
node_x = self._read_nc_attribute(ds["mesh2d_node_x"])
node_y = self._read_nc_attribute(ds["mesh2d_node_y"])
edge_nodes = self._read_nc_attribute(ds["mesh2d_edge_nodes"])
mesh2d._set_mesh2d(node_x=node_x, node_y=node_y, edge_nodes=edge_nodes)

# bathymetry
node_z = self._read_nc_attribute(ds["mesh2d_node_z"])
mesh2d.mesh2d_node_z = node_z

ds.close()

def read_link1d2d(self, link1d2d: Link1d2d) -> None:
Expand All @@ -116,14 +116,14 @@ def read_link1d2d(self, link1d2d: Link1d2d) -> None:
setattr(link1d2d, meshkey, self._read_nc_attribute(ds[nckey]))

# TODO: setting contacts is not possible in meshkernel (e.g. with contacts_set())
# so misalignment between link1d2d.link1d2d and
# so misalignment between link1d2d.link1d2d and
# empty _link1d2d.meshkernel.contacts_get().mesh2d_indices
# mesh1d_indices = link1d2d.link1d2d[:,0]
# mesh2d_indices = link1d2d.link1d2d[:,1]
# import meshkernel as mk
# contacts = mk.Contacts(mesh1d_indices=mesh1d_indices, mesh2d_indices=mesh2d_indices)
# link1d2d.meshkernel.contacts_set(contacts)

ds.close()

def _read_nc_attribute(self, attr: nc._netCDF4.Variable) -> np.ndarray:
Expand Down
42 changes: 22 additions & 20 deletions tests/dflowfm/test_net.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,15 +17,17 @@
from ..utils import test_input_dir, test_output_dir


def plot_network(network): #TODO: can be removed
def plot_network(network): # TODO: can be removed
_, ax = plt.subplots()
ax.set_aspect(1.0)
network.plot(ax=ax) #TODO: newly added but was already called
network.plot(ax=ax) # TODO: newly added but was already called
ax.autoscale()
plt.show()


def _plot_mesh2d(mesh2d, ax=None, **kwargs): #TODO, this can be removed meshkernel.mesh2d_get().plot_edges() does the trick
def _plot_mesh2d(
mesh2d, ax=None, **kwargs
): # TODO, this can be removed meshkernel.mesh2d_get().plot_edges() does the trick
from matplotlib.collections import LineCollection

if ax is None:
Expand Down Expand Up @@ -166,14 +168,14 @@ def test_create_1d_2d_1d2d():
assert network1_link1d2d.shape == (21, 2)
assert network1_con_m1d.size == 21
assert network1_con_m2d.size == 21

# Write to file
file_out = Path(test_output_dir / "test_net.nc")
network.to_file(file_out)
#read from written file

# read from written file
network2 = NetworkModel(file_out)

mesh2d_output = network2._mesh2d.get_mesh2d()
assert len(network2._mesh2d.mesh2d_face_x) == 152
mesh1d_output = network2._mesh1d._get_mesh1d()
Expand All @@ -183,10 +185,10 @@ def test_create_1d_2d_1d2d():
network2_con_m1d = network2._link1d2d.meshkernel.contacts_get().mesh1d_indices
network2_con_m2d = network2._link1d2d.meshkernel.contacts_get().mesh2d_indices
assert network2_link1d2d.shape == (21, 2)
#TODO: below asserts fail, since the meshkernel contacts are not set upon reading (not implemented in meshkernel)
# TODO: below asserts fail, since the meshkernel contacts are not set upon reading (not implemented in meshkernel)
# assert network2_con_m1d.size == 21
# assert network2_con_m2d.size == 21

# plot both networks
network.plot()
network2.plot()
Expand Down Expand Up @@ -332,25 +334,25 @@ def test_read_write_read_compare_nodes(filepath):

# Create network model
network1 = NetworkModel(filepath=filepath)

# Save to temporary location
save_path = (
test_output_dir
/ "test_read_write_read_compare"#.__name__
/ "test_read_write_read_compare" # .__name__
/ network1._generate_name()
)
network1.save(filepath=save_path)

# # Read a second network from this location
network2 = NetworkModel(filepath=network1.filepath)

network1_mesh1d_node_x = network1._mesh1d._get_mesh1d().node_x
network1_mesh2d_node_x = network1._mesh2d.get_mesh2d().node_x
network2_mesh1d_node_x = network2._mesh1d._get_mesh1d().node_x
network2_mesh2d_node_x = network2._mesh2d.get_mesh2d().node_x

netw1_nnodes = len(network1_mesh1d_node_x) + len(network1_mesh2d_node_x)

assert netw1_nnodes > 0
assert (network1_mesh1d_node_x == network2_mesh1d_node_x).all()
assert (network1_mesh2d_node_x == network2_mesh2d_node_x).all()
Expand Down Expand Up @@ -480,13 +482,13 @@ def test_read_net_nc_2d_without_faces(self):
network = NetworkModel(filepath=filepath)
assert network._mesh1d.is_empty()
assert not network._mesh2d.is_empty()

mesh2d_output = network._mesh2d.get_mesh2d()
#TODO: not empty anymore since meshkernel creates the face_node_connectivity

# TODO: not empty anymore since meshkernel creates the face_node_connectivity
assert len(mesh2d_output.face_x) == 208
assert len(mesh2d_output.face_y) == 208
assert len(network._mesh2d.mesh2d_face_z) == 0 #TODO: update face_z length
assert len(network._mesh2d.mesh2d_face_z) == 0 # TODO: update face_z length
assert network._mesh2d.mesh2d_face_nodes.shape == (208, 4)

assert len(mesh2d_output.node_x) == 238
Expand Down Expand Up @@ -646,7 +648,7 @@ def test_create_triangular():
)

network.mesh2d_create_triangular_within_polygon(polygon)

assert np.array_equiv(
network._mesh2d.mesh2d_node_x,
np.array([6.0, 4.0, 2.0, 0.0]),
Expand Down Expand Up @@ -676,7 +678,7 @@ def test_add_1d2d_links():
m2d_nnodes = network._mesh2d.meshkernel.mesh2d_get().node_x.size
assert m1d_nnodes == 12
assert m2d_nnodes == 121

# Get required arguments
node_mask = network._mesh1d.get_node_mask([branchid])
exterior = GeometryList(
Expand Down

0 comments on commit 395fa18

Please sign in to comment.