From 58fdba736074b49c458f63901542445644eaf401 Mon Sep 17 00:00:00 2001 From: Cori Pegliasco Date: Thu, 4 Feb 2021 18:35:19 +0100 Subject: [PATCH 1/4] - add scatter and event_map --- src/py_eddy_tracker/observations/network.py | 80 +++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/src/py_eddy_tracker/observations/network.py b/src/py_eddy_tracker/observations/network.py index f1e0a9a9..d9ef806b 100644 --- a/src/py_eddy_tracker/observations/network.py +++ b/src/py_eddy_tracker/observations/network.py @@ -407,6 +407,86 @@ def scatter_timeline( mappables["scatter"] = ax.scatter(self.time, y, **kwargs) return mappables + def event_map(self, ax, **kwargs): + """Add the merging and splitting events """ + j = 0 + mappables = dict() + symbol_kw = dict( + markersize=10, + color="k", + ) + symbol_kw.update(kwargs) + symbol_kw_split = symbol_kw.copy() + symbol_kw_split["markersize"] += 5 + for i, b0, b1 in self.iter_on("segment"): + nb = i.stop - i.start + if nb == 0: + continue + event_kw = dict(color=self.COLORS[j % self.NB_COLORS], ls="-", **kwargs) + i_n, i_p = ( + self.next_obs[i.stop - 1], + self.previous_obs[i.start], + ) + + if i_n != -1: + y0, y1 = self.lat[i.stop - 1], self.lat[i_n] + x0, x1 = self.lon[i.stop - 1], self.lon[i_n] + ax.plot((x0, x1), (y0, y1), **event_kw)[0] + ax.plot(x0, y0, marker="s", **symbol_kw)[0] + if i_p != -1: + y0, y1 = self.lat[i.start], self.lat[i_p] + x0, x1 = self.lon[i.start], self.lon[i_p] + ax.plot((x0, x1), (y0, y1), **event_kw)[0] + ax.plot(x0, y0, marker="*", **symbol_kw_split)[0] + + j += 1 + return mappables + + def scatter( + self, + ax, + name="time", + s=25, + factor=1, + ref=None, + edgecolor_cycle=None, + cmap="Spectral_r", + **kwargs, + ): + """ + This function will scatter the path of each trajectory + + :param matplotlib.axes.Axes ax: ax to draw + :param float,int ref: if defined, all coordinates will be wrapped with ref like west boundary + :param dict kwargs: keyword arguments for Axes.plot + :return: a list of matplotlib mappables + """ + mappables = dict() + nb_colors = len(edgecolor_cycle) if edgecolor_cycle != None else None + x = self.longitude + if ref is not None: + x = (x - ref) % 360 + ref + kwargs = kwargs.copy() + if nb_colors: + edgecolors = list() + seg_previous = self.segment[0] + j = 0 + for seg in self.segment: + if seg != seg_previous: + j += 1 + edgecolors.append(edgecolor_cycle[j % nb_colors]) + seg_previous = seg + mappables["edges"] = ax.scatter( + x, self.latitude, edgecolor=edgecolors, **kwargs + ) + kwargs.pop("linewidths", None) + kwargs["lw"] = 0 + if name is not None and "c" not in kwargs: + v = self.parse_varname(name) + kwargs["c"] = v * factor + mappables["scatter"] = ax.scatter(x, self.latitude, **kwargs) + return mappables + def insert_virtual(self): # TODO pass From 7be041e3cb8e15e48ef0e174e2a7ae6513f0dba1 Mon Sep 17 00:00:00 2001 From: Cori Pegliasco Date: Thu, 4 Feb 2021 18:50:40 +0100 Subject: [PATCH 2/4] add documentation --- src/py_eddy_tracker/observations/network.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/py_eddy_tracker/observations/network.py b/src/py_eddy_tracker/observations/network.py index d9ef806b..2babaf9e 100644 --- a/src/py_eddy_tracker/observations/network.py +++ b/src/py_eddy_tracker/observations/network.py @@ -446,23 +446,25 @@ def scatter( self, ax, name="time", - s=25, factor=1, ref=None, edgecolor_cycle=None, - cmap="Spectral_r", **kwargs, ): """ - This function will scatter the path of each trajectory - - :param matplotlib.axes.Axes ax: ax to draw - :param float,int ref: if defined, all coordinates will be wrapped with ref like west boundary - :param dict kwargs: keyword arguments for Axes.plot - :return: a list of matplotlib mappables + This function will scatter the path of each network, with the merging and splitting events + + :param matplotlib.axes.Axes ax: matplotlib axe used to draw + :param str,array,None name: + variable used to fill the contour, if None all elements have the same color + :param float,None ref: if define use like west bound + :param float factor: multiply value by + :param list edgecolor_cycle: list of colors + :param dict kwargs: look at :py:meth:`matplotlib.axes.Axes.scatter` + :return: a dict of scattered mappables """ mappables = dict() - nb_colors = len(edgecolor_cycle) if edgecolor_cycle != None else None + nb_colors = len(edgecolor_cycle) if edgecolor_cycle else None x = self.longitude if ref is not None: x = (x - ref) % 360 + ref From b7964ebf5ef3b07a291afd6fa6eaf5bb53bb9cb5 Mon Sep 17 00:00:00 2001 From: Cori Pegliasco Date: Fri, 5 Feb 2021 18:55:51 +0100 Subject: [PATCH 3/4] - Add documentation in examples - change marker --- examples/16_network/pet_atlas.py | 32 +++++++++++++--- examples/16_network/pet_ioannou_2017_case.py | 28 +++++++++++--- examples/16_network/pet_relative.py | 38 ++++++++++++++++--- .../16_network/pet_replay_segmentation.py | 20 ++++++---- src/py_eddy_tracker/appli/gui.py | 6 ++- src/py_eddy_tracker/observations/network.py | 6 +-- 6 files changed, 103 insertions(+), 27 deletions(-) diff --git a/examples/16_network/pet_atlas.py b/examples/16_network/pet_atlas.py index 687e5d24..83dac5b1 100644 --- a/examples/16_network/pet_atlas.py +++ b/examples/16_network/pet_atlas.py @@ -1,6 +1,6 @@ """ -Network analyse -=============== +Network Analysis +================ """ from matplotlib import pyplot as plt from numpy import ma @@ -21,7 +21,7 @@ # %% -# Function +# Functions def start_axes(title): fig = plt.figure(figsize=(13, 5)) ax = fig.add_axes([0.03, 0.03, 0.90, 0.94], projection="full_axes") @@ -41,6 +41,7 @@ def update_axes(ax, mappable=None): # %% # All # --- +# Display the % of time each pixel (1/10°) is within an anticyclonic network ax = start_axes("") g_all = n.grid_count(bins) m = g_all.display(ax, **kw_time, vmin=0, vmax=75) @@ -49,12 +50,16 @@ def update_axes(ax, mappable=None): # %% # Network longer than 10 days # --------------------------- +# Display the % of time each pixel (1/10°) is within an anticyclonic network +# which total lifetime in longer than 10 days ax = start_axes("") n10 = n.longer_than(10) g_10 = n10.grid_count(bins) m = g_10.display(ax, **kw_time, vmin=0, vmax=75) update_axes(ax, m).set_label("Pixel used in % of time") +# Display the ratio between the short and total presence. +# Red = mostly short networks ax = start_axes("") m = g_10.display( ax, @@ -67,12 +72,17 @@ def update_axes(ax, mappable=None): # %% # Network longer than 20 days # --------------------------- +# Display the % of time each pixel (1/10°) is within an anticyclonic network +# which total lifetime in longer than 20 days ax = start_axes("") n20 = n.longer_than(20) g_20 = n20.grid_count(bins) m = g_20.display(ax, **kw_time, vmin=0, vmax=75) update_axes(ax, m).set_label("Pixel used in % of time") +# %% +# Display the ratio between the short and total presence. +# Red = mostly short networks ax = start_axes("") m = g_20.display( ax, @@ -82,6 +92,11 @@ def update_axes(ax, mappable=None): name=g_20.vars["count"] * 100.0 / g_all.vars["count"] ) update_axes(ax, m).set_label("Pixel used in % all atlas") + +# %% +# Display the ratio between the short and total presence. +# Red = mostly short networks +# Networks shorter than 365 days are masked ax = start_axes("") m = g_20.display( ax, @@ -93,6 +108,11 @@ def update_axes(ax, mappable=None): ) ) update_axes(ax, m).set_label("Pixel used in % all atlas") +# %% +# Display the ratio between the short and total presence. +# Red = mostly short networks +# Networks longer than 365 days are masked +# -> Coastal areas are mostly populated by short networks ax = start_axes("") m = g_20.display( ax, @@ -110,11 +130,13 @@ def update_axes(ax, mappable=None): # %% # All merging # ----------- +# Display the occurence of merging events ax = start_axes("") g_all_merging = n.merging_event().grid_count(bins) m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=0.2) update_axes(ax, m).set_label("Pixel used in % of time") # %% +# Ratio merging events / eddy presence ax = start_axes("") m = g_all_merging.display( ax, @@ -126,8 +148,8 @@ def update_axes(ax, mappable=None): update_axes(ax, m).set_label("Pixel used in % all atlas") # %% -# Merging in network longer than 10 days -# -------------------------------------- +# Merging in networks longer than 10 days +# --------------------------------------- ax = start_axes("") g_10_merging = n10.merging_event().grid_count(bins) m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=0.2) diff --git a/examples/16_network/pet_ioannou_2017_case.py b/examples/16_network/pet_ioannou_2017_case.py index 906b90eb..93e0b4ef 100644 --- a/examples/16_network/pet_ioannou_2017_case.py +++ b/examples/16_network/pet_ioannou_2017_case.py @@ -4,12 +4,18 @@ Figure 10 from https://doi.org/10.1002/2017JC013158 """ + +# %% +# We want to find the Ierapetra Eddy described above in the networks + +# %% from datetime import datetime, timedelta import numpy as np from matplotlib import pyplot as plt from matplotlib.animation import FuncAnimation from matplotlib.ticker import FuncFormatter +from matplotlib import colors import py_eddy_tracker.gui from py_eddy_tracker.appli.gui import Anim @@ -62,6 +68,8 @@ def update_axes(ax, mappable=None): # %% +# We know the position and the time of a specific eddy +# `n.extract_with_mask` give us the corresponding network n = NetworkObservations.load_file( "med/Anticyclonic_seg.nc" ) @@ -77,6 +85,7 @@ def update_axes(ax, mappable=None): print(ioannou_case.infos()) # %% +# It seems that this network is huge! Our case is in purple... ax = start_axes() ioannou_case.plot(ax) update_axes(ax) @@ -84,6 +93,7 @@ def update_axes(ax, mappable=None): # %% # Full Timeline # ------------- +# The network span for many years... How to cut the interesting part? fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.05, 0.92, 0.92]) ax.xaxis.set_major_formatter(formatter), ax.grid() @@ -93,6 +103,7 @@ def update_axes(ax, mappable=None): # %% # Sub network and new numbering # ----------------------------- +# Here we chose to keep only the order 3 segments relatives to our chosen eddy i = np.where( (ioannou_case.lat > 33) * (ioannou_case.lat < 34) @@ -107,10 +118,14 @@ def update_axes(ax, mappable=None): # %% # Anim # ---- +# Quick movie to see better! +cmap = colors.ListedColormap( + list(close_to_i3.COLORS), name="from_list", N=close_to_i3.segment.max() +) a = Anim( close_to_i3, figsize=(12, 4), - cmap="Spectral_r", + cmap=cmap, nb_step=7, dpi=55, field_color="segment", @@ -136,8 +151,8 @@ def update_axes(ax, mappable=None): update_axes(ax) # %% -# Local Timeline -# -------------- +# Latitude Timeline +# ----------------- ax = timeline_axes(f"Close segments ({close_to_i3.infos()})") n_copy = close_to_i3.copy() n_copy.median_filter(15, "time", "latitude") @@ -146,6 +161,7 @@ def update_axes(ax, mappable=None): # %% # Local radius timeline # --------------------- +# Effective (bold) and Speed (thin) Radius together n_copy.median_filter(2, "time", "radius_e") n_copy.median_filter(2, "time", "radius_s") for b0, b1 in [ @@ -167,14 +183,16 @@ def update_axes(ax, mappable=None): ) # %% -# Parameter timeline -# ------------------ +# Parameters timeline +# ------------------- +# Effective Radius kw = dict(s=35, cmap=plt.get_cmap("Spectral_r", 8), zorder=10) ax = timeline_axes() m = close_to_i3.scatter_timeline(ax, "radius_e", factor=1e-3, vmin=20, vmax=100, **kw) cb = update_axes(ax, m["scatter"]) cb.set_label("Effective radius (km)") # %% +# Shape error ax = timeline_axes() m = close_to_i3.scatter_timeline(ax, "shape_error_e", vmin=14, vmax=70, **kw) cb = update_axes(ax, m["scatter"]) diff --git a/examples/16_network/pet_relative.py b/examples/16_network/pet_relative.py index 5fac0fc2..edb537d5 100644 --- a/examples/16_network/pet_relative.py +++ b/examples/16_network/pet_relative.py @@ -18,7 +18,10 @@ e.lon[:] = (e.lon + 180) % 360 - 180 e.contour_lon_e[:] = ((e.contour_lon_e.T - e.lon + 180) % 360 - 180 + e.lon).T e.contour_lon_s[:] = ((e.contour_lon_s.T - e.lon + 180) % 360 - 180 + e.lon).T -############## +# %% +# Do segmentation +# --------------- +# Segmentation based on maximum overlap, temporal window for candidates = 5 days n = NetworkObservations.from_split_network(e, e.split_network(intern=False, window=5)) # %% @@ -26,7 +29,9 @@ # -------- # %% -# Display timeline +# Display timeline with events +# A segment generated by a splitting is marked with a star +# A segment merging in another is marked with an exagon fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.04, 0.92, 0.92]) _ = n.display_timeline(ax) @@ -38,27 +43,34 @@ _ = n.display_timeline(ax, event=False) # %% -# Timeline by latitude mean +# Timeline by mean latitude # ------------------------- +# Display timeline with the mean latitude of the segments in yaxis fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.04, 0.92, 0.92]) +ax.set_ylabel("Latitude") _ = n.display_timeline(ax, field="latitude") # %% -# Timeline by radius mean +# Timeline by mean Effective Radius # ----------------------- +# The factor argument is applied on the chosen field fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.04, 0.92, 0.92]) -_ = n.display_timeline(ax, field="radius_e") +ax.set_ylabel("Effective Radius (km)") +_ = n.display_timeline(ax, field="radius_e", factor=1e-3) # %% # Timeline by latitude # -------------------- +# Use `method="all"` to display the consecutive values of the field fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.05, 0.92, 0.92]) +ax.set_ylabel("Latitude") _ = n.display_timeline(ax, field="lat", method="all") # %% +# You can filter the data, here with a time window of 15 days fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.05, 0.92, 0.92]) n_copy = n.copy() @@ -68,6 +80,8 @@ # %% # Parameters timeline # ------------------- +# Scatter is usefull to display the parameters' temporal evolution +# Effective Radius and Amplitude kw = dict(s=25, cmap="Spectral_r", zorder=10) fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.54, 0.90, 0.44]) @@ -85,6 +99,7 @@ cb.set_label("Amplitude (cm)") # %% +# Speed fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) m = n.scatter_timeline(ax, "speed_average", factor=100, vmin=0, vmax=40, **kw) @@ -94,6 +109,7 @@ cb.set_label("Maximum speed (cm/s)") # %% +# Speed Radius fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) m = n.scatter_timeline(ax, "radius_s", factor=1e-3, vmin=20, vmax=100, **kw) @@ -105,7 +121,7 @@ # %% # Remove dead branch # ------------------ -# Remove all tiny segment with less than N obs which didn't join two segments +# Remove all tiny segments with less than N obs which didn't join two segments # # .. warning:: # Must be explore, no solution to solve all the case @@ -122,6 +138,7 @@ # %% # Keep close relative # ------------------- +# When you want to investigate one particular observation and select only the closest segments # First choose an observation in the network fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) @@ -132,6 +149,7 @@ _ = ax.plot(*obs_args, **obs_kw) # %% +# Colors show the relative order of the segment with regards to the chosen one fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) m = n.scatter_timeline( @@ -143,18 +161,21 @@ ) cb.set_label("Relative order") # %% +# You want to keep only the segments at the order 1 fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) close_to_i1 = n.relative(i, order=1) ax.set_title(f"Close segments ({close_to_i1.infos()})") _ = close_to_i1.display_timeline(ax) # %% +# You want to keep the segments until order 2 fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) close_to_i2 = n.relative(i, order=2) ax.set_title(f"Close segments ({close_to_i2.infos()})") _ = close_to_i2.display_timeline(ax) # %% +# You want to keep the segments until order 3 fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) close_to_i3 = n.relative(i, order=3) @@ -164,6 +185,7 @@ # %% # Display track on map # -------------------- +# Only a map can be tricky to understand, with a timeline it's easier! fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.94, 0.88], projection="full_axes") close_to_i2.plot(ax, color_cycle=close_to_i2.COLORS) @@ -175,6 +197,7 @@ # %% # Get merging event # ----------------- +# Display the position of the eddies after a merging fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") merging = close_to_i2.merging_event() @@ -185,6 +208,7 @@ # %% # Get spliting event # ------------------ +# Display the position of the eddies before a splitting fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") spliting = close_to_i2.spliting_event() @@ -195,6 +219,7 @@ # %% # Get birth event # ------------------ +# Display the starting position of non-splitted eddies fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") birth = close_to_i2.birth_event() @@ -205,6 +230,7 @@ # %% # Get death event # ------------------ +# Display the last position of non-merged eddies fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") death = close_to_i2.death_event() diff --git a/examples/16_network/pet_replay_segmentation.py b/examples/16_network/pet_replay_segmentation.py index 67c89127..c90fca9d 100644 --- a/examples/16_network/pet_replay_segmentation.py +++ b/examples/16_network/pet_replay_segmentation.py @@ -4,12 +4,17 @@ Case from figure 10 from https://doi.org/10.1002/2017JC013158 """ + +# %% +# Again with the Ierapetra Eddy + from datetime import datetime, timedelta import numpy as np from matplotlib import pyplot as plt from matplotlib.animation import FuncAnimation from matplotlib.ticker import FuncFormatter +from matplotlib import colors import py_eddy_tracker.gui from py_eddy_tracker.appli.gui import Anim @@ -75,7 +80,7 @@ def follow_obs(cls, i_next, track_id, used, ids, *args, **kwargs): """ Method to overwrite behaviour in merging. - We will give the point to the older one + We will give the point to the older one instead of the maximum overlap ratio """ while i_next != -1: # Flag @@ -116,14 +121,14 @@ def get_obs(dataset): # %% # Get original network, we will isolate only relative at order *2* n = NetworkObservations.load_file( - "/tmp/Anticyclonic_seg.nc" - # "/data/adelepoulle/work/Eddies/20201217_network_build/tracking/med/Anticyclonic_seg.nc" + "/data/adelepoulle/work/Eddies/20201217_network_build/tracking/med/Anticyclonic_seg.nc" ) n = n.extract_with_mask(n.track == n.track[get_obs(n)]) n_ = n.relative(get_obs(n), order=2) # %% +# Display the default segmentation ax = start_axes(n_.infos()) n_.plot(ax, color_cycle=n.COLORS) update_axes(ax) @@ -143,8 +148,9 @@ def get_obs(dataset): n_.numbering_segment() # %% -# New version -# ----------- +# New segmentation +# ---------------- +# "The oldest wins" method produce a very long segment ax = start_axes(n_.infos()) n_.plot(ax, color_cycle=n_.COLORS) update_axes(ax) @@ -154,8 +160,8 @@ def get_obs(dataset): _ = n_.display_timeline(ax) # %% -# Parameter timeline -# ------------------ +# Parameters timeline +# ------------------- kw = dict(s=35, cmap=plt.get_cmap("Spectral_r", 8), zorder=10) ax = timeline_axes() n_.median_filter(15, "time", "latitude") diff --git a/src/py_eddy_tracker/appli/gui.py b/src/py_eddy_tracker/appli/gui.py index bbe5067c..e9ff7f7a 100644 --- a/src/py_eddy_tracker/appli/gui.py +++ b/src/py_eddy_tracker/appli/gui.py @@ -186,7 +186,11 @@ def update(self): # Update id txt for i in where(m)[0]: mappable = self.ax.text( - self.x_core[i], self.y_core[i], self.field_txt[i], fontsize=8 + self.x_core[i], + self.y_core[i], + self.field_txt[i], + fontsize=12, + fontweight="demibold", ) self.mappables.append(mappable) self.ax.draw_artist(mappable) diff --git a/src/py_eddy_tracker/observations/network.py b/src/py_eddy_tracker/observations/network.py index 2babaf9e..89c51b29 100644 --- a/src/py_eddy_tracker/observations/network.py +++ b/src/py_eddy_tracker/observations/network.py @@ -333,7 +333,7 @@ def event_timeline(self, ax, field=None, method=None, factor=1): ) ) ax.plot((x[-1], self.time[i_n]), (y0, y1), **event_kw)[0] - ax.plot(x[-1], y0, color="k", marker=">", markersize=10, zorder=-1)[0] + ax.plot(x[-1], y0, color="k", marker="H", markersize=10, zorder=-1)[0] if i_p != -1: seg_previous = self.segment[i_p] if field is not None and method == "all": @@ -417,7 +417,7 @@ def event_map(self, ax, **kwargs): ) symbol_kw.update(kwargs) symbol_kw_split = symbol_kw.copy() - symbol_kw_split["markersize"] += 5 + symbol_kw_split["markersize"] += 4 for i, b0, b1 in self.iter_on("segment"): nb = i.stop - i.start if nb == 0: @@ -432,7 +432,7 @@ def event_map(self, ax, **kwargs): y0, y1 = self.lat[i.stop - 1], self.lat[i_n] x0, x1 = self.lon[i.stop - 1], self.lon[i_n] ax.plot((x0, x1), (y0, y1), **event_kw)[0] - ax.plot(x0, y0, marker="s", **symbol_kw)[0] + ax.plot(x0, y0, marker="H", **symbol_kw)[0] if i_p != -1: y0, y1 = self.lat[i.start], self.lat[i_p] x0, x1 = self.lon[i.start], self.lon[i_p] From ed2848f89da9c9536cd8b0783f484ce1d75873b7 Mon Sep 17 00:00:00 2001 From: Cori Pegliasco Date: Wed, 10 Feb 2021 10:08:45 +0100 Subject: [PATCH 4/4] - add triplet in merging and splitting --- examples/16_network/pet_atlas.py | 27 ++++--- examples/16_network/pet_ioannou_2017_case.py | 1 + examples/16_network/pet_relative.py | 32 ++++++-- src/py_eddy_tracker/observations/network.py | 82 ++++++++++++++++---- 4 files changed, 111 insertions(+), 31 deletions(-) diff --git a/examples/16_network/pet_atlas.py b/examples/16_network/pet_atlas.py index 83dac5b1..cbb7b3e9 100644 --- a/examples/16_network/pet_atlas.py +++ b/examples/16_network/pet_atlas.py @@ -17,7 +17,7 @@ step = 1 / 10.0 bins = ((-10, 37, step), (30, 46, step)) kw_time = dict(cmap="terrain_r", factor=100.0 / n.nb_days, name="count") -kw_ratio = dict(cmap=plt.get_cmap("coolwarm_r", 10)) +kw_ratio = dict(cmap=plt.get_cmap("magma_r", 10)) # %% @@ -58,8 +58,10 @@ def update_axes(ax, mappable=None): m = g_10.display(ax, **kw_time, vmin=0, vmax=75) update_axes(ax, m).set_label("Pixel used in % of time") +# %% # Display the ratio between the short and total presence. -# Red = mostly short networks +# +# Light = mostly short networks ax = start_axes("") m = g_10.display( ax, @@ -73,7 +75,7 @@ def update_axes(ax, mappable=None): # Network longer than 20 days # --------------------------- # Display the % of time each pixel (1/10°) is within an anticyclonic network -# which total lifetime in longer than 20 days +# which total lifetime is longer than 20 days ax = start_axes("") n20 = n.longer_than(20) g_20 = n20.grid_count(bins) @@ -82,7 +84,8 @@ def update_axes(ax, mappable=None): # %% # Display the ratio between the short and total presence. -# Red = mostly short networks +# +# Light = mostly short networks ax = start_axes("") m = g_20.display( ax, @@ -95,7 +98,9 @@ def update_axes(ax, mappable=None): # %% # Display the ratio between the short and total presence. -# Red = mostly short networks +# +# Light = mostly short networks +# # Networks shorter than 365 days are masked ax = start_axes("") m = g_20.display( @@ -110,9 +115,10 @@ def update_axes(ax, mappable=None): update_axes(ax, m).set_label("Pixel used in % all atlas") # %% # Display the ratio between the short and total presence. -# Red = mostly short networks +# # Networks longer than 365 days are masked -# -> Coastal areas are mostly populated by short networks +# +# # -> Coastal areas are mostly populated by short networks ax = start_axes("") m = g_20.display( ax, @@ -133,8 +139,9 @@ def update_axes(ax, mappable=None): # Display the occurence of merging events ax = start_axes("") g_all_merging = n.merging_event().grid_count(bins) -m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=0.2) +m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=1) update_axes(ax, m).set_label("Pixel used in % of time") + # %% # Ratio merging events / eddy presence ax = start_axes("") @@ -152,7 +159,7 @@ def update_axes(ax, mappable=None): # --------------------------------------- ax = start_axes("") g_10_merging = n10.merging_event().grid_count(bins) -m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=0.2) +m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=1) update_axes(ax, m).set_label("Pixel used in % of time") # %% ax = start_axes("") @@ -160,7 +167,7 @@ def update_axes(ax, mappable=None): ax, **kw_ratio, vmin=0, - vmax=0.5, + vmax=2, name=ma.array( g_10_merging.vars["count"] * 100.0 / g_10.vars["count"], mask=g_10.vars["count"] < 365, diff --git a/examples/16_network/pet_ioannou_2017_case.py b/examples/16_network/pet_ioannou_2017_case.py index 93e0b4ef..2e7f807e 100644 --- a/examples/16_network/pet_ioannou_2017_case.py +++ b/examples/16_network/pet_ioannou_2017_case.py @@ -69,6 +69,7 @@ def update_axes(ax, mappable=None): # %% # We know the position and the time of a specific eddy +# # `n.extract_with_mask` give us the corresponding network n = NetworkObservations.load_file( "med/Anticyclonic_seg.nc" diff --git a/examples/16_network/pet_relative.py b/examples/16_network/pet_relative.py index 7b6a96f4..f743d3b4 100644 --- a/examples/16_network/pet_relative.py +++ b/examples/16_network/pet_relative.py @@ -32,6 +32,7 @@ # %% # Display timeline with events # A segment generated by a splitting is marked with a star +# # A segment merging in another is marked with an exagon fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.04, 0.92, 0.92]) @@ -82,6 +83,7 @@ # Parameters timeline # ------------------- # Scatter is usefull to display the parameters' temporal evolution +# # Effective Radius and Amplitude kw = dict(s=25, cmap="Spectral_r", zorder=10) fig = plt.figure(figsize=(15, 8)) @@ -125,7 +127,7 @@ # Remove all tiny segments with less than N obs which didn't join two segments # # .. warning:: -# Must be explore, no solution to solve all the case +# Must be explore, no solution to solve all cases n_clean = n.remove_dead_end(nobs=10) fig = plt.figure(figsize=(15, 8)) @@ -140,6 +142,7 @@ # Keep close relative # ------------------- # When you want to investigate one particular observation and select only the closest segments +# # First choose an observation in the network fig = plt.figure(figsize=(15, 5)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88]) @@ -201,10 +204,18 @@ # Display the position of the eddies after a merging fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") -merging = close_to_i2.merging_event() -merging.display(ax) +close_to_i2.plot(ax, color_cycle=close_to_i2.COLORS) +m1, m0, m0_stop = close_to_i2.merging_event(triplet=True) +m1.display(ax, color="violet", lw=2, label="Eddies after merging") +m0.display(ax, color="blueviolet", lw=2, label="Eddies before merging") +m0_stop.display(ax, color="black", lw=2, label="Eddies stopped by merging") +ax.plot(m1.lon, m1.lat, marker=".", color="purple", ls="") +ax.plot(m0.lon, m0.lat, marker=".", color="blueviolet", ls="") +ax.plot(m0_stop.lon, m0_stop.lat, marker=".", color="black", ls="") +ax.legend() ax.set_xlim(-10, 20), ax.set_ylim(-37, -21), ax.grid() -merging +m1 + # %% # Get spliting event @@ -212,10 +223,17 @@ # Display the position of the eddies before a splitting fig = plt.figure(figsize=(15, 8)) ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes") -spliting = close_to_i2.spliting_event() -spliting.display(ax) +close_to_i2.plot(ax, color_cycle=close_to_i2.COLORS) +s0, s1, s1_start = close_to_i2.spliting_event(triplet=True) +s0.display(ax, color="violet", lw=2, label="Eddies before splitting") +s1.display(ax, color="blueviolet", lw=2, label="Eddies after splitting") +s1_start.display(ax, color="black", lw=2, label="Eddies starting by splitting") +ax.plot(s0.lon, s0.lat, marker=".", color="purple", ls="") +ax.plot(s1.lon, s1.lat, marker=".", color="blueviolet", ls="") +ax.plot(s1_start.lon, s1_start.lat, marker=".", color="black", ls="") +ax.legend() ax.set_xlim(-10, 20), ax.set_ylim(-37, -21), ax.grid() -spliting +s1 # %% # Get birth event diff --git a/src/py_eddy_tracker/observations/network.py b/src/py_eddy_tracker/observations/network.py index 74bec98d..34a47c3a 100644 --- a/src/py_eddy_tracker/observations/network.py +++ b/src/py_eddy_tracker/observations/network.py @@ -113,8 +113,8 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1): """ Select network on time duration - :param int nb_day_min: Minimal number of day which must be covered by one network, if negative -> not used - :param int nb_day_max: Maximal number of day which must be covered by one network, if negative -> not used + :param int nb_day_min: Minimal number of day covered by one network, if negative -> not used + :param int nb_day_max: Maximal number of day covered by one network, if negative -> not used """ if nb_day_max < 0: nb_day_max = 1000000000000 @@ -132,7 +132,7 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1): @classmethod def from_split_network(cls, group_dataset, indexs, **kwargs): """ - Build a NetworkObservations object with Group dataset and indexs + Build a NetworkObservations object with Group dataset and indexes :param TrackEddiesObservations group_dataset: Group dataset :param indexs: result from split_network @@ -204,6 +204,9 @@ def __close_segment(cls, father, shift, connexions, distance): cls.__close_segment(son, shift, connexions, distance) def segment_relative_order(self, seg_origine): + """ + Compute the relative order of each segment to the chosen segment + """ i_s, i_e, i_ref = build_index(self.segment) segment_connexions = self.connexions() relative_tr = -ones(i_s.shape, dtype="i4") @@ -216,6 +219,9 @@ def segment_relative_order(self, seg_origine): return d def relative(self, i_obs, order=2, direct=True, only_past=False, only_future=False): + """ + Extract the segments at a certain order. + """ d = self.segment_relative_order(self.segment[i_obs]) m = (d <= order) * (d != -1) return self.extract_with_mask(m) @@ -234,7 +240,7 @@ def only_one_network(self): """ _, i_start, _ = self.index_network if len(i_start) > 1: - raise Exception("Several network") + raise Exception("Several networks") def position_filter(self, median_half_window, loess_half_window): self.median_filter(median_half_window, "time", "lon").loess_filter( @@ -266,7 +272,15 @@ def display_timeline( self, ax, event=True, field=None, method=None, factor=1, **kwargs ): """ - Must be call on only one network + Plot a timeline of a network. + Must be called on only one network. + + :param matplotlib.axes.Axes ax: matplotlib axe used to draw + :param bool event: if True, draw the splitting and merging events + :param str,array field: yaxis values, if None, segments are used + :param str method: if None, mean values are used + :param float factor: to multiply field + :return: plot mappable """ self.only_one_network() j = 0 @@ -516,6 +530,7 @@ def extract_event(self, indices): @property def segment_track_array(self): + """Return a unique segment id when multiple networks are considered""" return build_unique_array(self.segment, self.track) def birth_event(self): @@ -542,27 +557,65 @@ def death_event(self): indices.append(i.stop - 1) return self.extract_event(list(set(indices))) - def merging_event(self): - indices = list() + def merging_event(self, triplet=False): + """Return observation after a merging event. + + If `triplet=True` return the eddy after a merging event, the eddy before the merging event, + and the eddy stopped due to merging. + """ + idx_m1 = list() + if triplet: + idx_m0_stop = list() + idx_m0 = list() + for i, _, _ in self.iter_on(self.segment_track_array): nb = i.stop - i.start if nb == 0: continue i_n = self.next_obs[i.stop - 1] if i_n != -1: - indices.append(i.stop - 1) - return self.extract_event(list(set(indices))) + if triplet: + idx_m0_stop.append(i.stop - 1) + idx_m0.append(self.previous_obs[i_n]) + idx_m1.append(i_n) + + if triplet: + return ( + self.extract_event(list(idx_m1)), + self.extract_event(list(idx_m0)), + self.extract_event(list(idx_m0_stop)), + ) + else: + return self.extract_event(list(set(idx_m1))) - def spliting_event(self): - indices = list() + def spliting_event(self, triplet=False): + """Return observation before a splitting event. + + If `triplet=True` return the eddy before a splitting event, the eddy after the splitting event, + and the eddy starting due to splitting. + """ + idx_s0 = list() + if triplet: + idx_s1_start = list() + idx_s1 = list() for i, _, _ in self.iter_on(self.segment_track_array): nb = i.stop - i.start if nb == 0: continue i_p = self.previous_obs[i.start] if i_p != -1: - indices.append(i.start) - return self.extract_event(list(set(indices))) + if triplet: + idx_s1_start.append(i.start) + idx_s1.append(self.next_obs[i_p]) + idx_s0.append(i_p) + if triplet: + return ( + self.extract_event(list(idx_s0)), + self.extract_event(list(idx_s1)), + self.extract_event(list(idx_s1_start)), + ) + else: + return self.extract_event(list(set(idx_s0))) def dissociate_network(self): """ @@ -655,7 +708,7 @@ def remove_dead_end(self, nobs=3, recursive=0, mask=None): self.only_one_network() segments_keep = list() connexions = self.connexions() - for i, b0, b1 in self.iter_on("segment"): + for i, b0, _ in self.iter_on("segment"): nb = i.stop - i.start if mask and mask[i].any(): segments_keep.append(b0) @@ -852,6 +905,7 @@ def apply_replace(x, x0, x1): @njit(cache=True) def build_unique_array(id1, id2): + """Give a unique id for each (id1, id2) with id1 and id2 increasing monotonically""" k = 0 new_id = empty(id1.shape, dtype=id1.dtype) id1_previous = id1[0]