From c32a9597b53cfa14676fdab08a5b862107415145 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Wed, 17 Jan 2024 14:19:36 -0300 Subject: [PATCH 1/5] IMPRV: adding in functions to label boundary conditions. Credit to Dr. Alexandra Maskell --- oceanmesh/boundary.py | 415 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 413 insertions(+), 2 deletions(-) diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py index 8d66784..87d83f0 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -1,7 +1,16 @@ -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 +''' 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"] @@ -40,6 +49,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 +114,404 @@ 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, + ccw=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 + 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(l): + index_list = [None] + [i for i in range(1, len(l)) if l[i] - l[i - 1] > 1] + [None] + return [l[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] + + x = points[boundary_edges,0] + y = points[boundary_edges,1] + # x_mid = x.mean(axis=0) + # y_mid = y.mean(axis=0) + # eb_mid = np.row_stack((x_mid, y_mid)).T + + 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 = [] + nope = (eb_class[:-1] < eb_class[1:]).sum() + # find index of open boundaries + idx_ope = split_list(np.where(eb_class==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==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 == 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])==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==True)[0]],all_boundary_nodes_unmasked[np.where(ebi_class==True)[0][0]]])) + if plot == 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, + land=False, + ccw=True, + node_limit = 999, +): + """ output list of boundary nodes in c + + 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==False: + #flip node direction + list_nodes.reverse() + + boundary_nodes.append(list_nodes) + + return boundary_nodes From fad991782cb5979463c47a9bc8ab980a75646c4b Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Wed, 17 Jan 2024 14:32:14 -0300 Subject: [PATCH 2/5] IMPRV: fix flake8 issues --- oceanmesh/boundary.py | 325 +++++++++++++++++++++++------------------- 1 file changed, 176 insertions(+), 149 deletions(-) diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py index 87d83f0..24a78c7 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -1,11 +1,10 @@ -''' -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. +"""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 -''' +""" import numpy as np from sklearn.neighbors import NearestNeighbors @@ -152,27 +151,26 @@ def identify_land_boundary_sections( ] # 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]) + 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]) - + 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], @@ -182,30 +180,25 @@ def identify_land_boundary_sections( ) 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] + 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, - ccw=False -): + +def identify_island_boundary_sections(points, cells, plot=False, ccw=False): """Identify the contiguous sections on the land boundary based on ocean boundary Parameters @@ -234,44 +227,59 @@ def identify_island_boundary_sections( 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))) + + 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] + all_island_boundary_nodes = [ + element for element in all_island_boundary_nodes if element != item + ] island_boundary_nodes_winded = [] choice = 0 - while True: + 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) + 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): + + 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 + + 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] - + 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() @@ -285,24 +293,30 @@ def identify_island_boundary_sections( ) 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: - + + # 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]) - + 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(l): - index_list = [None] + [i for i in range(1, len(l)) if l[i] - l[i - 1] > 1] + [None] - return [l[index_list[j - 1]:index_list[j]] for j in range(1, len(index_list))] + +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, @@ -311,7 +325,6 @@ def identify_boundary_sections_knn( 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() @@ -319,69 +332,63 @@ def identify_boundary_sections_knn( boundary_nodes_unmasked = [ boundary_edges[unique_index] for unique_index in sorted(unique_indexes) ] - + boundary_points = points[boundary_nodes_unmasked] - - x = points[boundary_edges,0] - y = points[boundary_edges,1] - # x_mid = x.mean(axis=0) - # y_mid = y.mean(axis=0) - # eb_mid = np.row_stack((x_mid, y_mid)).T land = shoreline.mainland - land = land[~np.isnan(land[:,0])] - + 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) + 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 = [] - nope = (eb_class[:-1] < eb_class[1:]).sum() # find index of open boundaries - idx_ope = split_list(np.where(eb_class==True)[0]) + 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])) + 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==False)[0]) - nmland = len(idx_mland) + 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 == True: + 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') - + 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]) + 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]) - + 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], @@ -391,14 +398,13 @@ def identify_boundary_sections_knn( ) 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] @@ -407,10 +413,9 @@ def shift(seq, n): # 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) + ## check for islands - need to be ccw order (TO ADD) idx_islands = [] - if len(shoreline.inner) >0: + 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] @@ -418,28 +423,38 @@ def shift(seq, n): all_boundary_points = points[all_boundary_nodes_unmasked] islands = shoreline.inner - islands = np.split(islands,np.where(np.isnan(islands[:,0])==True)[0]) - islands = [y[~np.isnan(y[:,0])] for y in islands][:-1] - + 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) + 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==True)[0]],all_boundary_nodes_unmasked[np.where(ebi_class==True)[0][0]]])) - if plot == True: - ax.plot(all_boundary_points[ebi_class][:,0],all_boundary_points[ebi_class ][:,1],'go') - nislands = len(idx_islands) - + 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 + idx_islands[i] = [x + 1 for x in idx_islands[i]] + + return ocean_boundary, idx_islands # ,ocean boundary,idx_islands def boundary_node_list( @@ -447,9 +462,9 @@ def boundary_node_list( boundary_sections, land=False, ccw=True, - node_limit = 999, + node_limit=999, ): - """ output list of boundary nodes in c + """output list of boundary nodes in c Parameters ---------- @@ -475,43 +490,55 @@ def boundary_node_list( 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])) + ] + + 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)) + 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 - + 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): + 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==False: - #flip node direction + 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 From 6907d162f2b5107ccba8060b487e9f512b260954 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Tue, 28 Jan 2025 20:44:00 -0500 Subject: [PATCH 3/5] Improving auto boundary functions --- oceanmesh/__init__.py | 4 +++- oceanmesh/boundary.py | 16 ++++++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index 40d138e..e049b6d 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -8,7 +8,7 @@ ), "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 +100,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 24a78c7..0b77fb6 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -1,9 +1,11 @@ -"""Functions to automatically label boundary segments as either +""" +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 @@ -11,7 +13,9 @@ 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( @@ -198,8 +202,8 @@ def shift(seq, n): return mainland_boundary_cw -def identify_island_boundary_sections(points, cells, plot=False, ccw=False): - """Identify the contiguous sections on the land boundary based on ocean boundary +def identify_island_boundary_sections(points, cells, plot=False): + """Identify the contiguous sections on the island boundary Parameters ---------- @@ -460,11 +464,11 @@ def shift(seq, n): def boundary_node_list( cells, boundary_sections, - land=False, ccw=True, node_limit=999, ): - """output list of boundary nodes in c + """ + Output list of boundary nodes in ccw order Parameters ---------- From 3650544ce3dd3d56d52faf880ff36b33a9c9870c Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Tue, 28 Jan 2025 20:47:24 -0500 Subject: [PATCH 4/5] Linting --- oceanmesh/__init__.py | 6 +++++- oceanmesh/boundary.py | 10 ++++++---- oceanmesh/geodata.py | 6 +++--- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index e049b6d..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, identify_land_boundary_sections, identify_island_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, diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py index 0b77fb6..70295ee 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -13,9 +13,11 @@ from oceanmesh.edges import get_winded_boundary_edges, get_boundary_edges import matplotlib.pyplot as plt -__all__ = ["identify_ocean_boundary_sections", - "identify_land_boundary_sections", - "identify_island_boundary_sections"] +__all__ = [ + "identify_ocean_boundary_sections", + "identify_land_boundary_sections", + "identify_island_boundary_sections", +] def identify_ocean_boundary_sections( @@ -203,7 +205,7 @@ def shift(seq, n): def identify_island_boundary_sections(points, cells, plot=False): - """Identify the contiguous sections on the island boundary + """Identify the contiguous sections on the island boundary Parameters ---------- 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.") From 415703d5e9608ba611d15c2642c7cb92a4f1590d Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Tue, 28 Jan 2025 20:50:23 -0500 Subject: [PATCH 5/5] Adding scikit-learn --- setup.cfg | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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]