diff --git a/.readthedocs.yml b/.readthedocs.yml index ddfbb747..ba36f8ea 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -7,5 +7,7 @@ build: python: "mambaforge-latest" python: install: - - method: setuptools + - method: pip path: . +sphinx: + configuration: doc/conf.py \ No newline at end of file diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f8eee72f..6d6d6a30 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -11,14 +11,23 @@ and this project adheres to `Semantic Versioning =3.10, <3.13 - ffmpeg - pip - pip: + - -r ../requirements.txt + - git+https://github.com/AntSimi/py-eddy-tracker-sample-id.git - sphinx-gallery - sphinx_rtd_theme - sphinx>=3.1 diff --git a/environment.yml b/environment.yml index 12ce70e7..e94c7bc1 100644 --- a/environment.yml +++ b/environment.yml @@ -1,9 +1,8 @@ name: binder-pyeddytracker channels: - conda-forge - - defaults dependencies: - - python=3.10 + - python>=3.10, <3.13 - pip - ffmpeg - pip: diff --git a/examples/02_eddy_identification/pet_eddy_detection_ACC.py b/examples/02_eddy_identification/pet_eddy_detection_ACC.py index 3d3d4ac1..d12c62f3 100644 --- a/examples/02_eddy_identification/pet_eddy_detection_ACC.py +++ b/examples/02_eddy_identification/pet_eddy_detection_ACC.py @@ -7,6 +7,7 @@ Two detections are provided : with a filtered ADT and without filtering """ + from datetime import datetime from matplotlib import pyplot as plt, style @@ -80,7 +81,7 @@ def set_fancy_labels(fig, ticklabelsize=14, labelsize=14, labelweight="semibold" # Identification # ^^^^^^^^^^^^^^ # Run the identification step with slices of 2 mm -date = datetime(2016, 5, 15) +date = datetime(2019, 2, 23) kw_ident = dict( date=date, step=0.002, shape_error=70, sampling=30, uname="u", vname="v" ) diff --git a/requirements.txt b/requirements.txt index 4c8af099..556cabbf 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,11 +1,11 @@ -matplotlib +matplotlib < 3.8 # need an update of contour management opencv-python pint polygon3 pyyaml requests scipy -zarr +zarr < 3.0 netCDF4 numpy numba \ No newline at end of file diff --git a/src/py_eddy_tracker/appli/gui.py b/src/py_eddy_tracker/appli/gui.py index 4a8cdeb0..c3d7619b 100644 --- a/src/py_eddy_tracker/appli/gui.py +++ b/src/py_eddy_tracker/appli/gui.py @@ -11,7 +11,7 @@ from matplotlib.animation import FuncAnimation from matplotlib.axes import Axes from matplotlib.collections import LineCollection -from numpy import arange, where +from numpy import arange, where, nan from .. import EddyParser from ..gui import GUI @@ -58,7 +58,10 @@ def setup( self.kw_label["fontweight"] = kwargs.pop("fontweight", "demibold") # To text each visible eddy if field_txt: - self.field_txt = self.eddy[field_txt] + if isinstance(field_txt,str): + self.field_txt = self.eddy[field_txt] + else : + self.field_txt=field_txt if field_color: # To color each visible eddy self.field_color = self.eddy[field_color].astype("f4") diff --git a/src/py_eddy_tracker/appli/network.py b/src/py_eddy_tracker/appli/network.py index b8c2da51..0a3d06ca 100644 --- a/src/py_eddy_tracker/appli/network.py +++ b/src/py_eddy_tracker/appli/network.py @@ -283,7 +283,7 @@ def previous_obs(d, i_seg): def display_compare(ref, others): def display(value, ref=None): if ref: - outs = [f"{v/ref[k] * 100:.1f}% ({v})" for k, v in value.items()] + outs = [f"{v / ref[k] * 100:.1f}% ({v})" for k, v in value.items()] else: outs = value return "".join([f"{v:^18}" for v in outs]) diff --git a/src/py_eddy_tracker/data/__init__.py b/src/py_eddy_tracker/data/__init__.py index f14fee87..bf062983 100644 --- a/src/py_eddy_tracker/data/__init__.py +++ b/src/py_eddy_tracker/data/__init__.py @@ -8,6 +8,7 @@ 20160515 adt None None longitude latitude . \ --cut 800 --fil 1 """ + import io import lzma from os import path @@ -26,14 +27,20 @@ def get_remote_demo_sample(path): if path.endswith(".nc"): return io.BytesIO(content) else: - if path.endswith(".nc"): + try: + import py_eddy_tracker_sample_id + if path.endswith(".nc"): + return py_eddy_tracker_sample_id.get_remote_demo_sample(path) + content = open(py_eddy_tracker_sample_id.get_remote_demo_sample(f"{path}.tar.xz"), "rb").read() + except: + if path.endswith(".nc"): + content = requests.get( + f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}" + ).content + return io.BytesIO(content) content = requests.get( - f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}" + f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}.tar.xz" ).content - return io.BytesIO(content) - content = requests.get( - f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}.tar.xz" - ).content # Tar module could manage lzma tar, but it will apply uncompress for each extractfile tar = tarfile.open(mode="r", fileobj=io.BytesIO(lzma.decompress(content))) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index edb96bac..f15503b2 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -9,6 +9,7 @@ from matplotlib.path import Path as BasePath from netCDF4 import Dataset from numba import njit, prange, types as numba_types +import numpy as np from numpy import ( arange, array, @@ -35,7 +36,6 @@ percentile, pi, radians, - round_, sin, sinc, sqrt, @@ -1670,7 +1670,7 @@ def spectrum_lonlat(self, grid_name, area=None, ref=None, **kwargs): (lat_content[0], lat_content[1] / ref_lat_content[1]), ) - def compute_finite_difference(self, data, schema=1, mode="reflect", vertical=False): + def compute_finite_difference(self, data, schema=1, mode="reflect", vertical=False, second=False): if not isinstance(schema, int) and schema < 1: raise Exception("schema must be a positive int") @@ -1694,13 +1694,16 @@ def compute_finite_difference(self, data, schema=1, mode="reflect", vertical=Fal data2[:schema] = nan # Distance for one degree - d = self.EARTH_RADIUS * 2 * pi / 360 + d = self.EARTH_RADIUS * 2 * pi / 360 * 2 * schema # Mulitply by 2 step if vertical: - d *= self.ystep * 2 * schema + d *= self.ystep else: - d *= self.xstep * cos(deg2rad(self.y_c)) * 2 * schema - return (data1 - data2) / d + d *= self.xstep * cos(deg2rad(self.y_c)) + if second: + return (data1 + data2 - 2 * data) / (d ** 2 / 4) + else: + return (data1 - data2) / d def compute_stencil( self, data, stencil_halfwidth=4, mode="reflect", vertical=False @@ -1787,13 +1790,14 @@ def compute_stencil( ) return ma.array(g, mask=m) - def add_uv_lagerloef(self, grid_height, uname="u", vname="v", schema=15): - self.add_uv(grid_height, uname, vname) + def add_uv_lagerloef(self, grid_height, uname="u", vname="v", schema=15, **kwargs): + self.add_uv(grid_height, uname, vname, **kwargs) latmax = 5 - _, (i_start, i_end) = self.nearest_grd_indice((0, 0), (-latmax, latmax)) + _, i_start = self.nearest_grd_indice(0, -latmax) + _, i_end = self.nearest_grd_indice(0, latmax) sl = slice(i_start, i_end) # Divide by sideral day - lat = self.y_c[sl] + lat = self.y_c gob = ( cos(deg2rad(lat)) * ones((self.x_c.shape[0], 1)) @@ -1807,39 +1811,26 @@ def add_uv_lagerloef(self, grid_height, uname="u", vname="v", schema=15): mode = "wrap" if self.is_circular() else "reflect" # fill data to compute a finite difference on all point - data = self.convolve_filter_with_dynamic_kernel( - grid_height, - self.kernel_bessel, - lat_max=10, - wave_length=500, - order=1, - extend=0.1, - ) - data = self.convolve_filter_with_dynamic_kernel( - data, self.kernel_bessel, lat_max=10, wave_length=500, order=1, extend=0.1 - ) - data = self.convolve_filter_with_dynamic_kernel( - data, self.kernel_bessel, lat_max=10, wave_length=500, order=1, extend=0.1 - ) + kw_filter = dict(kernel_func=self.kernel_bessel, order=1, extend=.1) + data = self.convolve_filter_with_dynamic_kernel(grid_height, wave_length=500, **kw_filter, lat_max=6+5+2+3) v_lagerloef = ( self.compute_finite_difference( - self.compute_finite_difference(data, mode=mode, schema=schema), - mode=mode, - schema=schema, - )[:, sl] - * gob - ) - u_lagerloef = ( - -self.compute_finite_difference( - self.compute_finite_difference(data, vertical=True, schema=schema), - vertical=True, - schema=schema, - )[:, sl] + self.compute_finite_difference(data, mode=mode, schema=1), + vertical=True, schema=1 + ) * gob ) - w = 1 - exp(-((lat / 2.2) ** 2)) - self.vars[vname][:, sl] = self.vars[vname][:, sl] * w + v_lagerloef * (1 - w) - self.vars[uname][:, sl] = self.vars[uname][:, sl] * w + u_lagerloef * (1 - w) + u_lagerloef = -self.compute_finite_difference(data, vertical=True, schema=schema, second=True) * gob + + v_lagerloef = self.convolve_filter_with_dynamic_kernel(v_lagerloef, wave_length=195, **kw_filter, lat_max=6 + 5 +2) + v_lagerloef = self.convolve_filter_with_dynamic_kernel(v_lagerloef, wave_length=416, **kw_filter, lat_max=6 + 5) + v_lagerloef = self.convolve_filter_with_dynamic_kernel(v_lagerloef, wave_length=416, **kw_filter, lat_max=6) + u_lagerloef = self.convolve_filter_with_dynamic_kernel(u_lagerloef, wave_length=195, **kw_filter, lat_max=6 + 5 +2) + u_lagerloef = self.convolve_filter_with_dynamic_kernel(u_lagerloef, wave_length=416, **kw_filter, lat_max=6 + 5) + u_lagerloef = self.convolve_filter_with_dynamic_kernel(u_lagerloef, wave_length=416, **kw_filter, lat_max=6) + w = 1 - exp(-((lat[sl] / 2.2) ** 2)) + self.vars[vname][:, sl] = self.vars[vname][:, sl] * w + v_lagerloef[:, sl] * (1 - w) + self.vars[uname][:, sl] = self.vars[uname][:, sl] * w + u_lagerloef[:, sl] * (1 - w) def add_uv(self, grid_height, uname="u", vname="v", stencil_halfwidth=4): r"""Compute a u and v grid @@ -2260,12 +2251,11 @@ def compute_pixel_path(x0, y0, x1, y1, x_ori, y_ori, x_step, y_step, nb_x): i_x1 = empty(nx, dtype=numba_types.int_) i_y0 = empty(nx, dtype=numba_types.int_) i_y1 = empty(nx, dtype=numba_types.int_) - # Because round_ is not accepted with array in numba for i in range(nx): - i_x0[i] = round_(((x0[i] - x_ori) % 360) / x_step) - i_x1[i] = round_(((x1[i] - x_ori) % 360) / x_step) - i_y0[i] = round_((y0[i] - y_ori) / y_step) - i_y1[i] = round_((y1[i] - y_ori) / y_step) + i_x0[i] = np.round(((x0[i] - x_ori) % 360) / x_step) + i_x1[i] = np.round(((x1[i] - x_ori) % 360) / x_step) + i_y0[i] = np.round((y0[i] - y_ori) / y_step) + i_y1[i] = np.round((y1[i] - y_ori) / y_step) # Delta index of x d_x = i_x1 - i_x0 d_x = (d_x + nb_x // 2) % nb_x - (nb_x // 2) @@ -2950,7 +2940,7 @@ def compute_stencil(x, y, h, m, earth_radius, vertical=False, stencil_halfwidth= h_3, h_2, h_1, h0 = h[-4, j], h[-3, j], h[-2, j], h[-1, j] m_3, m_2, m_1, m0 = m[-4, j], m[-3, j], m[-2, j], m[-1, j] else: - m_3, m_2, m_1, m0 = False, False, False, False + m_3, m_2, m_1, m0 = True, True, True, True h1, h2, h3, h4 = h[0, j], h[1, j], h[2, j], h[3, j] m1, m2, m3, m4 = m[0, j], m[1, j], m[2, j], m[3, j] for i in range(nb_x): diff --git a/src/py_eddy_tracker/generic.py b/src/py_eddy_tracker/generic.py index 612def68..2fdb737a 100644 --- a/src/py_eddy_tracker/generic.py +++ b/src/py_eddy_tracker/generic.py @@ -615,7 +615,6 @@ def build_circle(x0, y0, r): return x_norm * r + x0, y_norm * r + y0 -@njit(cache=True) def window_index(x, x0, half_window=1): """ Give for a fixed half_window each start and end index for each x0, in @@ -626,7 +625,12 @@ def window_index(x, x0, half_window=1): :param float half_window: half window """ # Sort array, bounds will be sort also - i_ordered = x.argsort() + i_ordered = x.argsort(kind="mergesort") + return window_index_(x, i_ordered, x0, half_window) + + +@njit(cache=True) +def window_index_(x, i_ordered, x0, half_window=1): nb_x, nb_pt = x.size, x0.size first_index = empty(nb_pt, dtype=i_ordered.dtype) last_index = empty(nb_pt, dtype=i_ordered.dtype) diff --git a/src/py_eddy_tracker/observations/network.py b/src/py_eddy_tracker/observations/network.py index a2e2daed..f0b9d7cc 100644 --- a/src/py_eddy_tracker/observations/network.py +++ b/src/py_eddy_tracker/observations/network.py @@ -5,16 +5,19 @@ from glob import glob import logging import time - +from datetime import timedelta, datetime +import os import netCDF4 from numba import njit, types as nb_types from numba.typed import List +import numpy as np from numpy import ( arange, array, bincount, bool_, concatenate, + empty, nan, ones, @@ -119,13 +122,15 @@ def __repr__(self): period = (self.period[1] - self.period[0]) / 365.25 nb_by_network = self.network_size() nb_trash = 0 if self.ref_index != 0 else nb_by_network[0] + lifetime=self.lifetime big = 50_000 infos = [ f"Atlas with {self.nb_network} networks ({self.nb_network / period:0.0f} networks/year)," f" {self.nb_segment} segments ({self.nb_segment / period:0.0f} segments/year), {len(self)} observations ({len(self) / period:0.0f} observations/year)", f" {m_event.size} merging ({m_event.size / period:0.0f} merging/year), {s_event.size} splitting ({s_event.size / period:0.0f} splitting/year)", - f" with {(nb_by_network > big).sum()} network with more than {big} obs and the biggest have {nb_by_network.max()} observations ({nb_by_network[nb_by_network> big].sum()} observations cumulate)", + f" with {(nb_by_network > big).sum()} network with more than {big} obs and the biggest have {nb_by_network.max()} observations ({nb_by_network[nb_by_network > big].sum()} observations cumulate)", f" {nb_trash} observations in trash", + f" {lifetime.max()} days max of lifetime", ] return "\n".join(infos) @@ -200,6 +205,13 @@ def ref_segment_track_index(self): @property def ref_index(self): return self.index_network[2] + + @property + def lifetime(self): + """Return lifetime for each observation""" + lt=self.networks_period.astype("int") + nb_by_network=self.network_size() + return lt.repeat(nb_by_network) def network_segment_size(self, id_networks=None): """Get number of segment by network @@ -225,6 +237,15 @@ def network_size(self, id_networks=None): i = id_networks - self.index_network[2] return self.index_network[1][i] - self.index_network[0][i] + @property + def networks_period(self): + """ + Return period for each network + """ + return get_period_with_index(self.time, *self.index_network[:2]) + + + def unique_segment_to_id(self, id_unique): """Return id network and id segment for a unique id @@ -274,7 +295,7 @@ def astype(self, cls): new[k][:] = self[k][:] new.sign_type = self.sign_type return new - + def longer_than(self, nb_day_min=-1, nb_day_max=-1): """ Select network on time duration @@ -1125,23 +1146,29 @@ def segment_track_array(self): self._segment_track_array = build_unique_array(self.segment, self.track) return self._segment_track_array - def birth_event(self): + def birth_event(self, only_index=False): """Extract birth events.""" i_start, _, _ = self.index_segment_track indices = i_start[self.previous_obs[i_start] == -1] if self.first_is_trash(): indices = indices[1:] - return self.extract_event(indices) - + if only_index: + return indices + else : + return self.extract_event(indices) + generation_event = birth_event - def death_event(self): + def death_event(self, only_index=False): """Extract death events.""" _, i_stop, _ = self.index_segment_track indices = i_stop[self.next_obs[i_stop - 1] == -1] - 1 if self.first_is_trash(): indices = indices[1:] - return self.extract_event(indices) + if only_index: + return indices + else : + return self.extract_event(indices) dissipation_event = death_event @@ -1452,7 +1479,7 @@ def plot(self, ax, ref=None, color_cycle=None, **kwargs): j += 1 return mappables - def remove_dead_end(self, nobs=3, ndays=0, recursive=0, mask=None): + def remove_dead_end(self, nobs=3, ndays=0, recursive=0, mask=None, return_mask=False): """ Remove short segments that don't connect several segments @@ -1478,6 +1505,8 @@ def remove_dead_end(self, nobs=3, ndays=0, recursive=0, mask=None): ) # get mask for selected obs m = ~self.segment_mask(segments_keep) + if return_mask: + return ~m self.track[m] = 0 self.segment[m] = 0 self.previous_obs[m] = -1 @@ -1495,6 +1524,8 @@ def remove_dead_end(self, nobs=3, ndays=0, recursive=0, mask=None): self.sort() if recursive > 0: self.remove_dead_end(nobs, ndays, recursive - 1) + + def extract_segment(self, segments, absolute=False): """Extract given segments @@ -1788,8 +1819,8 @@ def date2file(julian_day): ) logger.info( ( - f"coherence {_t} / {range_end-1} ({(_t - range_start) / (range_end - range_start-1):.1%})" - f" : {time.time()-_timestamp:5.2f}s" + f"coherence {_t} / {range_end - 1} ({(_t - range_start) / (range_end - range_start - 1):.1%})" + f" : {time.time() - _timestamp:5.2f}s" ) ) @@ -1865,8 +1896,8 @@ def date2file(julian_day): ) logger.info( ( - f"coherence {_t} / {range_end-1} ({(_t - range_start) / (range_end - range_start-1):.1%})" - f" : {time.time()-_timestamp:5.2f}s" + f"coherence {_t} / {range_end - 1} ({(_t - range_start) / (range_end - range_start - 1):.1%})" + f" : {time.time() - _timestamp:5.2f}s" ) ) return itf_final, ptf_final @@ -2035,6 +2066,7 @@ def group_observations(self, min_overlap=0.2, minimal_area=False, **kwargs): results, nb_obs = list(), list() # To display print only in INFO display_iteration = logger.getEffectiveLevel() == logging.INFO + for i, filename in enumerate(self.filenames): if display_iteration: print(f"{filename} compared to {self.window} next", end="\r") @@ -2065,7 +2097,7 @@ def group_observations(self, min_overlap=0.2, minimal_area=False, **kwargs): nb_alone, nb_obs, nb_gr = (gr == self.NOGROUP).sum(), len(gr), len(unique(gr)) logger.info( f"{nb_alone} alone / {nb_obs} obs, {nb_gr} groups, " - f"{nb_alone *100./nb_obs:.2f} % alone, {(nb_obs - nb_alone) / (nb_gr - 1):.1f} obs/group" + f"{nb_alone * 100. / nb_obs:.2f} % alone, {(nb_obs - nb_alone) / (nb_gr - 1):.1f} obs/group" ) return gr @@ -2316,3 +2348,21 @@ def mask_follow_obs(m, next_obs, time, indexs, dt=3): m[i_next] = True i_next = next_obs[i_next] dt_ = abs(time[i_next] - t0) + + +@njit(cache=True) +def get_period_with_index(t, i0, i1): + """Return peek to peek cover by each slice define by i0 and i1 + + :param array t: array which contain values to estimate spread + :param array i0: index which determine start of slice + :param array i1: index which determine end of slice + :return array: Peek to peek of t + """ + periods = np.empty(i0.size, t.dtype) + for i in range(i0.size): + if i1[i] == i0[i]: + periods[i] = 0 + continue + periods[i] = t[i0[i] : i1[i]].ptp() + return periods diff --git a/src/py_eddy_tracker/observations/tracking.py b/src/py_eddy_tracker/observations/tracking.py index 164f9724..fa1c1f93 100644 --- a/src/py_eddy_tracker/observations/tracking.py +++ b/src/py_eddy_tracker/observations/tracking.py @@ -380,7 +380,7 @@ def extract_toward_direction(self, west=True, delta_lon=None): d_lon = lon[i1] - lon[i0] m = d_lon < 0 if west else d_lon > 0 if delta_lon is not None: - m *= delta_lon < d_lon + m *= delta_lon < abs(d_lon) m = m.repeat(nb) return self.extract_with_mask(m)