Skip to content

Commit 6921457

Browse files
Merge pull request AntSimi#96 from ludwigVonKoopa/validation_particle
update coherence forward/backward
2 parents bfd6e6a + 2f20ca6 commit 6921457

File tree

8 files changed

+199
-155
lines changed

8 files changed

+199
-155
lines changed

CHANGELOG.rst

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,9 +12,12 @@ Changed
1212
^^^^^^^
1313
Fixed
1414
^^^^^
15+
- GridCollection get_next_time_step & get_previous_time_step needed more files to work in the dataset list.
16+
The loop needed explicitly self.dataset[i+-1] even when i==0, therefore indice went out of range
1517
Added
1618
^^^^^
1719

20+
1821
[3.4.0] - 2021-03-29
1922
--------------------
2023
Changed

examples/16_network/pet_follow_particle.py

Lines changed: 2 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,6 @@
1616
from py_eddy_tracker.dataset.grid import GridCollection
1717
from py_eddy_tracker.observations.groups import particle_candidate
1818
from py_eddy_tracker.observations.network import NetworkObservations
19-
from py_eddy_tracker.poly import group_obs
2019

2120
start_logger().setLevel("ERROR")
2221

@@ -128,25 +127,19 @@ def update(frame):
128127
# ^^^^^^^^^^^^^^^^^^
129128
step = 1 / 60.0
130129

131-
x, y = meshgrid(arange(24, 36, step), arange(31, 36, step))
132-
x0, y0 = x.reshape(-1), y.reshape(-1)
133-
# Pre-order to speed up
134-
_, i = group_obs(x0, y0, 1, 360)
135-
x0, y0 = x0[i], y0[i]
136-
137130
t_start, t_end = n.period
138131
dt = 14
139132

140133
shape = (n.obs.size, 2)
141134
# Forward run
142135
i_target_f, pct_target_f = -ones(shape, dtype="i4"), zeros(shape, dtype="i1")
143136
for t in range(t_start, t_end - dt):
144-
particle_candidate(x0, y0, c, n, t, i_target_f, pct_target_f, n_days=dt)
137+
particle_candidate(c, n, step, t, i_target_f, pct_target_f, n_days=dt)
145138

146139
# Backward run
147140
i_target_b, pct_target_b = -ones(shape, dtype="i4"), zeros(shape, dtype="i1")
148141
for t in range(t_start + dt, t_end):
149-
particle_candidate(x0, y0, c, n, t, i_target_b, pct_target_b, n_days=-dt)
142+
particle_candidate(c, n, step, t, i_target_b, pct_target_b, n_days=-dt)
150143

151144
# %%
152145
fig = plt.figure(figsize=(10, 10))

requirements.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
matplotlib
22
netCDF4
3-
numba
3+
numba>=0.53
44
numpy
55
opencv-python
66
pint

src/py_eddy_tracker/__init__.py

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -32,9 +32,7 @@
3232

3333

3434
def start_logger():
35-
FORMAT_LOG = (
36-
"%(levelname)-8s %(asctime)s %(module)s.%(funcName)s :\n\t\t\t\t\t%(message)s"
37-
)
35+
FORMAT_LOG = "%(levelname)-8s %(asctime)s %(module)s.%(funcName)s :\n\t%(message)s"
3836
logger = logging.getLogger("pet")
3937
if len(logger.handlers) == 0:
4038
# set up logging to CONSOLE

src/py_eddy_tracker/dataset/grid.py

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -2260,11 +2260,13 @@ def from_netcdf_cube(cls, filename, x_name, y_name, t_name, heigth=None):
22602260
@classmethod
22612261
def from_netcdf_list(cls, filenames, t, x_name, y_name, indexs=None, heigth=None):
22622262
new = cls()
2263-
for i, t in enumerate(t):
2264-
d = RegularGridDataset(filenames[i], x_name, y_name, indexs=indexs)
2263+
for i, _t in enumerate(t):
2264+
filename = filenames[i]
2265+
logger.debug(f"load file {i:02d}/{len(t)} t={_t} : {filename}")
2266+
d = RegularGridDataset(filename, x_name, y_name, indexs=indexs)
22652267
if heigth is not None:
22662268
d.add_uv(heigth)
2267-
new.datasets.append((t, d))
2269+
new.datasets.append((_t, d))
22682270
return new
22692271

22702272
def shift_files(self, t, filename, heigth=None, **rgd_kwargs):
@@ -2276,6 +2278,7 @@ def shift_files(self, t, filename, heigth=None, **rgd_kwargs):
22762278
if heigth is not None:
22772279
d.add_uv(heigth)
22782280
self.datasets.append((t, d))
2281+
logger.debug(f"shift and adding i={len(self.datasets)} t={t} : {filename}")
22792282

22802283
def interp(self, grid_name, t, lons, lats, method="bilinear"):
22812284
"""
@@ -2436,6 +2439,7 @@ def advect(
24362439
else:
24372440
mask_particule += isnan(x) + isnan(y)
24382441
while True:
2442+
logger.debug(f"advect : t={t}")
24392443
if (backward and t <= t1) or (not backward and t >= t1):
24402444
t0, u0, v0, m0 = t1, u1, v1, m1
24412445
t1, d1 = generator.__next__()
@@ -2462,25 +2466,21 @@ def advect(
24622466
yield t, x, y
24632467

24642468
def get_next_time_step(self, t_init):
2465-
first = True
24662469
for i, (t, dataset) in enumerate(self.datasets):
24672470
if t < t_init:
24682471
continue
2469-
if first:
2470-
first = False
2471-
yield self.datasets[i - 1]
2472+
2473+
logger.debug(f"i={i}, t={t}, dataset={dataset}")
24722474
yield t, dataset
24732475

24742476
def get_previous_time_step(self, t_init):
2475-
first = True
24762477
i = len(self.datasets)
24772478
for t, dataset in reversed(self.datasets):
24782479
i -= 1
24792480
if t > t_init:
24802481
continue
2481-
if first:
2482-
first = False
2483-
yield self.datasets[i + 1]
2482+
2483+
logger.debug(f"i={i}, t={t}, dataset={dataset}")
24842484
yield t, dataset
24852485

24862486

src/py_eddy_tracker/observations/groups.py

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -87,11 +87,9 @@ def advect(x, y, c, t0, n_days):
8787
return t, x, y
8888

8989

90-
def particle_candidate(x, y, c, eddies, t_start, i_target, pct, **kwargs):
90+
def particle_candidate(c, eddies, step_mesh, t_start, i_target, pct, **kwargs):
9191
"""Select particles within eddies, advect them, return target observation and associated percentages
9292
93-
:param np.array(float) x: longitude of particles
94-
:param np.array(float) y: latitude of particles
9593
:param `~py_eddy_tracker.dataset.grid.GridCollection` c: GridCollection with speed for particles
9694
:param GroupEddiesObservations eddies: GroupEddiesObservations considered
9795
:param int t_start: julian day of the advection
@@ -102,24 +100,26 @@ def particle_candidate(x, y, c, eddies, t_start, i_target, pct, **kwargs):
102100
"""
103101
# Obs from initial time
104102
m_start = eddies.time == t_start
105-
106103
e = eddies.extract_with_mask(m_start)
104+
107105
# to be able to get global index
108106
translate_start = where(m_start)[0]
109-
# Identify particle in eddies (only in core)
110-
i_start = e.contains(x, y, intern=True)
111-
m = i_start != -1
112107

113-
x, y, i_start = x[m], y[m], i_start[m]
114-
# Advect
108+
x, y, i_start = e.create_particles(step_mesh)
109+
110+
# Advection
115111
t_end, x, y = advect(x, y, c, t_start, **kwargs)
112+
116113
# eddies at last date
117114
m_end = eddies.time == t_end / 86400
118115
e_end = eddies.extract_with_mask(m_end)
116+
119117
# to be able to get global index
120118
translate_end = where(m_end)[0]
119+
121120
# Id eddies for each alive particle (in core and extern)
122121
i_end = e_end.contains(x, y)
122+
123123
# compute matrix and fill target array
124124
get_matrix(i_start, i_end, translate_start, translate_end, i_target, pct)
125125

0 commit comments

Comments
 (0)