update network example to be build remotely

Author: AntSimi

examples/01_general_things/pet_storage.py

@@ -9,8 +9,10 @@
 There are 3 class of files:
 
 - **Eddies collections** : contain a list of eddies without link between them
-- **Track eddies collections** : manage eddies associated in trajectories, the ```track``` field allows to separate each trajectory
-- **Network eddies collections** : manage eddies associated in networks, the ```track``` and ```segment``` fields allow to separate observations
+- **Track eddies collections** :
+   manage eddies associated in trajectories, the ```track``` field allows to separate each trajectory
+- **Network eddies collections** :
+   manage eddies associated in networks, the ```track``` and ```segment``` fields allow to separate observations
 """
 
 import py_eddy_tracker_sample
@@ -62,7 +64,8 @@
 # Tracks eddies collections add several fields :
 #
 # - **track** : Trajectory number
-# - **observation_flag** : Flag indicating if the value is interpolated between two observations or not (0: observed eddy, 1: interpolated eddy)"
+# - **observation_flag** : Flag indicating if the value is interpolated between two observations or not
+#                          (0: observed eddy, 1: interpolated eddy)"
 # - **observation_number** : Eddy temporal index in a trajectory, days starting at the eddy first detection
 # - **cost_association** : result of the cost function to associate the eddy with the next observation
 eddies_tracks = TrackEddiesObservations.load_file(
@@ -78,7 +81,8 @@
 #
 # - track :  ID of network (ID 0 correspond to lonely eddies)
 # - segment :  ID of a segment within a network (from 1 to N)
-# - previous_obs : Index of the previous observation in the full dataset, if -1 there are no previous observation (the segment starts)
+# - previous_obs : Index of the previous observation in the full dataset,
+#                  if -1 there are no previous observation (the segment starts)
 # - next_obs : Index of the next observation in the full dataset, if -1 there are no next observation (the segment ends)
 # - previous_cost : Result of the cost function (1 is a good association, 0 is bad) with previous observation
 # - next_cost : Result of the cost function (1 is a good association, 0 is bad) with next observation

examples/07_cube_manipulation/pet_fsle_med.py

@@ -70,7 +70,7 @@ def check_p(x, y, flse, theta, m_set, m, dt, dist_init=0.02, dist_max=0.6):
             theta[i] = arctan2(at1, at2 + s2) * 180 / pi
             # To know where value are set
             m_set[i] = False
-            # To stop partcile advection
+            # To stop particle advection
             m[i0], m[i_n], m[i_e] = True, True, True
 
 

examples/16_network/pet_atlas.py

@@ -6,16 +6,20 @@
 from numpy import ma
 
 import py_eddy_tracker.gui
-from py_eddy_tracker.data import get_path
+from py_eddy_tracker.data import get_remote_sample
 from py_eddy_tracker.observations.network import NetworkObservations
 
-n = NetworkObservations.load_file(get_path("Anticyclonic_seg.nc"))
+n = NetworkObservations.load_file(
+    get_remote_sample(
+        "eddies_med_adt_allsat_dt2018_err70_filt500_order1/Anticyclonic_network.nc"
+    )
+)
 # %%
 # Parameters
 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("viridis", 10))
+kw_ratio = dict(cmap=plt.get_cmap("YlGnBu_r", 10))
 
 
 # %%
@@ -57,19 +61,17 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m).set_label("Pixel used in % of time")
 
 # %%
-# Display the ratio between the short and total presence.
-#
-# Light = mostly short networks
-ax = start_axes("")
-m = g_10.display(
-    ax,
-    **kw_ratio,
-    vmin=50,
-    vmax=100,
-    name=g_10.vars["count"] * 100.0 / g_all.vars["count"]
-)
+# Ratio
+# ^^^^^
+# Ratio between the longer and total presence
+ax = start_axes("")
+g_ = g_10.vars["count"] * 100.0 / g_all.vars["count"]
+m = g_10.display(ax, **kw_ratio, vmin=50, vmax=100, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
+
 # %%
+# Blue = mostly short networks
+#
 # Network longer than 20 days
 # ---------------------------
 # Display the % of time each pixel (1/10°) is within an anticyclonic network
@@ -81,57 +83,34 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m).set_label("Pixel used in % of time")
 
 # %%
-# Display the ratio between the short and total presence.
-#
-# Light = mostly short networks
-ax = start_axes("")
-m = g_20.display(
-    ax,
-    **kw_ratio,
-    vmin=50,
-    vmax=100,
-    name=g_20.vars["count"] * 100.0 / g_all.vars["count"]
-)
+# Ratio
+# ^^^^^
+# Ratio between the longer and total presence
+ax = start_axes("")
+g_ = g_20.vars["count"] * 100.0 / g_all.vars["count"]
+m = g_20.display(ax, **kw_ratio, vmin=50, vmax=100, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
 # %%
-# Display the ratio between the short and total presence.
-#
-# Light = mostly short networks
-#
-# Networks shorter than 365 days are masked
-ax = start_axes("")
-m = g_20.display(
-    ax,
-    **kw_ratio,
-    vmin=50,
-    vmax=100,
-    name=ma.array(
-        g_20.vars["count"] * 100.0 / g_all.vars["count"], mask=g_all.vars["count"] < 365
-    )
+# Now we will hide pixel which are used less than 365 times
+g_ = ma.array(
+    g_20.vars["count"] * 100.0 / g_all.vars["count"], mask=g_all.vars["count"] < 365
 )
+ax = start_axes("")
+m = g_20.display(ax, **kw_ratio, vmin=50, vmax=100, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 # %%
-# Display the ratio between the short and total presence.
-#
-# Networks longer than 365 days are masked
-#
-# # -> Coastal areas are mostly populated by short networks
-ax = start_axes("")
-m = g_20.display(
-    ax,
-    **kw_ratio,
-    vmin=50,
-    vmax=100,
-    name=ma.array(
-        g_20.vars["count"] * 100.0 / g_all.vars["count"],
-        mask=g_all.vars["count"] >= 365,
-    )
+# Now we will hide pixel which are used more than 365 times
+ax = start_axes("")
+g_ = ma.array(
+    g_20.vars["count"] * 100.0 / g_all.vars["count"], mask=g_all.vars["count"] >= 365
 )
+m = g_20.display(ax, **kw_ratio, vmin=50, vmax=100, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
-
 # %%
+# Coastal areas are mostly populated by short networks
+#
 # All merging
 # -----------
 # Display the occurence of merging events
@@ -143,18 +122,13 @@ def update_axes(ax, mappable=None):
 # %%
 # Ratio merging events / eddy presence
 ax = start_axes("")
-m = g_all_merging.display(
-    ax,
-    **kw_ratio,
-    vmin=0,
-    vmax=5,
-    name=g_all_merging.vars["count"] * 100.0 / g_all.vars["count"]
-)
+g_ = g_all_merging.vars["count"] * 100.0 / g_all.vars["count"]
+m = g_all_merging.display(ax, **kw_ratio, vmin=0, vmax=5, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
 # %%
-# Merging in networks longer than 10 days, with dead end remove (shorter than 5 days)
-# -----------------------------------------------------------------------------------
+# Merging in networks longer than 10 days, with dead end remove (shorter than 10 observations)
+# --------------------------------------------------------------------------------------------
 ax = start_axes("")
 merger = n10.remove_dead_end(nobs=10).merging_event()
 g_10_merging = merger.grid_count(bins)
@@ -170,17 +144,13 @@ def update_axes(ax, mappable=None):
 m = g_10_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("")
-m = g_10_merging.display(
-    ax,
-    **kw_ratio,
-    vmin=0,
-    vmax=5,
-    name=ma.array(
-        g_10_merging.vars["count"] * 100.0 / g_10.vars["count"],
-        mask=g_10.vars["count"] < 365,
-    )
+g_ = ma.array(
+    g_10_merging.vars["count"] * 100.0 / g_10.vars["count"],
+    mask=g_10.vars["count"] < 365,
 )
+m = g_10_merging.display(ax, **kw_ratio, vmin=0, vmax=5, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
 # %%
@@ -193,15 +163,10 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m).set_label("Pixel used in % of time")
 
 # %%
-# Ratio merging events / eddy presence
+# Ratio spliting events / eddy presence
 ax = start_axes("")
-m = g_all_spliting.display(
-    ax,
-    **kw_ratio,
-    vmin=0,
-    vmax=5,
-    name=g_all_spliting.vars["count"] * 100.0 / g_all.vars["count"]
-)
+g_ = g_all_spliting.vars["count"] * 100.0 / g_all.vars["count"]
+m = g_all_spliting.display(ax, **kw_ratio, vmin=0, vmax=5, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 
 # %%
@@ -213,14 +178,9 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m).set_label("Pixel used in % of time")
 # %%
 ax = start_axes("")
-m = g_10_spliting.display(
-    ax,
-    **kw_ratio,
-    vmin=0,
-    vmax=5,
-    name=ma.array(
-        g_10_spliting.vars["count"] * 100.0 / g_10.vars["count"],
-        mask=g_10.vars["count"] < 365,
-    )
+g_ = ma.array(
+    g_10_spliting.vars["count"] * 100.0 / g_10.vars["count"],
+    mask=g_10.vars["count"] < 365,
 )
+m = g_10_spliting.display(ax, **kw_ratio, vmin=0, vmax=5, name=g_)
 update_axes(ax, m).set_label("Pixel used in % all atlas")

examples/16_network/pet_group_anim.py

@@ -24,7 +24,7 @@ def _repr_html_(self, *args, **kwargs):
         """To get video in html and have a player"""
         content = self.to_html5_video()
         return re.sub(
-            'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
+            r'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
         )
 
     def save(self, *args, **kwargs):
@@ -93,7 +93,7 @@ def get_group_array(self, results, nb_obs):
 
 # %%
 # Get data from period and area
-e = EddiesObservations.load_file(data.get_path("Anticyclonic_seg.nc"))
+e = EddiesObservations.load_file(data.get_path("network_med.nc"))
 e = e.extract_with_mask((e.time >= t0) * (e.time < t1)).extract_with_area(
     dict(llcrnrlon=25, urcrnrlon=35, llcrnrlat=31, urcrnrlat=37.5)
 )
@@ -134,7 +134,7 @@ def update(frame):
 # ----
 fig = plt.figure(figsize=(16, 9), dpi=50)
 ax = fig.add_axes([0, 0, 1, 1])
-ax.set_aspect("equal"), ax.grid(), ax.set_xlim(26, 34), ax.set_ylim(32, 36.5)
+ax.set_aspect("equal"), ax.grid(), ax.set_xlim(26, 34), ax.set_ylim(31, 35.5)
 groups = ax.scatter(
     e.lon,
     e.lat,
@@ -145,6 +145,6 @@ def update(frame):
 )
 current_contour = ax.plot([], [], "k", lw=2, label="Current contour")[0]
 matched_contour = ax.plot([], [], "r", lw=1, ls="--", label="Candidate contour")[0]
-txt = ax.text(31.5, 35.5, "", fontsize=25)
+txt = ax.text(29, 35, "", fontsize=25)
 ax.legend(fontsize=25)
 ani = VideoAnimation(fig, update, frames=len(NETWORK_GROUPS), interval=220)

examples/16_network/pet_ioannou_2017_case.py

@@ -3,25 +3,22 @@
 ============
 Figure 10 from https://doi.org/10.1002/2017JC013158
 
+We want to find the Ierapetra Eddy described above in a network demonstration run.
 """
 
 # %%
-# We want to find the Ierapetra Eddy described above in the networks
-
 import re
-
-# %%
 from datetime import datetime, timedelta
 
-import numpy as np
 from matplotlib import colors
 from matplotlib import pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.ticker import FuncFormatter
+from numpy import arange, where
 
 import py_eddy_tracker.gui
 from py_eddy_tracker.appli.gui import Anim
-from py_eddy_tracker.data import get_remote_sample
+from py_eddy_tracker.data import get_path
 from py_eddy_tracker.observations.network import NetworkObservations
 
 
@@ -74,23 +71,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(
-    get_remote_sample(
-        "eddies_med_adt_allsat_dt2018_err70_filt500_order1/Anticyclonic_network.nc"
-    )
-)
-i = np.where(
-    (n.lat > 33)
-    * (n.lat < 34)
-    * (n.lon > 22)
-    * (n.lon < 23)
-    * (n.time > 20630)
-    * (n.time < 20650)
-)[0][0]
-ioannou_case = n.extract_with_mask(n.track == n.track[i])
+# We know the network ID, we will get directly
+ioannou_case = NetworkObservations.load_file(get_path("network_med.nc")).network(651)
 print(ioannou_case.infos())
 
 # %%
@@ -113,7 +95,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(
+i = where(
     (ioannou_case.lat > 33)
     * (ioannou_case.lat < 34)
     * (ioannou_case.lon > 22)
@@ -128,22 +110,22 @@ 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=cmap,
+    cmap=colors.ListedColormap(
+        list(close_to_i3.COLORS), name="from_list", N=close_to_i3.segment.max() + 1
+    ),
     nb_step=7,
-    dpi=80,
+    dpi=70,
     field_color="segment",
     field_txt="segment",
 )
 a.ax.set_xlim(19, 30), a.ax.set_ylim(32, 35.25)
 a.ax.update_env()
 a.txt.set_position((21.5, 32.7))
-kwargs = dict(frames=np.arange(*a.period), interval=100)
+# We display in video only from the 100th day to the 500th
+kwargs = dict(frames=arange(*a.period)[100:501], interval=100)
 ani = VideoAnimation(a.fig, a.func_animation, **kwargs)
 
 # %%