Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 73 additions & 34 deletions examples/02_eddy_identification/pet_eddy_detection_ACC.py
Original file line number Diff line number Diff line change
@@ -1,35 +1,58 @@
"""
Eddy detection : Antartic circum polar
======================================
Eddy detection : Antartic Circumpolar Current
=============================================

Script will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)
This script detect eddies on the ADT field, and compute u,v with the method add_uv (use it only if the Equator is avoided)

Two ones with filtering adt and another without
Two detections are provided : with a filtered ADT and without filtering

"""
from datetime import datetime

from matplotlib import pyplot as plt
from matplotlib import style

from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset

pos_cb = [0.1, 0.52, 0.83, 0.015]
pos_cb2 = [0.1, 0.07, 0.4, 0.015]


def quad_axes(title):
fig = plt.figure(figsize=(13, 8.5))
fig.suptitle(title, weight="bold")
style.use("default")
fig = plt.figure(figsize=(13, 10))
fig.suptitle(title, weight="bold", fontsize=14)
axes = list()
for position in (
[0.05, 0.53, 0.44, 0.44],
[0.53, 0.53, 0.44, 0.44],
[0.05, 0.03, 0.44, 0.44],
[0.53, 0.03, 0.44, 0.44],
):

ax_pos = dict(
topleft=[0.1, 0.54, 0.4, 0.38],
topright=[0.53, 0.54, 0.4, 0.38],
botleft=[0.1, 0.09, 0.4, 0.38],
botright=[0.53, 0.09, 0.4, 0.38],
)

for key, position in ax_pos.items():
ax = fig.add_axes(position)
ax.set_xlim(5, 45), ax.set_ylim(-60, -37)
ax.set_aspect("equal"), ax.grid(True)
axes.append(ax)
return axes
if "right" in key:
ax.set_yticklabels("")
return fig, axes


def set_fancy_labels(fig, ticklabelsize=14, labelsize=14, labelweight="semibold"):
for ax in fig.get_axes():
ax.grid()
ax.grid(which="major", linestyle="-", linewidth="0.5", color="black")
if ax.get_ylabel() != "":
ax.set_ylabel(ax.get_ylabel(), fontsize=labelsize, fontweight=labelweight)
if ax.get_xlabel() != "":
ax.set_xlabel(ax.get_xlabel(), fontsize=labelsize, fontweight=labelweight)
if ax.get_title() != "":
ax.set_title(ax.get_title(), fontsize=labelsize, fontweight=labelweight)
ax.tick_params(labelsize=ticklabelsize)


# %%
Expand Down Expand Up @@ -69,20 +92,25 @@ def quad_axes(title):
# %%
# Figures
# -------
axs = quad_axes("General properties field")
m = g_raw.display(axs[0], "adt", vmin=-1, vmax=1, cmap="RdBu_r")
axs[0].set_title("ADT(m)")
m = g.display(axs[1], "adt_low", vmin=-1, vmax=1, cmap="RdBu_r")
axs[1].set_title("ADT (m) large scale with cut at 700 km")
m = g.display(axs[2], "adt", vmin=-1, vmax=1, cmap="RdBu_r")
axs[2].set_title("ADT (m) high scale with cut at 700 km")
cb = plt.colorbar(
m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation="horizontal"
)
cb.set_label("ADT(m)", labelpad=-2)
kw_adt = dict(vmin=-1.5, vmax=1.5, cmap=plt.get_cmap("RdBu_r", 30))
fig, axs = quad_axes("General properties field")
m = g_raw.display(axs[0], "adt", **kw_adt)
axs[0].set_title("Total ADT (m)")
m = g.display(axs[1], "adt_low", **kw_adt)
axs[1].set_title("ADT (m) large scale, cutoff at 700 km")
m2 = g.display(axs[2], "adt", cmap=plt.get_cmap("RdBu_r", 20), vmin=-0.5, vmax=0.5)
axs[2].set_title("ADT (m) high-pass filtered, a cutoff at 700 km")
cb = plt.colorbar(m, cax=axs[0].figure.add_axes(pos_cb), orientation="horizontal")
cb.set_label("ADT (m)", labelpad=0)
cb2 = plt.colorbar(m2, cax=axs[2].figure.add_axes(pos_cb2), orientation="horizontal")
cb2.set_label("ADT (m)", labelpad=0)
set_fancy_labels(fig)

# %%
# The large-scale North-South gradient is removed by the filtering step.

# %%
axs = quad_axes("")
fig, axs = quad_axes("")
axs[0].set_title("Without filter")
axs[0].set_ylabel("Contours used in eddies")
axs[1].set_title("With filter")
Expand All @@ -91,13 +119,18 @@ def quad_axes(title):
g.contours.display(axs[1], lw=0.5, only_used=True)
g_raw.contours.display(axs[2], lw=0.5, only_unused=True)
g.contours.display(axs[3], lw=0.5, only_unused=True)
set_fancy_labels(fig)

# %%
# Removing the large-scale North-South gradient reveals closed contours in the
# South-Western corner of the ewample region.

# %%
kw = dict(ref=-10, linewidth=0.75)
kw_a = dict(color="r", label="Anticyclonic ({nb_obs} eddies)")
kw_c = dict(color="b", label="Cyclonic ({nb_obs} eddies)")
kw_filled = dict(vmin=0, vmax=100, cmap="Spectral_r", lut=20, intern=True, factor=100)
axs = quad_axes("Comparison between two detection")
fig, axs = quad_axes("Comparison between two detections")
# Match with intern/inner contour
i_a, j_a, s_a = a_.match(a, intern=True, cmin=0.15)
i_c, j_c, s_c = c_.match(c, intern=True, cmin=0.15)
Expand All @@ -107,10 +140,8 @@ def quad_axes(title):
c_.index(i_c).filled(axs[0], s_c, **kw_filled)
m = c.index(j_c).filled(axs[1], s_c, **kw_filled)

cb = plt.colorbar(
m, cax=axs[0].figure.add_axes([0.03, 0.51, 0.94, 0.01]), orientation="horizontal"
)
cb.set_label("Similarity index", labelpad=-5)
cb = plt.colorbar(m, cax=axs[0].figure.add_axes(pos_cb), orientation="horizontal")
cb.set_label("Similarity index (%)", labelpad=-5)
a_.display(axs[0], **kw, **kw_a), c_.display(axs[0], **kw, **kw_c)
a.display(axs[1], **kw, **kw_a), c.display(axs[1], **kw, **kw_c)

Expand All @@ -124,6 +155,12 @@ def quad_axes(title):

for ax in axs:
ax.legend()

set_fancy_labels(fig)

# %%
# Very similar eddies have Similarity Indexes >= 40%

# %%
# Criteria for rejecting a contour :
# 0. Accepted (green)
Expand All @@ -140,10 +177,10 @@ def quad_axes(title):

for i, (label, field, factor, stop) in enumerate(
(
("speed radius (km)", "radius_s", 0.001, 120),
("outter radius (km)", "radius_e", 0.001, 120),
("amplitude (cm)", "amplitude", 100, 25),
("speed max (cm/s)", "speed_average", 100, 25),
("Speed radius (km)", "radius_s", 0.001, 120),
("Effective radius (km)", "radius_e", 0.001, 120),
("Amplitude (cm)", "amplitude", 100, 25),
("Speed max (cm/s)", "speed_average", 100, 25),
)
):
ax = fig.add_subplot(2, 2, i + 1, title=label)
Expand All @@ -166,3 +203,5 @@ def quad_axes(title):
ax.plot((0, 1000), (0, 1000), "g")
ax.set_xlim(0, stop), ax.set_ylim(0, stop)
ax.legend()

set_fancy_labels(fig)