Skip to content

Commit e5c56c2

Browse files
AntSimiAntSimi
committed
Add an entry point to track eddy like network
1 parent 19b6fd1 commit e5c56c2

File tree

2 files changed

+245
-1
lines changed

2 files changed

+245
-1
lines changed

‎setup.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,10 @@
2828
],
2929
zip_safe=False,
3030
entry_points=dict(
31-
console_scripts=["MergeEddies = py_eddy_tracker.appli:merge_eddies"]
31+
console_scripts=[
32+
"MergeEddies = py_eddy_tracker.appli:merge_eddies",
33+
"EddyNetworkGroup = py_eddy_tracker.network:build_network",
34+
]
3235
),
3336
package_data={
3437
"py_eddy_tracker.featured_tracking": ["*.nc"],

‎src/py_eddy_tracker/network.py‎

Lines changed: 241 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,241 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
===========================================================================
4+
This file is part of py-eddy-tracker.
5+
6+
py-eddy-tracker is free software: you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation, either version 3 of the License, or
9+
(at your option) any later version.
10+
11+
py-eddy-tracker is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with py-eddy-tracker. If not, see <http://www.gnu.org/licenses/>.
18+
19+
Copyright (c) 2014-2020 by Evan Mason
20+
Email: [email protected]
21+
===========================================================================
22+
"""
23+
24+
import logging
25+
from glob import glob
26+
from netCDF4 import Dataset
27+
from Polygon import Polygon
28+
from numba import njit, uint32
29+
from numpy import empty, array, arange, unique, bincount
30+
from py_eddy_tracker.poly import bbox_intersection, create_vertice
31+
from py_eddy_tracker import EddyParser
32+
33+
34+
logger = logging.getLogger("pet")
35+
36+
NOGROUP = 0
37+
# To be used like a buffer
38+
DATA = dict()
39+
FLIST = list()
40+
41+
42+
def load_contour(filename, xname, yname):
43+
"""
44+
To avoid multiple load of same file we store last 20 result
45+
"""
46+
if filename not in DATA:
47+
with Dataset(filename) as h:
48+
if len(FLIST) >= 20:
49+
DATA.pop(FLIST.pop(0))
50+
FLIST.append(filename)
51+
DATA[filename] = h.variables[xname][:], h.variables[yname][:]
52+
return DATA[filename]
53+
54+
55+
@njit(cache=True)
56+
def get_wrap_vertice(x0, y0, x1, y1, i):
57+
x0_, x1_ = x0[i], x1[i]
58+
if abs(x0_[0] - x1_[0]) > 180:
59+
ref = x0_[0] - x0.dtype.type(180)
60+
x1_ = x1[i]
61+
x1_ = (x1[i] - ref) % 360 + ref
62+
return create_vertice(x0_, y0[i]), create_vertice(x1_, y1[i])
63+
64+
65+
def common_area(x0, y0, x1, y1, minimal_area=False):
66+
nb, _ = x0.shape
67+
cost = empty((nb))
68+
for i in range(nb):
69+
# Get wrapped vertice for index i
70+
v0, v1 = get_wrap_vertice(x0, y0, x1, y1, i)
71+
p0 = Polygon(v0)
72+
p1 = Polygon(v1)
73+
# Area of intersection
74+
intersection = (p0 & p1).area()
75+
# we divide intersection with the little one result from 0 to 1
76+
if minimal_area:
77+
p0_area = p0.area()
78+
p1_area = p1.area()
79+
cost[i] = intersection / min(p0_area, p1_area)
80+
# we divide intersection with polygon merging result from 0 to 1
81+
else:
82+
union = (p0 + p1).area()
83+
cost[i] = intersection / union
84+
return cost
85+
86+
87+
def get_group_array(results, nb_obs):
88+
"""With a loop on all pair of index, we will label each obs with a group
89+
number
90+
"""
91+
nb_obs = array(nb_obs)
92+
day_start = nb_obs.cumsum() - nb_obs
93+
gr = empty(nb_obs.sum(), dtype="u4")
94+
gr[:] = NOGROUP
95+
96+
next_id_group = 1
97+
for i, j, i_ref, i_etu in results:
98+
sl_ref = slice(day_start[i], day_start[i] + nb_obs[i])
99+
sl_etu = slice(day_start[j], day_start[j] + nb_obs[j])
100+
# obs with no groups
101+
m = (gr[sl_ref][i_ref] == NOGROUP) * (gr[sl_etu][i_etu] == NOGROUP)
102+
nb_no_groups = m.sum()
103+
gr[sl_ref][i_ref[m]] = gr[sl_etu][i_etu[m]] = arange(
104+
next_id_group, next_id_group + nb_no_groups
105+
)
106+
next_id_group += nb_no_groups
107+
# associate obs with no group with obs with group
108+
m = (gr[sl_ref][i_ref] != NOGROUP) * (gr[sl_etu][i_etu] == NOGROUP)
109+
gr[sl_etu][i_etu[m]] = gr[sl_ref][i_ref[m]]
110+
m = (gr[sl_ref][i_ref] == NOGROUP) * (gr[sl_etu][i_etu] != NOGROUP)
111+
gr[sl_ref][i_ref[m]] = gr[sl_etu][i_etu[m]]
112+
# case where 2 obs have a different group
113+
m = gr[sl_ref][i_ref] != gr[sl_etu][i_etu]
114+
if m.any():
115+
# Merge of group, ref over etu
116+
for i_, j_ in zip(i_ref[m], i_etu[m]):
117+
g_ref, g_etu = gr[sl_ref][i_], gr[sl_etu][j_]
118+
gr[gr == g_ref] = g_etu
119+
return gr
120+
121+
122+
def save(filenames, gr, out):
123+
"""Function to saved group output
124+
"""
125+
new_i = get_next_index(gr)
126+
nb = gr.shape[0]
127+
dtype = list()
128+
with Dataset(out, "w") as h_out:
129+
with Dataset(filenames[0]) as h_model:
130+
for name, dim in h_model.dimensions.items():
131+
h_out.createDimension(name, len(dim) if name != "obs" else nb)
132+
v = h_out.createVariable(
133+
"track", "u4", ("obs",), True, chunksizes=(min(250000, nb),)
134+
)
135+
d = v[:].copy()
136+
d[new_i] = gr
137+
v[:] = d
138+
for k in h_model.ncattrs():
139+
h_out.setncattr(k, h_model.getncattr(k))
140+
for name, v in h_model.variables.items():
141+
dtype.append(
142+
(
143+
name,
144+
(v.datatype, 50) if "NbSample" in v.dimensions else v.datatype,
145+
)
146+
)
147+
new_v = h_out.createVariable(
148+
name,
149+
v.datatype,
150+
v.dimensions,
151+
True,
152+
chunksizes=(min(25000, nb), 50)
153+
if "NbSample" in v.dimensions
154+
else (min(250000, nb),),
155+
)
156+
for k in v.ncattrs():
157+
if k in ("min", "max",):
158+
continue
159+
new_v.setncattr(k, v.getncattr(k))
160+
data = empty(nb, dtype)
161+
i = 0
162+
debug_active = logger.getEffectiveLevel() == logging.DEBUG
163+
for filename in filenames:
164+
if debug_active:
165+
print(f'Load {filename} to copy', end="\r")
166+
with Dataset(filename) as h_in:
167+
stop = i + len(h_in.dimensions["obs"])
168+
sl = slice(i, stop)
169+
for name, _ in dtype:
170+
v = h_in.variables[name]
171+
v.set_auto_maskandscale(False)
172+
data[name][new_i[sl]] = v[:]
173+
i = stop
174+
if debug_active:
175+
print()
176+
for name, _ in dtype:
177+
v = h_out.variables[name]
178+
v.set_auto_maskandscale(False)
179+
v[:] = data[name]
180+
181+
182+
@njit(cache=True)
183+
def get_next_index(gr):
184+
"""Return for each obs index the new position to join all group
185+
"""
186+
nb_obs_gr = bincount(gr)
187+
i_gr = nb_obs_gr.cumsum() - nb_obs_gr
188+
new_index = empty(gr.shape, dtype=uint32)
189+
for i, g in enumerate(gr):
190+
new_index[i] = i_gr[g]
191+
i_gr[g] += 1
192+
return new_index
193+
194+
195+
def build_network():
196+
parser = EddyParser("Merge eddies")
197+
parser.add_argument(
198+
"identification_regex",
199+
help="Give an expression which will use with glob, currently only netcdf file",
200+
)
201+
parser.add_argument("out", help="output file, currently only netcdf file")
202+
parser.add_argument(
203+
"--window", "-w", type=int, help="Half time window to search eddy", default=1
204+
)
205+
parser.add_argument(
206+
"--intern",
207+
action="store_true",
208+
help="Use intern contour instead of outter contour",
209+
)
210+
args = parser.parse_args()
211+
network(args.identification_regex, args.out, window=args.window, intern=args.intern)
212+
213+
214+
def network(regex, filename_out, window=1, intern=False):
215+
filenames = glob(regex)
216+
filenames.sort()
217+
nb_in = len(filenames)
218+
if intern:
219+
coord = "speed_contour_longitude", "speed_contour_latitude"
220+
else:
221+
coord = "effective_contour_longitude", "effective_contour_latitude"
222+
results, nb_obs = list(), list()
223+
debug_active = logger.getEffectiveLevel() == logging.DEBUG
224+
for i, filename in enumerate(filenames):
225+
if debug_active:
226+
print(f'{filename} compared to {window} next', end="\r")
227+
xi, yi = load_contour(filename, *coord)
228+
nb_obs.append(xi.shape[0])
229+
for j in range(i + 1, min(window + i + 1, nb_in)):
230+
xj, yj = load_contour(filenames[j], *coord)
231+
ii, ij = bbox_intersection(xi, yi, xj, yj)
232+
m = common_area(xi[ii], yi[ii], xj[ij], yj[ij], minimal_area=True) > 0.2
233+
results.append((i, j, ii[m], ij[m]))
234+
if debug_active:
235+
print()
236+
237+
gr = get_group_array(results, nb_obs)
238+
logger.info(
239+
f"{(gr == NOGROUP).sum()} alone / {len(gr)} obs, {len(unique(gr))} groups"
240+
)
241+
save(filenames, gr, filename_out)

0 commit comments

Comments
 (0)