Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
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
Next Next commit
- add ellipse fit
- minor english modifs
  • Loading branch information
CoriPegliasco committed Apr 1, 2021
commit d81604d8bbedd833cf4a108028ea02033c6c9779
44 changes: 44 additions & 0 deletions examples/16_network/pet_ioannou_2017_case.py
Original file line number Diff line number Diff line change
Expand Up @@ -188,3 +188,47 @@ def update_axes(ax, mappable=None):
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
# --------------
from py_eddy_tracker.generic import coordinates_to_local, local_to_coordinates
from py_eddy_tracker.poly import fit_ellips
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_ellips(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')
cb = 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')
cb = 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')
cb = 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')
cb = update_axes(ax, m["scatter"])

13 changes: 9 additions & 4 deletions src/py_eddy_tracker/dataset/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -609,17 +609,17 @@ def eddy_identification(
:param str vname: Grid name of v speed component
:param datetime.datetime date: Date which will be stored in object to date data
:param float,int step: Height between two layers in m
:param float,int shape_error: Maximal error allowed for outter contour in %
:param float,int shape_error: Maximal error allowed for outermost contour in %
:param int sampling: Number of points to store contours and speed profile
:param str sampling_method: Method to resample 'uniform' or 'visvalingam'
:param (int,int),None pixel_limit:
Min and max number of pixels inside the inner and the outer contour to be considered as an eddy
Min and max number of pixels inside the inner and the outermost contour to be considered as an eddy
:param float,None precision: Truncate values at the defined precision in m
:param str force_height_unit: Unit used for height unit
:param str force_speed_unit: Unit used for speed unit
:param dict kwargs: Argument given to amplitude

:return: Return a list of 2 elements: Anticyclone and Cyclone
:return: Return a list of 2 elements: Anticyclones and Cyclones
:rtype: py_eddy_tracker.observations.observation.EddiesObservations

.. minigallery:: py_eddy_tracker.GridDataset.eddy_identification
Expand Down Expand Up @@ -729,7 +729,8 @@ def eddy_identification(
for contour in contour_paths:
if contour.used:
continue
# FIXME : center could be not in contour and fit on raw sampling
# FIXME : center could be outside the contour due to the fit
# FIXME : warning : the fit is made on raw sampling
_, _, _, aerr = contour.fit_circle()

# Filter for shape
Expand All @@ -752,6 +753,7 @@ def eddy_identification(
):
continue

# Test the number of pixels within the outermost contour
# FIXME : Maybe limit max must be replace with a maximum of surface
if (
contour.nb_pixel < pixel_limit[0]
Expand All @@ -760,6 +762,9 @@ def eddy_identification(
contour.reject = 3
continue

# Here the considered contour passed shape_error test, masked_pixels test,
# values strictly above (AEs) or below (CEs) the contour, number_pixels test)

# Compute amplitude
reset_centroid, amp = self.get_amplitude(
contour,
Expand Down
26 changes: 14 additions & 12 deletions src/py_eddy_tracker/eddy_feature.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
class Amplitude(object):
"""
Class to calculate *amplitude* and counts of *local maxima/minima*
within a closed region of a sea level anomaly field.
within a closed region of a sea surface height field.
"""

EPSILON = 1e-8
Expand Down Expand Up @@ -66,8 +66,8 @@ def __init__(
:param array data:
:param float interval:
:param int mle: maximum number of local maxima in contour
:param int nb_step_min: number of interval to consider like an eddy
:param int nb_step_to_be_mle: number of interval to be consider like another maxima
:param int nb_step_min: number of intervals to consider an eddy
:param int nb_step_to_be_mle: number of intervals to be considered as an another maxima
"""

# Height of the contour
Expand Down Expand Up @@ -102,8 +102,9 @@ def __init__(
self.nb_pixel = i_x.shape[0]

# Only pixel in contour
# FIXME : change sla by ssh as the grid can be adt?
self.sla = data[contour.pixels_index]
# Amplitude which will be provide
# Amplitude which will be provided
self.amplitude = 0
# Maximum local extrema accepted
self.mle = mle
Expand All @@ -115,9 +116,10 @@ def within_amplitude_limits(self):
def all_pixels_below_h0(self, level):
"""
Check CSS11 criterion 1: The SSH values of all of the pixels
are below a given SSH threshold for cyclonic eddies.
are below (above) a given SSH threshold for cyclonic (anticyclonic)
eddies.
"""
# In some case pixel value must be very near of contour bounds
# In some cases pixel value may be very close to the contour bounds
if self.sla.mask.any() or ((self.sla.data - self.h_0) > self.EPSILON).any():
return False
else:
Expand Down Expand Up @@ -293,10 +295,10 @@ class Contours(object):

Attributes:
contour:
A matplotlib contour object of high-pass filtered SLA
A matplotlib contour object of high-pass filtered SSH

eddy:
A tracklist object holding the SLA data
A tracklist object holding the SSH data

grd:
A grid object
Expand Down Expand Up @@ -406,7 +408,7 @@ def __init__(self, x, y, z, levels, wrap_x=False, keep_unclose=False):
fig = Figure()
ax = fig.add_subplot(111)
if wrap_x:
logger.debug("wrapping activate to compute contour")
logger.debug("wrapping activated to compute contour")
x = concatenate((x, x[:1] + 360))
z = ma.concatenate((z, z[:1]))
logger.debug("X shape : %s", x.shape)
Expand Down Expand Up @@ -602,8 +604,8 @@ def display(
Must be 'shape_error', 'x', 'y' or 'radius'.
If define display_criterion is not use.
bins argument must be define
:param array bins: bins use to colorize contour
:param str cmap: Name of cmap to use for field display
:param array bins: bins used to colorize contour
:param str cmap: Name of cmap for field display
:param dict kwargs: look at :py:meth:`matplotlib.collections.LineCollection`

.. minigallery:: py_eddy_tracker.Contours.display
Expand Down Expand Up @@ -688,7 +690,7 @@ def display(
ax.autoscale_view()

def label_contour_unused_which_contain_eddies(self, eddies):
"""Select contour which contain several eddies"""
"""Select contour containing several eddies"""
if eddies.sign_type == 1:
# anticyclonic
sl = slice(None, -1)
Expand Down
10 changes: 5 additions & 5 deletions src/py_eddy_tracker/poly.py
Original file line number Diff line number Diff line change
Expand Up @@ -474,13 +474,13 @@ def polygon_overlap(p0, p1, minimal_area=False):
return cost


# FIXME: only one function are needed
# FIXME: only one function is needed
@njit(cache=True)
def fit_circle(x, y):
"""
From a polygon, function will fit a circle.

Must be call with local coordinates (in m, to get a radius in m).
Must be called with local coordinates (in m, to get a radius in m).

:param array x: x of polygon
:param array y: y of polygon
Expand Down Expand Up @@ -510,11 +510,11 @@ def fit_circle(x, y):
radius **= 0.5
x0 *= scale
y0 *= scale
# radius of fitted circle
# radius of fit circle
radius *= scale
# center X-position of fitted circle
# center X-position of fit circle
x0 += x_mean
# center Y-position of fitted circle
# center Y-position of fit circle
y0 += y_mean

err = shape_error(x, y, x0, y0, radius)
Expand Down