first version to fit an ellips

Author: AntSimi

examples/02_eddy_identification/pet_eddy_detection.py

@@ -143,18 +143,15 @@ def update_axes(ax, mappable=None):
 # %%
 # Display the speed radius of the detected eddies
 ax = start_axes("Speed Radius (km)")
-a.scatter(ax, "radius_s", vmin=10, vmax=50, s=80, ref=-10, cmap="magma_r", factor=0.001)
-m = c.scatter(
-    ax, "radius_s", vmin=10, vmax=50, s=80, ref=-10, cmap="magma_r", factor=0.001
-)
+kwargs = dict(vmin=10, vmax=50, s=80, ref=-10, cmap="magma_r", factor=0.001)
+a.scatter(ax, "radius_s", **kwargs)
+m = c.scatter(ax, "radius_s", **kwargs)
 update_axes(ax, m)
 
 # %%
 # Filling the effective radius contours with the effective radius values
 ax = start_axes("Effective Radius (km)")
 kwargs = dict(vmin=10, vmax=80, cmap="magma_r", factor=0.001, lut=14, ref=-10)
 a.filled(ax, "effective_radius", **kwargs)
-m = c.filled(
-    ax, "radius_e", vmin=10, vmax=80, cmap="magma_r", factor=0.001, lut=14, ref=-10
-)
+m = c.filled(ax, "radius_e", **kwargs)
 update_axes(ax, m)

examples/14_generic_tools/pet_fit_contour.py

@@ -253,7 +253,7 @@
       },
       "outputs": [],
       "source": [
-        "ax = start_axes(\"Speed Radius (km)\")\na.scatter(ax, \"radius_s\", vmin=10, vmax=50, s=80, ref=-10, cmap=\"magma_r\", factor=0.001)\nm = c.scatter(\n    ax, \"radius_s\", vmin=10, vmax=50, s=80, ref=-10, cmap=\"magma_r\", factor=0.001\n)\nupdate_axes(ax, m)"
+        "ax = start_axes(\"Speed Radius (km)\")\nkwargs = dict(vmin=10, vmax=50, s=80, ref=-10, cmap=\"magma_r\", factor=0.001)\na.scatter(ax, \"radius_s\", **kwargs)\nm = c.scatter(\n    ax, \"radius_s\", **kwargs\n)\nupdate_axes(ax, m)"
       ]
     },
     {
@@ -271,7 +271,7 @@
       },
       "outputs": [],
       "source": [
-        "ax = start_axes(\"Effective Radius (km)\")\nkwargs = dict(vmin=10, vmax=80, cmap=\"magma_r\", factor=0.001, lut=14, ref=-10)\na.filled(ax, \"effective_radius\", **kwargs)\nm = c.filled(\n    ax, \"radius_e\", vmin=10, vmax=80, cmap=\"magma_r\", factor=0.001, lut=14, ref=-10\n)\nupdate_axes(ax, m)"
+        "ax = start_axes(\"Effective Radius (km)\")\nkwargs = dict(vmin=10, vmax=80, cmap=\"magma_r\", factor=0.001, lut=14, ref=-10)\na.filled(ax, \"effective_radius\", **kwargs)\nm = c.filled(\n    ax, \"radius_e\", **kwargs\n)\nupdate_axes(ax, m)"
       ]
     }
   ],

notebooks/python_module/02_eddy_identification/pet_eddy_detection.ipynb

@@ -619,7 +619,7 @@ def split_network(self, intern=True, **kwargs):
             dtype=[
                 ("group", self.tracks.dtype),
                 ("time", self.time.dtype),
-                ("track", "u2"),
+                ("track", "u4"),
                 ("previous_cost", "f4"),
                 ("next_cost", "f4"),
                 ("previous_obs", "i4"),

src/py_eddy_tracker/observations/tracking.py

@@ -7,7 +7,7 @@
 
 from numba import njit, prange
 from numba import types as numba_types
-from numpy import array, concatenate, empty, nan, ones, pi, where
+from numpy import arctan, array, concatenate, empty, nan, ones, pi, where
 from numpy.linalg import lstsq
 from Polygon import Polygon
 
@@ -471,6 +471,7 @@ def polygon_overlap(p0, p1, minimal_area=False):
     return cost
 
 
+# FIXME: only one function are needed
 @njit(cache=True)
 def fit_circle(x, y):
     """
@@ -517,6 +518,62 @@ def fit_circle(x, y):
     return x0, y0, radius, err
 
 
+@njit(cache=True)
+def fit_ellips(x, y):
+    r"""
+    From a polygon, function will fit an ellips.
+
+    Must be call with local coordinates (in m, to get a radius in m).
+
+    .. math:: (\frac{x - x_0}{a})^2 + (\frac{y - y_0}{b})^2 = 1
+
+    .. math:: (\frac{x^2 - 2 * x * x_0 + x_0 ^2}{a^2}) + \frac{y^2 - 2 * y * y_0 + y_0 ^2}{b^2}) = 1
+
+    In case of angle
+    https://en.wikipedia.org/wiki/Ellipse
+
+    """
+    # x,y = x[1:],y[1:]
+    nb = x.shape[0]
+    datas = ones((nb, 5), dtype=x.dtype)
+    datas[:, 0] = x ** 2
+    datas[:, 1] = x * y
+    datas[:, 2] = y ** 2
+    datas[:, 3] = x
+    datas[:, 4] = y
+    (a, b, c, d, e), _, _, _ = lstsq(datas, ones(nb, dtype=x.dtype))
+    det = b ** 2 - 4 * a * c
+    if det > 0:
+        print(det)
+    x0 = (2 * c * d - b * e) / det
+    y0 = (2 * a * e - b * d) / det
+
+    A = (
+        -(
+            (
+                2
+                * (a * e ** 2 + c * d ** 2 - b * d * e - det)
+                * (a + c + ((a - c) ** 2 + b ** 2) ** 0.5)
+            )
+            ** 0.5
+        )
+        / det
+    )
+    B = (
+        -(
+            (
+                2
+                * (a * e ** 2 + c * d ** 2 - b * d * e - det)
+                * (a + c - ((a - c) ** 2 + b ** 2) ** 0.5)
+            )
+            ** 0.5
+        )
+        / det
+    )
+    theta = arctan((c - a - ((a - c) ** 2 + b ** 2) ** 0.5) / b)
+    return x0, y0, A, B, theta
+
+
 @njit(cache=True)
 def fit_circle_(x, y):
     r"""