Change segment color in video and add comments (AntSimi#50)

- add scatter and event_map, change marker for event
- Add documentation in examples

Author: CoriPegliasco

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)

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,13 +85,15 @@ 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)
 
 # %%
 # 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"])

examples/16_network/pet_relative.py

@@ -19,15 +19,20 @@
 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))
 
 # %%
 # Timeline
 # --------
 
 # %%
-# 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)
@@ -39,27 +44,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()
@@ -69,6 +81,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])
@@ -86,6 +100,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)
@@ -95,6 +110,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)
@@ -106,7 +122,10 @@
 # %%
 # 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
 
 n_clean = n.remove_dead_end(nobs=10)
 fig = plt.figure(figsize=(15, 8))
@@ -120,6 +139,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])
@@ -130,6 +150,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(
@@ -141,18 +162,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)
@@ -162,6 +186,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)
@@ -173,6 +198,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()
@@ -183,6 +209,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()
@@ -193,6 +220,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()
@@ -203,6 +231,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()