full example fit contour

Author: AntSimi

examples/14_generic_tools/pet_fit_contour.py

@@ -0,0 +1,61 @@
+"""
+Contour fit
+===========
+
+Two type of fit :
+ - Ellipse
+ - Circle
+
+In the two case we use a least square algorithm
+"""
+
+from matplotlib import pyplot as plt
+from numpy import cos, linspace, radians, sin
+
+from py_eddy_tracker import data
+from py_eddy_tracker.generic import coordinates_to_local, local_to_coordinates
+from py_eddy_tracker.observations.observation import EddiesObservations
+from py_eddy_tracker.poly import fit_circle_, fit_ellips
+
+# %%
+# Load example identification file
+a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
+
+# %%
+# Function to draw circle or ellips from parameter
+def build_circle(x0, y0, r):
+    angle = radians(linspace(0, 360, 50))
+    x_norm, y_norm = cos(angle), sin(angle)
+    return local_to_coordinates(x_norm * r, y_norm * r, x0, y0)
+
+
+def build_ellips(x0, y0, a, b, theta):
+    angle = radians(linspace(0, 360, 50))
+    x = a * cos(theta) * cos(angle) - b * sin(theta) * sin(angle)
+    y = a * sin(theta) * cos(angle) + b * cos(theta) * sin(angle)
+    return local_to_coordinates(x, y, x0, y0)
+
+
+# %%
+# Plot fitted circle or ellips on stored contour
+xs, ys = a.contour_lon_s, a.contour_lat_s
+
+fig = plt.figure(figsize=(15, 15))
+
+j = 1
+for i in range(0, 800, 30):
+    x, y = xs[i], ys[i]
+    x0_, y0_ = x.mean(), y.mean()
+    x_, y_ = coordinates_to_local(x, y, x0_, y0_)
+    ax = fig.add_subplot(4, 4, j)
+    ax.grid(), ax.set_aspect("equal")
+    ax.plot(x, y, label="store", color="black")
+    x0, y0, a, b, theta = fit_ellips(x_, y_)
+    x0, y0 = local_to_coordinates(x0, y0, x0_, y0_)
+    ax.plot(*build_ellips(x0, y0, a, b, theta), label="ellips", color="green")
+    x0, y0, radius, shape_error = fit_circle_(x_, y_)
+    x0, y0 = local_to_coordinates(x0, y0, x0_, y0_)
+    ax.plot(*build_circle(x0, y0, radius), label="circle", color="red", lw=0.5)
+    if j == 16:
+        break
+    j += 1

notebooks/python_module/14_generic_tools/pet_fit_contour.ipynb

@@ -0,0 +1,65 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Contour fit\n\nTwo type of fit :\n - Ellipse\n - Circle\n\nIn the two case we use a least square algorithm\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "from py_eddy_tracker.poly import fit_circle_, fit_ellips, fit_circle\nfrom py_eddy_tracker.generic import local_to_coordinates, coordinates_to_local\nfrom matplotlib import pyplot as plt\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.observations.observation import EddiesObservations\nfrom numpy import radians, linspace, cos, sin\n\na = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20190223.nc\"))\n\n\ndef build_circle(x0, y0, r):\n    angle = radians(linspace(0, 360, 50))\n    x_norm, y_norm = cos(angle), sin(angle)\n    return local_to_coordinates(x_norm * r, y_norm * r, x0, y0)\n\ndef build_ellips(x0, y0, a, b, theta):\n    angle = radians(linspace(0, 360, 50))\n    x = a * cos(theta) * cos(angle) - b * sin(theta) * sin(angle)\n    y = a * sin(theta) * cos(angle) + b * cos(theta) * sin(angle)\n    return local_to_coordinates(x, y, x0, y0)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "xs,ys=a.contour_lon_s,a.contour_lat_s\n\nfor i in range(20):\n    x, y = xs[i], ys[i]\n    x0_, y0_ = x.mean(), y.mean()\n    x_, y_ = coordinates_to_local(x,y, x0_, y0_)\n    fig = plt.figure()\n    ax = fig.add_subplot(111)\n    ax.grid(), ax.set_aspect('equal')\n    ax.plot(x, y, label='store')\n    x0, y0, a,b, theta = fit_ellips(x_,y_)\n    x0, y0 = local_to_coordinates(x0, y0, x0_, y0_)\n    ax.plot(*build_ellips(x0, y0, a,b, theta), label='ellips')\n    x0, y0, radius, shape_error = fit_circle_(x_,y_)\n    x0,y0 = local_to_coordinates(x0,y0, x0_, y0_)\n    ax.plot(*build_circle(x0,y0, radius), label='circle')\n    ax.legend()"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.7.7"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

src/py_eddy_tracker/poly.py

@@ -548,29 +548,12 @@ def fit_ellips(x, y):
     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)
+    AB1 = 2 * (a * e ** 2 + c * d ** 2 - b * d * e - det)
+    AB2 = a + c
+    AB3 = ((a - c) ** 2 + b ** 2) ** 0.5
+    A = -((AB1 * (AB2 + AB3)) ** 0.5) / det
+    B = -((AB1 * (AB2 - AB3)) ** 0.5) / det
+    theta = arctan((c - a - AB3) / b)
     return x0, y0, A, B, theta