Skip to content

Commit 6c05fd7

Browse files
bjorkqvibjorkqvi
committed
Added function to calculate coordinate for the SWAN SEGMENT command based on the boundary mask of a grid.
1 parent 2ebda09 commit 6c05fd7

File tree

2 files changed

+304
-0
lines changed

2 files changed

+304
-0
lines changed

‎dnora/aux_funcs.py

+118
Original file line numberDiff line numberDiff line change
@@ -504,3 +504,121 @@ def determine_patch_periods(times, start_time, end_time):
504504
inds[0:len(block)] = []
505505

506506
return patch_start, patch_end
507+
508+
def identify_boundary_edges(boundary_mask: np.ndarray) -> list[str]:
509+
"""Identifies which edges has some boundary points
510+
511+
North = [-1,:]
512+
South = [0,:]
513+
East = [:,-1]
514+
West = [:,0]
515+
"""
516+
edges = []
517+
518+
if boundary_mask.shape[1] > 2:
519+
n0 = 1
520+
n1 = -1
521+
else:
522+
n0 = 0
523+
n1 = None
524+
525+
if np.any(boundary_mask[-1,n0:n1]):
526+
edges.append('N')
527+
528+
if np.any(boundary_mask[0,n0:n1]):
529+
edges.append('S')
530+
531+
if boundary_mask.shape[0] > 2:
532+
n0 = 1
533+
n1 = -1
534+
else:
535+
n0 = 0
536+
n1 = None
537+
538+
if np.any(boundary_mask[n0:n1,-1]):
539+
edges.append('E')
540+
541+
if np.any(boundary_mask[n0:n1,0]):
542+
edges.append('W')
543+
544+
return edges
545+
546+
def create_ordered_boundary_list(edge_list):
547+
"""Gets all edges, but given in a continuous clockwise direction.
548+
If this is not possible (e.g. ['N','S']), then ampty list is returned."""
549+
full_list = ['N', 'E', 'S', 'W']
550+
for ind, edge in enumerate(full_list):
551+
if edge not in edge_list:
552+
full_list[ind] = ''
553+
full_array = np.array(full_list)
554+
555+
ct = 0
556+
while (np.where(full_array == '')[0][-1] != len(np.where(full_array == '')[0])-1) and ct < 5:
557+
full_array = np.roll(full_array,1)
558+
ct += 1
559+
560+
if ct > 4:
561+
print(f'No continuous boundary can be found for edges {edge_list}. Returning empy list.')
562+
return []
563+
564+
full_array = full_array[full_array != '']
565+
566+
return full_array.tolist()
567+
568+
def get_coords_for_boundary_edges(edges: list, lon_edges: tuple[float, float], lat_edges: tuple[float, float]) -> tuple[np.ndarray, np.ndarray]:
569+
"""Create coordinate vectors for clockwise running edges.
570+
571+
Assumes that edges are clockwise and continuous, which is imposed by the
572+
function create_ordered_boundary_list.
573+
574+
Empty list return empty arrays.
575+
"""
576+
lon = []
577+
lat = []
578+
for edge in edges:
579+
if edge == 'N':
580+
lon.append(lon_edges[0])
581+
lat.append(lat_edges[1])
582+
if edge == 'S':
583+
lon.append(lon_edges[1])
584+
lat.append(lat_edges[0])
585+
if edge == 'W':
586+
lon.append(lon_edges[0])
587+
lat.append(lat_edges[0])
588+
if edge == 'E':
589+
lon.append(lon_edges[1])
590+
lat.append(lat_edges[1])
591+
592+
if edges:
593+
edge = edges[-1] # Close the loop
594+
if edge == 'N':
595+
lon.append(lon_edges[1])
596+
lat.append(lat_edges[1])
597+
if edge == 'S':
598+
lon.append(lon_edges[0])
599+
lat.append(lat_edges[0])
600+
if edge == 'W':
601+
lon.append(lon_edges[0])
602+
lat.append(lat_edges[1])
603+
if edge == 'E':
604+
lon.append(lon_edges[1])
605+
lat.append(lat_edges[0])
606+
607+
608+
return np.array(lon), np.array(lat)
609+
610+
611+
def create_swan_segment_coords(boundary_mask, lon_edges, lat_edges):
612+
"""Createsa longitude and latitude arrays for the SWAN BOUND SEGEMENT
613+
command based on boundary mask.
614+
615+
Identifies edges (north, south etc.) and sets all the points on the edges
616+
as boundary points.
617+
618+
If no continuous boundary can be identified, it returns empty list.
619+
"""
620+
621+
edge_list = identify_boundary_edges(boundary_mask)
622+
clean_edge_list = create_ordered_boundary_list(edge_list)
623+
lon, lat = get_coords_for_boundary_edges(clean_edge_list, lon_edges, lat_edges)
624+
return lon, lat

‎tests/unittests/utest_boundary_points.py

+186
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
import unittest
2+
from dnora.grd import Grid
3+
from dnora.grd.boundary import EdgesAsBoundary
4+
from dnora.aux_funcs import identify_boundary_edges, create_ordered_boundary_list, get_coords_for_boundary_edges, create_swan_segment_coords
5+
import numpy as np
6+
from copy import copy
7+
8+
class IdentifyBnd(unittest.TestCase):
9+
def test_all_edges(self):
10+
grid = Grid(lon=(0,1), lat=(0,1))
11+
grid.set_spacing(nx=4, ny =4)
12+
grid.set_boundary(EdgesAsBoundary(edges=['N','W','S','E']))
13+
14+
self.assertEqual(set(identify_boundary_edges(grid.boundary_mask())),set(['N','W','S','E']))
15+
16+
def test_one_missing(self):
17+
grid = Grid(lon=(0,1), lat=(0,1))
18+
grid.set_spacing(nx=4, ny =5)
19+
20+
full_list = ['N','W','S','E']
21+
22+
for not_this in full_list:
23+
incomplete_list = copy(full_list)
24+
incomplete_list.remove(not_this)
25+
26+
grid.set_boundary(EdgesAsBoundary(edges=incomplete_list))
27+
self.assertEqual(set(identify_boundary_edges(grid.boundary_mask())),set(incomplete_list))
28+
29+
class OrderBnd(unittest.TestCase):
30+
def test_order_full_list(self):
31+
edges = ['S','E']
32+
self.assertEqual(create_ordered_boundary_list(edges), ['E','S'])
33+
34+
edges = ['N','E']
35+
self.assertEqual(create_ordered_boundary_list(edges), ['N','E'])
36+
37+
edges = ['W','E','S']
38+
self.assertEqual(create_ordered_boundary_list(edges), ['E','S','W'])
39+
40+
edges = ['W','E','N']
41+
self.assertEqual(create_ordered_boundary_list(edges), ['W','N','E'])
42+
43+
edges = ['W','E']
44+
self.assertEqual(create_ordered_boundary_list(edges), [])
45+
46+
47+
class BoundaryCoordinates(unittest.TestCase):
48+
def test_single_edge(self):
49+
grid = Grid(lon=(0,1), lat=(2,3))
50+
grid.set_spacing(nx=4, ny =5)
51+
52+
edges = ['N']
53+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
54+
self.assertEqual(np.array_equal(lons,np.array([0,1])), True)
55+
self.assertEqual(np.array_equal(lats,np.array([3,3])), True)
56+
57+
edges = ['S']
58+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
59+
self.assertEqual(np.array_equal(lons,np.array([1,0])), True)
60+
self.assertEqual(np.array_equal(lats,np.array([2,2])), True)
61+
62+
edges = ['W']
63+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
64+
self.assertEqual(np.array_equal(lons,np.array([0,0])), True)
65+
self.assertEqual(np.array_equal(lats,np.array([2,3])), True)
66+
67+
edges = ['E']
68+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
69+
self.assertEqual(np.array_equal(lons,np.array([1,1])), True)
70+
self.assertEqual(np.array_equal(lats,np.array([3,2])), True)
71+
72+
73+
def test_two_edge(self):
74+
grid = Grid(lon=(0,1), lat=(2,3))
75+
grid.set_spacing(nx=4, ny =5)
76+
77+
edges = ['N','E']
78+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
79+
self.assertEqual(np.array_equal(lons,np.array([0,1,1])), True)
80+
self.assertEqual(np.array_equal(lats,np.array([3,3,2])), True)
81+
82+
edges = ['S','W']
83+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
84+
self.assertEqual(np.array_equal(lons,np.array([1,0,0])), True)
85+
self.assertEqual(np.array_equal(lats,np.array([2,2,3])), True)
86+
87+
edges = ['W','N']
88+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
89+
self.assertEqual(np.array_equal(lons,np.array([0,0,1])), True)
90+
self.assertEqual(np.array_equal(lats,np.array([2,3,3])), True)
91+
92+
edges = ['E','S']
93+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
94+
self.assertEqual(np.array_equal(lons,np.array([1,1,0])), True)
95+
self.assertEqual(np.array_equal(lats,np.array([3,2,2])), True)
96+
97+
def test_three_edge(self):
98+
grid = Grid(lon=(0,1), lat=(2,3))
99+
grid.set_spacing(nx=4, ny =5)
100+
101+
edges = ['N','E','S']
102+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
103+
self.assertEqual(np.array_equal(lons,np.array([0,1,1,0])), True)
104+
self.assertEqual(np.array_equal(lats,np.array([3,3,2,2])), True)
105+
106+
edges = ['S','W','N']
107+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
108+
self.assertEqual(np.array_equal(lons,np.array([1,0,0,1])), True)
109+
self.assertEqual(np.array_equal(lats,np.array([2,2,3,3])), True)
110+
111+
edges = ['W','N','E']
112+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
113+
self.assertEqual(np.array_equal(lons,np.array([0,0,1,1])), True)
114+
self.assertEqual(np.array_equal(lats,np.array([2,3,3,2])), True)
115+
116+
edges = ['E','S','W']
117+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
118+
self.assertEqual(np.array_equal(lons,np.array([1,1,0,0])), True)
119+
self.assertEqual(np.array_equal(lats,np.array([3,2,2,3])), True)
120+
121+
122+
def test_four_edge(self):
123+
grid = Grid(lon=(0,1), lat=(2,3))
124+
grid.set_spacing(nx=4, ny =5)
125+
126+
edges = ['N','E','S','W']
127+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
128+
self.assertEqual(np.array_equal(lons,np.array([0,1,1,0,0])), True)
129+
self.assertEqual(np.array_equal(lats,np.array([3,3,2,2,3])), True)
130+
131+
def test_invalid_edge(self):
132+
grid = Grid(lon=(0,1), lat=(2,3))
133+
grid.set_spacing(nx=4, ny =5)
134+
135+
edges = []
136+
lons, lats = get_coords_for_boundary_edges(edges, grid.lon_edges(), grid.lat_edges())
137+
self.assertEqual(np.array_equal(lons,np.array([])), True)
138+
self.assertEqual(np.array_equal(lats,np.array([])), True)
139+
140+
class FullPipeline(unittest.TestCase):
141+
def test_single_edge(self):
142+
grid = Grid(lon=(0,1), lat=(2,3))
143+
grid.set_spacing(nx=4, ny =5)
144+
145+
grid.set_boundary(EdgesAsBoundary(edges=['N']))
146+
lons, lats = create_swan_segment_coords(grid.boundary_mask(), grid.lon_edges(), grid.lat_edges())
147+
self.assertEqual(np.array_equal(lons,np.array([0,1])), True)
148+
self.assertEqual(np.array_equal(lats,np.array([3,3])), True)
149+
150+
def test_two_edges(self):
151+
grid = Grid(lon=(0,1), lat=(2,3))
152+
grid.set_spacing(nx=4, ny =5)
153+
154+
grid.set_boundary(EdgesAsBoundary(edges=['N','W']))
155+
lons, lats = create_swan_segment_coords(grid.boundary_mask(), grid.lon_edges(), grid.lat_edges())
156+
self.assertEqual(np.array_equal(lons,np.array([0,0,1])), True)
157+
self.assertEqual(np.array_equal(lats,np.array([2,3,3])), True)
158+
159+
def test_three_edges(self):
160+
grid = Grid(lon=(0,1), lat=(2,3))
161+
grid.set_spacing(nx=4, ny =5)
162+
163+
grid.set_boundary(EdgesAsBoundary(edges=['N','W','E']))
164+
lons, lats = create_swan_segment_coords(grid.boundary_mask(), grid.lon_edges(), grid.lat_edges())
165+
self.assertEqual(np.array_equal(lons,np.array([0,0,1,1])), True)
166+
self.assertEqual(np.array_equal(lats,np.array([2,3,3,2])), True)
167+
168+
def test_no_edges(self):
169+
grid = Grid(lon=(0,1), lat=(2,3))
170+
grid.set_spacing(nx=4, ny =5)
171+
172+
lons, lats = create_swan_segment_coords(grid.boundary_mask(), grid.lon_edges(), grid.lat_edges())
173+
self.assertEqual(np.array_equal(lons,np.array([])), True)
174+
self.assertEqual(np.array_equal(lats,np.array([])), True)
175+
176+
def test_invalid_edges(self):
177+
grid = Grid(lon=(0,1), lat=(2,3))
178+
grid.set_spacing(nx=4, ny =5)
179+
grid.set_boundary(EdgesAsBoundary(edges=['W','E']))
180+
181+
lons, lats = create_swan_segment_coords(grid.boundary_mask(), grid.lon_edges(), grid.lat_edges())
182+
self.assertEqual(np.array_equal(lons,np.array([])), True)
183+
self.assertEqual(np.array_equal(lats,np.array([])), True)
184+
185+
if __name__ == '__main__':
186+
unittest.main()

0 commit comments

Comments
 (0)