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