Ioannou case

Figure 10 from https://doi.org/10.1002/2017JC013158

We want to find the Ierapetra Eddy described above in a network demonstration run.

from datetime import datetime, timedelta

import re

from matplotlib import colors, pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.ticker import FuncFormatter
from numpy import arange, array, pi, where

from py_eddy_tracker.appli.gui import Anim
from py_eddy_tracker.data import get_demo_path
from py_eddy_tracker.generic import coordinates_to_local
from py_eddy_tracker.gui import GUI_AXES
from py_eddy_tracker.observations.network import NetworkObservations
from py_eddy_tracker.poly import fit_ellipse

class VideoAnimation(FuncAnimation):
    def _repr_html_(self, *args, **kwargs):
        """To get video in html and have a player"""
        content = self.to_html5_video()
        return re.sub(
            r'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
        )

    def save(self, *args, **kwargs):
        if args[0].endswith("gif"):
            # In this case gif is used to create thumbnail which is not used but consume same time than video
            # So we create an empty file, to save time
            with open(args[0], "w") as _:
                pass
            return
        return super().save(*args, **kwargs)


@FuncFormatter
def formatter(x, pos):
    return (timedelta(x) + datetime(1950, 1, 1)).strftime("%d/%m/%Y")


def start_axes(title=""):
    fig = plt.figure(figsize=(13, 6))
    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94], projection=GUI_AXES)
    ax.set_xlim(19, 29), ax.set_ylim(31, 35.5)
    ax.set_aspect("equal")
    ax.set_title(title, weight="bold")
    return ax


def timeline_axes(title=""):
    fig = plt.figure(figsize=(15, 5))
    ax = fig.add_axes([0.03, 0.06, 0.90, 0.88])
    ax.set_title(title, weight="bold")
    ax.xaxis.set_major_formatter(formatter), ax.grid()
    return ax


def update_axes(ax, mappable=None):
    ax.grid(True)
    if mappable:
        return plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))

We know the network ID, we will get directly

ioannou_case = NetworkObservations.load_file(get_demo_path("network_med.nc")).network(
    651
)
print(ioannou_case.infos())

11948 obs 244 segments

It seems that this network is huge! Our case is visible at 22E 33.5N

ax = start_axes()
ioannou_case.plot(ax, color_cycle=ioannou_case.COLORS)
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()
_ = ioannou_case.display_timeline(ax)

Sub network and new numbering

Here we chose to keep only the order 3 segments relatives to our chosen eddy

i = where(
    (ioannou_case.lat > 33)
    * (ioannou_case.lat < 34)
    * (ioannou_case.lon > 22)
    * (ioannou_case.lon < 23)
    * (ioannou_case.time > 20630)
    * (ioannou_case.time < 20650)
)[0][0]
close_to_i3 = ioannou_case.relative(i, order=3)
close_to_i3.numbering_segment()

Anim

Quick movie to see better!

a = Anim(
    close_to_i3,
    figsize=(12, 4),
    cmap=colors.ListedColormap(
        list(close_to_i3.COLORS), name="from_list", N=close_to_i3.segment.max() + 1
    ),
    nb_step=7,
    dpi=70,
    field_color="segment",
    field_txt="segment",
)
a.ax.set_xlim(19, 30), a.ax.set_ylim(32, 35.25)
a.txt.set_position((21.5, 32.7))
# 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)

Classic display

ax = timeline_axes()
_ = close_to_i3.display_timeline(ax)

ax = start_axes("")
n_copy = close_to_i3.copy()
n_copy.position_filter(2, 4)
n_copy.plot(ax, color_cycle=n_copy.COLORS)
update_axes(ax)

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")
_ = n_copy.display_timeline(ax, field="lat", method="all")

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 [
    (datetime(i, 1, 1), datetime(i, 12, 31)) for i in (2004, 2005, 2006, 2007)
]:
    ref, delta = datetime(1950, 1, 1), 20
    b0_, b1_ = (b0 - ref).days, (b1 - ref).days
    ax = timeline_axes()
    ax.set_xlim(b0_ - delta, b1_ + delta)
    ax.set_ylim(10, 115)
    ax.axvline(b0_, color="k", lw=1.5, ls="--"), ax.axvline(
        b1_, color="k", lw=1.5, ls="--"
    )
    n_copy.display_timeline(
        ax, field="radius_e", method="all", lw=4, markersize=8, factor=1e-3
    )
    n_copy.display_timeline(
        ax, field="radius_s", method="all", lw=1, markersize=3, factor=1e-3
    )

• 
• 
• 
• 

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"])
cb.set_label("Effective shape error")

Rotation angle

For each obs, fit an ellipse to the contour, with theta the angle from the x-axis, a the semi ax in x direction and b the semi ax in y dimension

theta_ = list()
a_ = list()
b_ = list()
for obs in close_to_i3:
    x, y = obs["contour_lon_s"], obs["contour_lat_s"]
    x0_, y0_ = x.mean(), y.mean()
    x_, y_ = coordinates_to_local(x, y, x0_, y0_)
    x0, y0, a, b, theta = fit_ellipse(x_, y_)
    theta_.append(theta)
    a_.append(a)
    b_.append(b)
a_ = array(a_)
b_ = array(b_)

Theta

ax = timeline_axes()
m = close_to_i3.scatter_timeline(ax, theta_, vmin=-pi / 2, vmax=pi / 2, cmap="hsv")
_ = update_axes(ax, m["scatter"])

a

ax = timeline_axes()
m = close_to_i3.scatter_timeline(ax, a_ * 1e-3, vmin=0, vmax=80, cmap="Spectral_r")
_ = update_axes(ax, m["scatter"])

b

ax = timeline_axes()
m = close_to_i3.scatter_timeline(ax, b_ * 1e-3, vmin=0, vmax=80, cmap="Spectral_r")
_ = update_axes(ax, m["scatter"])

a/b

ax = timeline_axes()
m = close_to_i3.scatter_timeline(ax, a_ / b_, vmin=1, vmax=2, cmap="Spectral_r")
_ = update_axes(ax, m["scatter"])