import re
import warnings
import numpy as np
import xarray as xr
from matplotlib.path import Path
from requireit import require_array
from requireit import require_length_at_least
from scipy.spatial import Delaunay
from scipy.spatial import Voronoi
from landlab.core.utils import as_id_array
from landlab.graph.sort.intpair import IntPairCollection
from landlab.graph.sort.intpair import IntPairMapping
from landlab.graph.sort.sort import reverse_one_to_one
from landlab.utils import jaggedarray
[docs]
class VoronoiDelaunay:
[docs]
def __init__(self, xy_of_node):
# What we need:
# * [x] xy_of_node
# * [x] nodes_at_link
# * [x] links_at_patch
# And then for the dual graph:
# * [x] xy_of_corner
# * [x] corners_at_face
# * [x] faces_at_cell
# And the to link the graphs:
# * [x] node_at_cell
# * [x] nodes_at_face
# points == xy_of_node
# vertices == xy_of_corner
# regions == corners_at_cell
# ridge_vertices == corners_at_face
# ridge_points == nodes_at_face
# point_region == node_at_cell
delaunay = Delaunay(xy_of_node)
voronoi = Voronoi(xy_of_node)
voronoi.regions, voronoi.point_region = VoronoiDelaunay._remove_empty_regions(
voronoi.regions, voronoi.point_region
)
mesh = xr.Dataset(
{
"node": xr.DataArray(
data=np.arange(len(voronoi.points)),
coords={
"x_of_node": xr.DataArray(voronoi.points[:, 0], dims=("node",)),
"y_of_node": xr.DataArray(voronoi.points[:, 1], dims=("node",)),
},
dims=("node",),
),
"corner": xr.DataArray(
data=np.arange(len(voronoi.vertices)),
coords={
"x_of_corner": xr.DataArray(
voronoi.vertices[:, 0], dims=("corner",)
),
"y_of_corner": xr.DataArray(
voronoi.vertices[:, 1], dims=("corner",)
),
},
dims=("corner",),
),
}
)
mesh.update(
{
"nodes_at_link": xr.DataArray(
as_id_array(voronoi.ridge_points), dims=("link", "Two")
),
"nodes_at_patch": xr.DataArray(
np.asarray(delaunay.simplices, dtype=int), dims=("patch", "Three")
),
"corners_at_face": xr.DataArray(
voronoi.ridge_vertices, dims=("face", "Two")
),
"corners_at_cell": xr.DataArray(
self._corners_at_cell(voronoi.regions),
dims=("cell", "max_corners_per_cell"),
),
"n_corners_at_cell": xr.DataArray(
[len(cell) for cell in voronoi.regions], dims=("cell",)
),
"nodes_at_face": xr.DataArray(
np.asarray(voronoi.ridge_points, dtype=int), dims=("face", "Two")
),
"cell_at_node": xr.DataArray(
np.asarray(voronoi.point_region, dtype=int), dims=("node",)
),
}
)
self._mesh = mesh
@staticmethod
def _corners_at_cell(regions):
jagged = jaggedarray.JaggedArray(regions)
return np.asarray(
jaggedarray.unravel(jagged.array, jagged.offset, pad=-1), dtype=int
)
@staticmethod
def _remove_empty_regions(regions, point_region):
"""Remove regions that have no points.
Parameters
----------
regions : list of list of int
Lists of each region's vertices.
point_regions : list in int
The region associated with each point.
Returns
-------
regions, point_regions
Copies of the input arrays with empty regions removed.
"""
size_of_region = np.array([len(region) for region in regions], dtype=int)
empty_regions = np.where(size_of_region == 0)[0]
if len(empty_regions):
point_region = np.asarray(point_region, dtype=int)
regions = list(regions)
for region in empty_regions[::-1]:
regions.pop(region)
point_region[point_region >= region] -= 1
return regions, point_region
@property
def number_of_nodes(self):
return self._mesh.sizes["node"]
@property
def number_of_links(self):
return self._mesh.sizes["link"]
@property
def number_of_patches(self):
return self._mesh.sizes["patch"]
@property
def number_of_corners(self):
return self._mesh.sizes["corner"]
@property
def number_of_faces(self):
return self._mesh.sizes["face"]
@property
def number_of_cells(self):
return self._mesh.sizes["cell"]
@property
def x_of_node(self):
return self._mesh["x_of_node"].values
@property
def y_of_node(self):
return self._mesh["y_of_node"].values
@property
def x_of_corner(self):
return self._mesh["x_of_corner"].values
@property
def y_of_corner(self):
return self._mesh["y_of_corner"].values
@property
def nodes_at_patch(self):
return self._mesh["nodes_at_patch"].values
@property
def nodes_at_link(self):
return self._mesh["nodes_at_link"].values
@property
def nodes_at_face(self):
return self._mesh["nodes_at_face"].values
@property
def corners_at_face(self):
return self._mesh["corners_at_face"].values
@property
def corners_at_cell(self):
return self._mesh["corners_at_cell"].values
@property
def n_corners_at_cell(self):
return self._mesh["n_corners_at_cell"].values
@property
def cell_at_node(self):
return self._mesh["cell_at_node"].values
[docs]
class VoronoiDelaunayToGraph(VoronoiDelaunay):
[docs]
def __init__(self, xy_of_node, perimeter_links=None):
xy_of_node = require_array(
np.asarray(xy_of_node, dtype=float),
dtype=np.floating,
shape=("n_nodes", 2),
name="xy_of_node",
)
require_length_at_least(xy_of_node, 3, name="xy_of_node")
if perimeter_links is not None:
perimeter_links = require_array(
np.asarray(perimeter_links),
dtype=np.integer,
shape=("n_links", 2),
name="perimeter_links",
)
require_length_at_least(perimeter_links, 3, name="perimeter_links")
super().__init__(xy_of_node)
self._perimeter_links = perimeter_links
mesh = self._mesh
mesh.update(
{
"links_at_patch": xr.DataArray(
self._links_at_patch(
mesh["nodes_at_link"].data, mesh["nodes_at_patch"].data
),
dims=("patch", "Three"),
),
"node_at_cell": xr.DataArray(
reverse_one_to_one(mesh["cell_at_node"].data), dims=("cell",)
),
}
)
mesh.update(
{
"faces_at_cell": xr.DataArray(
self._links_at_patch(
mesh["corners_at_face"].data,
mesh["corners_at_cell"].data,
n_links_at_patch=self.n_corners_at_cell,
),
dims=("cell", "max_faces_per_cell"),
)
}
)
is_extra_link = self._is_extra_link()
self.drop_corners(
self.unbound_corners(is_extra_link), is_extra_link=is_extra_link
)
self.drop_perimeter_faces()
self.drop_perimeter_cells()
self._trim_faces_at_cell()
def _trim_faces_at_cell(self):
"""Remove trailing columns that contain no valid face IDs as -1."""
faces_at_cell = self._mesh["faces_at_cell"]
valid_columns = np.any(faces_at_cell.values != -1, axis=0)
valid_column_ids = np.flatnonzero(valid_columns)
if valid_column_ids.size == 0:
return
width = valid_column_ids[-1] + 1
if width < faces_at_cell.shape[1]:
trimmed = faces_at_cell.values[:, :width].copy()
self._mesh = self._mesh.drop_vars("faces_at_cell")
self._mesh["faces_at_cell"] = xr.DataArray(
trimmed, dims=("cell", "max_faces_per_cell")
)
@staticmethod
def _links_at_patch(nodes_at_link, nodes_at_patch, n_links_at_patch=None):
"""Construct links_at_path from nodes_at_link/nodes_at_path."""
pairs = IntPairMapping(nodes_at_link, values=np.arange(len(nodes_at_link)))
return pairs.get_items(nodes_at_patch, wraparound=True)
[docs]
def is_perimeter_face(self):
return np.any(self.corners_at_face == -1, axis=1)
[docs]
def is_perimeter_cell(self):
"""Check for invalid cells that are on the perimeter.
Sometimes cells on the perimeter of a graph structure are invalid
for landlab's purposes. These are cells with missing corners or
are unbounded.
"""
is_missing_corner = jaggedarray.padded_row_contains(
self.corners_at_cell, self.n_corners_at_cell, value=-1
)
return is_missing_corner | (self.n_corners_at_cell < 3)
[docs]
def is_perimeter_link(self):
if self._perimeter_links is not None:
pairs = IntPairCollection(self._perimeter_links)
is_perimeter_link = pairs.contains_pairs(self.nodes_at_link)
else:
is_perimeter_link = self.is_perimeter_face()
return is_perimeter_link
def _is_extra_link(self):
is_extra_link = np.zeros(self.nodes_at_link.shape[0], dtype=bool)
if self._perimeter_links is not None:
idx, non_perimeter_links = (
self._find_non_perimeter_links_between_perimeter_nodes()
)
outside_links_mask = self._get_outside_links_mask(non_perimeter_links)
idx_outside_links = idx[outside_links_mask]
is_extra_link[idx_outside_links] = True
return is_extra_link
def _find_non_perimeter_links_between_perimeter_nodes(self):
# get perimeter nodes ID
perimeter_nodes = np.unique(self._perimeter_links.ravel())
# sort nodes at links
links_pair = np.sort(self.nodes_at_link, axis=1)
perim_pair = np.sort(self._perimeter_links, axis=1)
# links with tail and head nodes on perimeter nodes
both_ends_on_perim_node = np.isin(links_pair[:, 0], perimeter_nodes) & np.isin(
links_pair[:, 1], perimeter_nodes
)
# identify non perimeter links
dtype = np.dtype([("a", links_pair.dtype), ("b", links_pair.dtype)])
links_view = links_pair.view(dtype).reshape(-1)
perim_view = perim_pair.view(dtype).reshape(-1)
is_perimeter = np.isin(links_view, perim_view)
mask = both_ends_on_perim_node & (~is_perimeter)
idx = np.where(mask)[0]
return idx, links_pair[idx]
def _get_outside_links_mask(self, non_perimeter_links, eps=1e-12):
# get adjacency links for perimeter nodes
edges = np.asarray(self._perimeter_links, dtype=int)
adj = {}
for a, b in edges:
adj.setdefault(a, []).append(b)
adj.setdefault(b, []).append(a)
# check perimeter nodes has only two neighbor links
degrees = {int(n): len(neighbors) for n, neighbors in adj.items()}
bad_nodes = {n: d for n, d in degrees.items() if d != 2}
if bad_nodes:
warnings.warn(
"Invalid perimeter_links:\n"
"Perimeter links must form one closed boundary loop.\n"
"Each boundary node must have exactly two boundary neighbors.\n"
f"Invalid nodes (node -> neighbor count): {bad_nodes}\n"
"Links outside the boundary defined by perimeter_links "
"may not be fully removed.",
UserWarning,
stacklevel=2,
)
return np.zeros(len(non_perimeter_links), dtype=bool)
# order perimeter nodes as a single closed ring
start = next(iter(adj))
ring = [start]
prev, cur = None, start
while True:
n1, n2 = adj[cur]
nxt = n1 if n1 != prev else n2
if nxt == start:
break
ring.append(nxt)
prev, cur = cur, nxt
if len(ring) != len(adj):
warnings.warn(
"Invalid perimeter_links:\n"
"Perimeter links must form one closed boundary loop.\n"
"Perimeter links appear to form multiple disconnected loops.\n"
"Links outside the boundary defined by perimeter_links "
"may not be fully removed.",
UserWarning,
stacklevel=2,
)
return np.zeros(len(non_perimeter_links), dtype=bool)
# check if midpoints of links are within the ring
polygon = Path(np.column_stack([self.x_of_node[ring], self.y_of_node[ring]]))
non_perimeter_links = np.asarray(non_perimeter_links, dtype=int)
mid_points = np.column_stack(
[
0.5
* (
self.x_of_node[non_perimeter_links[:, 0]]
+ self.x_of_node[non_perimeter_links[:, 1]]
),
0.5
* (
self.y_of_node[non_perimeter_links[:, 0]]
+ self.y_of_node[non_perimeter_links[:, 1]]
),
]
)
inside = polygon.contains_points(mid_points, radius=eps)
return ~inside
[docs]
def unbound_corners(self, is_extra_link=None):
if is_extra_link is None:
is_extra_link = self._is_extra_link()
faces_to_drop = np.where(
(self.is_perimeter_face() & ~self.is_perimeter_link()) | is_extra_link
)
unbound_corners = self.corners_at_face[faces_to_drop].reshape((-1,))
return np.unique(unbound_corners[unbound_corners >= 0])
[docs]
def is_bound_corner(self):
corners = np.full(self._mesh.sizes["corner"], True)
corners[self.unbound_corners()] = False
return corners
[docs]
def drop_corners(self, corners, is_extra_link=None):
if len(corners) == 0:
return
# Remove the corners
corners_to_drop = np.asarray(corners, dtype=int)
self.drop_element(corners_to_drop, at="corner")
# Remove bad links
if is_extra_link is None:
is_extra_link = self._is_extra_link()
is_a_link = np.any(self._mesh["corners_at_face"].data != -1, axis=1)
is_a_link[is_extra_link] = False
self.drop_element(np.where(~is_a_link)[0], at="link")
# Remove the bad patches
is_a_patch = np.all(self._mesh["links_at_patch"] >= 0, axis=1)
self.drop_element(np.where(~is_a_patch)[0], at="patch")
[docs]
def drop_perimeter_faces(self):
self.drop_element(np.where(self.is_perimeter_face())[0], at="face")
[docs]
def drop_perimeter_cells(self):
self.drop_element(np.where(self.is_perimeter_cell())[0], at="cell")
[docs]
def ids_with_prefix(self, at):
matches = set()
if at == "patch":
prefix = re.compile(f"^{at}(es)?_at_")
else:
prefix = re.compile(f"^{at}(s)?_at_")
for name in self._mesh.variables:
if prefix.search(name):
matches.add(name)
return matches
[docs]
def ids_with_suffix(self, at):
matches = set()
suffix = re.compile(f"at_{at}$")
for name in self._mesh.variables:
if suffix.search(name):
matches.add(name)
return matches
[docs]
def drop_element(self, ids, at="node"):
dropped_ids = np.asarray(ids, dtype=int)
dropped_ids.sort()
is_a_keeper = np.full(self._mesh.sizes[at], True)
is_a_keeper[dropped_ids] = False
at_ = {}
if at in self._mesh.coords:
x = self._mesh[f"x_of_{at}"].values[is_a_keeper]
y = self._mesh[f"y_of_{at}"].values[is_a_keeper]
data = np.arange(len(x))
at_[at] = xr.DataArray(
data=data,
coords={
f"x_of_{at}": xr.DataArray(x, dims=(at,)),
f"y_of_{at}": xr.DataArray(y, dims=(at,)),
},
dims=(at,),
)
self._mesh = self._mesh.drop_vars([f"x_of_{at}", f"y_of_{at}"])
for name in self.ids_with_suffix(at):
var = self._mesh[name]
at_[name] = xr.DataArray(var.values[is_a_keeper], dims=var.sizes)
self._mesh = self._mesh.drop_vars(list(at_))
self._mesh.update(at_)
for name in self.ids_with_prefix(at):
var = self._mesh[name]
array = var.values.reshape((-1,))
array[np.isin(array, dropped_ids)] = -1
for id_ in dropped_ids[::-1]:
array[array > id_] -= 1
@property
def links_at_patch(self):
return self._mesh["links_at_patch"].values
@property
def node_at_cell(self):
return self._mesh["node_at_cell"].values
@property
def faces_at_cell(self):
return self._mesh["faces_at_cell"].values