diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index 40d138e..27398da 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -8,7 +8,11 @@ ), "The environment variable CGAL_BIN must be set." os.add_dll_directory(os.environ["CGAL_BIN"]) -from oceanmesh.boundary import identify_ocean_boundary_sections +from oceanmesh.boundary import ( + identify_ocean_boundary_sections, + identify_land_boundary_sections, + identify_island_boundary_sections, +) from oceanmesh.clean import ( delete_boundary_faces, delete_exterior_faces, @@ -100,6 +104,8 @@ "Difference", "Intersection", "identify_ocean_boundary_sections", + "identify_land_boundary_sections", + "identify_island_boundary_sections", "Shoreline", "generate_multiscale_mesh", "get_polygon_coordinates", diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py index 8d66784..70295ee 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -1,9 +1,23 @@ -import matplotlib.pyplot as plt +""" +Functions to automatically label boundary segments as either +elevation-specified (i.e., ocean) or no-flux (i.e., land) +based on geometric and topobathymetric aspects. + +Author: Dr. Alexandra Maskell +Date: 2024-01-17 +Edits by: Keith Roberts, 2025-01-28 +""" import numpy as np +from sklearn.neighbors import NearestNeighbors -from .edges import get_winded_boundary_edges +from oceanmesh.edges import get_winded_boundary_edges, get_boundary_edges +import matplotlib.pyplot as plt -__all__ = ["identify_ocean_boundary_sections"] +__all__ = [ + "identify_ocean_boundary_sections", + "identify_land_boundary_sections", + "identify_island_boundary_sections", +] def identify_ocean_boundary_sections( @@ -40,6 +54,7 @@ def identify_ocean_boundary_sections( """ # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) boundary_edges = boundary_edges.flatten() unique_indexes = np.unique(boundary_edges, return_index=True)[1] @@ -104,3 +119,432 @@ def identify_ocean_boundary_sections( for s1, s2 in boundary_sections ] return boundary_sections + + +def identify_land_boundary_sections( + points, + cells, + ocean_boundary, + plot=False, +): + """Identify the contiguous sections on the land boundary based on ocean boundary + + Parameters + ---------- + points: numpy.ndarray + Array of points (x,y) + cells : numpy.ndarray + Array of cells + ocean_boundary : numpy.ndarray + List of tuples of the nodes that define the ocean boundary sections + Note these map back into the points array. + plot : bool, optional + Plot the mesh and the identified boundary sections, by default False + + Returns + -------- + boundary_sections : list + List of tuples of the nodes that define the ocean boundary sections + Note these map back into the points array. + + """ + # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) + boundary_edges = boundary_edges.flatten() + unique_indexes = np.unique(boundary_edges, return_index=True)[1] + boundary_nodes_unmasked = [ + boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + # Define a boolean array of valid nodes + boundary_sections = [] + first = True + for idx, (s1, s2) in enumerate(ocean_boundary): + start_idx = np.where(boundary_nodes_unmasked == s1)[0][0] + end_idx = np.where(boundary_nodes_unmasked == s2)[0][0] + + if (first) & (start_idx > 0): + boundary_sections.append([0, start_idx]) + st_idx = end_idx + first = False + elif idx <= (len(ocean_boundary) - 1): + boundary_sections.append([st_idx, start_idx]) + + if idx == (len(ocean_boundary) - 1): + boundary_sections.append([end_idx, len(boundary_nodes_unmasked) - 1]) + + # Plot the mesh + if plot: + fig, ax = plt.subplots() + ax.triplot(points[:, 0], points[:, 1], cells, color="k", lw=0.1) + + for s1, s2 in boundary_sections: + ax.scatter( + points[boundary_nodes_unmasked[s1:s2], 0], + points[boundary_nodes_unmasked[s1:s2], 1], + 5, + c="g", + ) + ax.set_title("Identified land boundary sections") + plt.show() + + # Map back to the original node indices associated with the points array + boundary_sections = [ + (boundary_nodes_unmasked[s1], boundary_nodes_unmasked[s2]) + for s1, s2 in boundary_sections + ] + + def shift(seq, n): + n = n % len(seq) + return seq[n:] + seq[:n] + + ## Sort land boundaries to occur in cw order + start_lb = np.where(boundary_sections == ocean_boundary[0][1])[0][0] + mainland_boundary_cw = shift(boundary_sections, start_lb) + + return mainland_boundary_cw + + +def identify_island_boundary_sections(points, cells, plot=False): + """Identify the contiguous sections on the island boundary + + Parameters + ---------- + points: numpy.ndarray + Array of points (x,y) + cells : numpy.ndarray + Array of cells + plot : bool, optional + Plot the mesh and the identified boundary sections, by default False + + Returns + -------- + boundary_nodes : list + List of nodes that define the island boundary sections (cw direction) + Note these map back into the points array. + boundary_sections : list + List of tuples of the nodes that define the island boundary sections + Note these map back into the points array. + + """ + # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) + boundary_edges = boundary_edges.flatten() + unique_indexes = np.unique(boundary_edges, return_index=True)[1] + boundary_nodes_unmasked = [ + boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + + all_boundary_edges = get_boundary_edges(cells) + all_boundary_edges = all_boundary_edges.flatten() + unique_indexes = np.unique(all_boundary_edges, return_index=True)[1] + + all_boundary_nodes_unmasked = [ + all_boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + + common_elements = list( + set(boundary_nodes_unmasked).intersection(set(all_boundary_nodes_unmasked)) + ) + all_island_boundary_nodes = boundary_nodes_unmasked + all_boundary_nodes_unmasked + + for item in common_elements: + all_island_boundary_nodes = [ + element for element in all_island_boundary_nodes if element != item + ] + + island_boundary_nodes_winded = [] + choice = 0 + while True: + idx = all_island_boundary_nodes[choice] + islands_boundary_edges = get_winded_boundary_edges(cells, vFirst=idx) + islands_boundary_edges = islands_boundary_edges.flatten() + unique_indexes = np.unique(islands_boundary_edges, return_index=True)[1] + + island_boundary_nodes_unmasked = [ + islands_boundary_edges[unique_index] + for unique_index in sorted(unique_indexes) + ] + island_boundary_nodes_unmasked.reverse() + island_boundary_nodes_winded.append(island_boundary_nodes_unmasked) + + if sum(len(ele) for ele in island_boundary_nodes_winded) == len( + all_island_boundary_nodes + ): + break + + common_elements = list( + set(island_boundary_nodes_unmasked).intersection( + set(all_island_boundary_nodes) + ) + ) + remaining_island_nodes = ( + island_boundary_nodes_unmasked + all_island_boundary_nodes + ) + for item in common_elements: + remaining_island_nodes = [ + element for element in remaining_island_nodes if element != item + ] + choice = np.where(all_island_boundary_nodes == remaining_island_nodes[0])[0][0] + + # Plot the mesh + if plot: + fig, ax = plt.subplots() + ax.triplot(points[:, 0], points[:, 1], cells, color="k", lw=0.1) + for i in range(len(island_boundary_nodes_winded)): + ax.scatter( + points[island_boundary_nodes_winded[i], 0], + points[island_boundary_nodes_winded[i], 1], + 5, + c="r", + ) + ax.set_title("Identified land boundary sections") + plt.show() + + # Map back to the original node indices associated with the points array + island_boundary_sections = [ + (island_boundary_nodes_winded[i][0], island_boundary_nodes_winded[i][-1]) + for i in range(len(island_boundary_nodes_winded)) + ] + + # append first node to end of list: + + for i in range(len(island_boundary_nodes_winded)): + island_boundary_nodes_winded[i] = [ + x + 1 for x in island_boundary_nodes_winded[i] + ] + island_boundary_nodes_winded[i].append(island_boundary_nodes_winded[i][0]) + + return [island_boundary_nodes_winded, island_boundary_sections] + + +def split_list(ll): + index_list = ( + [None] + [i for i in range(1, len(ll)) if ll[i] - ll[i - 1] > 1] + [None] + ) + return [ll[index_list[j - 1] : index_list[j]] for j in range(1, len(index_list))] + + +def identify_boundary_sections_knn( + points, + cells, + shoreline, + edge_length, + plot=False, +): + # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) + boundary_edges = boundary_edges.flatten() + unique_indexes = np.unique(boundary_edges, return_index=True)[1] + boundary_nodes_unmasked = [ + boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + + boundary_points = points[boundary_nodes_unmasked] + + land = shoreline.mainland + land = land[~np.isnan(land[:, 0])] + + knn = NearestNeighbors(n_neighbors=5) + knn.fit(land) + dist_lim = 25 * edge_length.values.min() + ldst, _ = knn.kneighbors(boundary_points, return_distance=True) + ldst = ldst.min(axis=1) + eb_class = ldst > dist_lim + + # count open boundaries + ocean_boundary = [] + # find index of open boundaries + idx_ope = split_list(np.where(eb_class is True)[0]) + for j in range(len(idx_ope)): + if len(idx_ope[j]) > 3: + idx_ope[j] = boundary_nodes_unmasked[idx_ope[j][0] : idx_ope[j][-1]] + ocean_boundary.append((idx_ope[j][0], idx_ope[j][-1])) + # find index of mainland boundaries + mainland_boundary = [] + idx_mland = split_list(np.where(eb_class is False)[0]) + # nmland = len(idx_mland) + for j in range(len(idx_mland)): + idx_mland[j] = boundary_nodes_unmasked[idx_mland[j][0] : idx_mland[j][-1]] + mainland_boundary.append((idx_mland[j][0], idx_mland[j][-1])) + if plot is True: + fig, ax = plt.subplots() + ax.plot( + boundary_points[~eb_class][:, 0], boundary_points[~eb_class][:, 1], "ko" + ) + ax.plot(boundary_points[eb_class][:, 0], boundary_points[eb_class][:, 1], "bo") + + # Define a boolean array of valid nodes + boundary_sections = [] + first = True + for idx, (s1, s2) in enumerate(ocean_boundary): + start_idx = np.where(boundary_nodes_unmasked == s1)[0][0] + end_idx = np.where(boundary_nodes_unmasked == s2)[0][0] + + if (first) & (start_idx > 0): + boundary_sections.append([0, start_idx]) + st_idx = end_idx + first = False + elif idx < (len(ocean_boundary) - 2): + boundary_sections.append([st_idx, start_idx]) + + if idx == (len(ocean_boundary) - 1): + boundary_sections.append([end_idx, len(boundary_nodes_unmasked) - 1]) + + # Plot the mesh + if plot: + fig, ax = plt.subplots() + ax.triplot(points[:, 0], points[:, 1], cells, color="k", lw=0.1) + + for s1, s2 in boundary_sections: + ax.scatter( + points[boundary_nodes_unmasked[s1:s2], 0], + points[boundary_nodes_unmasked[s1:s2], 1], + 5, + c="g", + ) + ax.set_title("Identified land boundary sections") + plt.show() + + # Map back to the original node indices associated with the points array + boundary_sections = [ + (boundary_nodes_unmasked[s1], boundary_nodes_unmasked[s2]) + for s1, s2 in boundary_sections + ] + + def shift(seq, n): + n = n % len(seq) + return seq[n:] + seq[:n] + + ## Sort land boundaries to occur in cw order + # start_lb = np.where(boundary_sections==ocean_boundary[0][1])[0][0] + # mainland_boundary_cw = shift(boundary_sections, start_lb) + + ## check for islands - need to be ccw order (TO ADD) + idx_islands = [] + if len(shoreline.inner) > 0: + all_boundary_edges = get_boundary_edges(cells) + all_boundary_edges = all_boundary_edges.flatten() + unique_indexes = np.unique(all_boundary_edges, return_index=True)[1] + all_boundary_nodes_unmasked = all_boundary_edges[np.sort(unique_indexes)] + all_boundary_points = points[all_boundary_nodes_unmasked] + + islands = shoreline.inner + islands = np.split(islands, np.where(np.isnan(islands[:, 0]) is True)[0]) + islands = [y[~np.isnan(y[:, 0])] for y in islands][:-1] + + for i in islands: + knn = NearestNeighbors(n_neighbors=3) + knn.fit(i) + dist_lim = 5 * edge_length.dx + idst, _ = knn.kneighbors(all_boundary_points, return_distance=True) + idst = idst.min(axis=1) + ebi_class = idst < dist_lim + idx_islands.append( + np.hstack( + [ + all_boundary_nodes_unmasked[np.where(ebi_class is True)[0]], + all_boundary_nodes_unmasked[np.where(ebi_class is True)[0][0]], + ] + ) + ) + if plot is True: + ax.plot( + all_boundary_points[ebi_class][:, 0], + all_boundary_points[ebi_class][:, 1], + "go", + ) + # nislands = len(idx_islands) + + # change index to start at 1 + for i in range(len(idx_islands)): + idx_islands[i] = np.flip(idx_islands[i]) + idx_islands[i] = [x + 1 for x in idx_islands[i]] + + return ocean_boundary, idx_islands # ,ocean boundary,idx_islands + + +def boundary_node_list( + cells, + boundary_sections, + ccw=True, + node_limit=999, +): + """ + Output list of boundary nodes in ccw order + + Parameters + ---------- + cells : numpy.ndarray + Array of cells + boundary_sections : list + List of tuples of the nodes that define the ocean boundary sections + Note these map back into the points array. + ccw : bool, optional + output boundary nodes in ccw order, by default True + node_limit : bool, optional + output boundary nodes in list length of 1000 max, by default 999 + Returns + -------- + boundary_nodes : list of nodes + List of tuples of the nodes that define the ocean boundary sections + Note these map back into the points array. + + """ + # Identify the nodes on the boundary of the mesh + boundary_edges = get_winded_boundary_edges(cells) + boundary_edges = boundary_edges.flatten() + unique_indexes = np.unique(boundary_edges, return_index=True)[1] + boundary_nodes_unmasked = [ + boundary_edges[unique_index] for unique_index in sorted(unique_indexes) + ] + + boundary_sections_limit = [] + for idx, (s1, s2) in enumerate(boundary_sections): + start_idx = np.where(boundary_nodes_unmasked == s1)[0][0] + end_idx = np.where(boundary_nodes_unmasked == s2)[0][0] + + if (end_idx - start_idx) < node_limit: + boundary_sections_limit.append( + (boundary_nodes_unmasked[start_idx], boundary_nodes_unmasked[end_idx]) + ) + else: + n = int(np.ceil((end_idx - start_idx) / node_limit)) + splt_idx = start_idx.copy() + for i in range(0, n): + if i < (n - 1): + boundary_sections_limit.append( + ( + boundary_nodes_unmasked[splt_idx], + boundary_nodes_unmasked[splt_idx + node_limit], + ) + ) + elif i == (n - 1): + boundary_sections_limit.append( + ( + boundary_nodes_unmasked[splt_idx], + boundary_nodes_unmasked[end_idx], + ) + ) + splt_idx = splt_idx + node_limit + + boundary_nodes = [] + for idx, (s1, s2) in enumerate(boundary_sections_limit): + start_idx = np.where(boundary_nodes_unmasked == s1)[0][0] + end_idx = np.where(boundary_nodes_unmasked == s2)[0][0] + if end_idx < (len(boundary_nodes_unmasked) - 1): + list_nodes = boundary_nodes_unmasked[start_idx : end_idx + 1] + elif end_idx == (len(boundary_nodes_unmasked) - 1): + list_nodes = boundary_nodes_unmasked[start_idx:end_idx] + list_nodes.append(boundary_nodes_unmasked[0]) + list_nodes = [x + 1 for x in list_nodes] + + # if (land==True) & (idx>0): + # list_nodes.append(boundary_nodes[idx-1][0]) + + if ccw is False: + # flip node direction + list_nodes.reverse() + + boundary_nodes.append(list_nodes) + + return boundary_nodes diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py index 00af47d..0d89ba7 100644 --- a/oceanmesh/geodata.py +++ b/oceanmesh/geodata.py @@ -794,9 +794,9 @@ def __init__(self, dem, crs="EPSG:4326", bbox=None, extrapolate=False): topobathy = np.transpose(topobathy, (1, 0)) # Ensure its a floating point array topobathy = topobathy.astype(np.float64) - topobathy[topobathy == nodata_value] = ( - np.nan - ) # set the no-data value to nan + topobathy[ + topobathy == nodata_value + ] = np.nan # set the no-data value to nan elif not dem.exists(): raise FileNotFoundError(f"File {dem} could not be located.") diff --git a/setup.cfg b/setup.cfg index 6d2d3ce..d57bd81 100644 --- a/setup.cfg +++ b/setup.cfg @@ -12,9 +12,9 @@ classifiers = Development Status :: 4 - Beta License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+) Operating System :: OS Independent - Programming Language :: Python :: 3.7 Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 + Programming Language :: Python :: 3.10 Topic :: Scientific/Engineering Topic :: Scientific/Engineering :: Mathematics Topic :: Scientific/Engineering :: Physics @@ -35,6 +35,7 @@ install_requires = rioxarray scikit-fmm scikit-image + scikit-learn python_requires = >=3.0 [versioneer]