Skip to content

Commit 6f5ff2b

Browse files
committed
Solution for issue AntSimi#11, wrap longitude if asked
1 parent ac0188b commit 6f5ff2b

File tree

3 files changed

+62
-17
lines changed

3 files changed

+62
-17
lines changed

src/py_eddy_tracker/generic.py

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
linspace,
1515
interp,
1616
where,
17+
isnan,
1718
)
1819
from numba import njit, prange
1920
from numpy.linalg import lstsq
@@ -287,3 +288,50 @@ def split_line(x, y, i):
287288
new_y[new_j] = nan
288289
new_j += 1
289290
return new_x, new_y
291+
292+
293+
@njit(cache=True)
294+
def wrap_longitude(x, y, ref, cut=False):
295+
if cut:
296+
indexs = list()
297+
nb = x.shape[0]
298+
new_previous = (x[0] - ref) % 360
299+
x_previous = x[0]
300+
for i in range(1, nb):
301+
x_ = x[i]
302+
new_x = (x_ - ref) % 360
303+
if not isnan(x_) and not isnan(x_previous):
304+
d_new = new_x - new_previous
305+
d = x_ - x_previous
306+
if d != d_new:
307+
indexs.append(i)
308+
x_previous, new_previous = x_, new_x
309+
310+
nb_indexs = len(indexs)
311+
new_size = nb + nb_indexs * 3
312+
out_x = empty(new_size, dtype=x.dtype)
313+
out_y = empty(new_size, dtype=y.dtype)
314+
i_ = 0
315+
j = 0
316+
for i in range(nb):
317+
if j < nb_indexs and i == indexs[j]:
318+
j += 1
319+
cor = 360 if x[i - 1] > x[i] else -360
320+
out_x[i + i_] = (x[i] - ref) % 360 + ref - cor
321+
out_y[i + i_] = y[i]
322+
out_x[i + i_ + 1] = nan
323+
out_y[i + i_ + 1] = nan
324+
out_x[i + i_ + 2] = (x[i - 1] - ref) % 360 + ref + cor
325+
out_y[i + i_ + 2] = y[i - 1]
326+
i_ += 3
327+
out_x[i + i_] = (x[i] - ref) % 360 + ref
328+
out_y[i + i_] = y[i]
329+
return out_x, out_y
330+
331+
else:
332+
nb = x.shape[0]
333+
out = empty(nb, dtype=x.dtype)
334+
for i in range(nb):
335+
out[i] = (x[i] - ref) % 360 + ref
336+
return out, y
337+

src/py_eddy_tracker/observations/observation.py

Lines changed: 11 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@
5858
from tokenize import TokenError
5959
from matplotlib.path import Path as BasePath
6060
from .. import VAR_DESCR, VAR_DESCR_inv
61-
from ..generic import distance_grid, distance, flatten_line_matrix
61+
from ..generic import distance_grid, distance, flatten_line_matrix, wrap_longitude
6262
from ..poly import bbox_intersection, common_area
6363

6464
logger = logging.getLogger("pet")
@@ -1242,7 +1242,7 @@ def set_global_attr_netcdf(self, h_nc):
12421242
for key, item in self.global_attr.items():
12431243
h_nc.setncattr(key, item)
12441244

1245-
def scatter(self, ax, name, ref=None, factor=1, ** kwargs):
1245+
def scatter(self, ax, name, ref=None, factor=1, **kwargs):
12461246
x = self.longitude
12471247
if ref is not None:
12481248
x = (x - ref) % 360 + ref
@@ -1260,17 +1260,15 @@ def display(self, ax, ref=None, extern_only=False, intern_only=False, **kwargs):
12601260
kwargs_e = kwargs.copy()
12611261
if not extern_only:
12621262
kwargs_e.pop("label", None)
1263-
if ref is None:
1264-
if not extern_only:
1265-
ax.plot(lon_s, lat_s, **kwargs)
1266-
if not intern_only:
1267-
ax.plot(lon_e, lat_e, linestyle="-.", **kwargs_e)
1268-
else:
1269-
# FIXME : ref could split eddies
1270-
if not extern_only:
1271-
ax.plot((lon_s - ref) % 360 + ref, lat_s, **kwargs)
1272-
if not intern_only:
1273-
ax.plot((lon_e - ref) % 360 + ref, lat_e, linestyle="-.", **kwargs_e)
1263+
1264+
if not extern_only:
1265+
if ref is not None:
1266+
lon_s, lat_s = wrap_longitude(lon_s, lat_s, ref, cut=True)
1267+
ax.plot(lon_s, lat_s, **kwargs)
1268+
if not intern_only:
1269+
if ref is not None:
1270+
lon_e, lat_e = wrap_longitude(lon_e, lat_e, ref, cut=True)
1271+
ax.plot(lon_e, lat_e, linestyle="-.", **kwargs_e)
12741272

12751273
def grid_count(self, bins, intern=False, center=False):
12761274
x_name, y_name = self.intern(intern)

src/py_eddy_tracker/observations/tracking.py

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@
4444
from datetime import datetime, timedelta
4545
from .observation import EddiesObservations
4646
from numba import njit
47-
from ..generic import split_line
47+
from ..generic import split_line, wrap_longitude
4848

4949
logger = logging.getLogger("pet")
5050

@@ -365,10 +365,9 @@ def __extract_with_mask(
365365
def plot(self, ax, ref=None, **kwargs):
366366
if "label" in kwargs:
367367
kwargs["label"] += " (%s eddies)" % (self.nb_obs_by_track != 0).sum()
368-
x = self.longitude
368+
x, y = split_line(self.longitude, self.latitude, self.tracks)
369369
if ref is not None:
370-
x = (x - ref) % 360 + ref
371-
x, y = split_line(x, self.latitude, self.tracks)
370+
x, y = wrap_longitude(x, y, ref, cut=True)
372371
return ax.plot(x, y, **kwargs)
373372

374373

0 commit comments

Comments
 (0)