Skip to content

Commit c483827

Browse files
AntSimiAntSimi
committed
example of visvalingam use
1 parent 50bc33a commit c483827

File tree

6 files changed

+209
-8
lines changed

6 files changed

+209
-8
lines changed

‎examples/14_generic_tools/README.rst‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Polygon tools
2+
=============
3+

‎examples/14_generic_tools/pet_visvalingam.py‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
"""
2+
Visvalingam algorithm
3+
=====================
4+
"""
5+
import matplotlib.animation as animation
6+
from matplotlib import pyplot as plt
7+
from numba import njit
8+
from numpy import array, empty
9+
10+
from py_eddy_tracker import data
11+
from py_eddy_tracker.generic import uniform_resample
12+
from py_eddy_tracker.observations.observation import EddiesObservations
13+
from py_eddy_tracker.poly import vertice_overlap, visvalingam
14+
15+
16+
@njit(cache=True)
17+
def visvalingam_polys(x, y, nb_pt):
18+
nb = x.shape[0]
19+
x_new = empty((nb, nb_pt), dtype=x.dtype)
20+
y_new = empty((nb, nb_pt), dtype=y.dtype)
21+
for i in range(nb):
22+
x_new[i], y_new[i] = visvalingam(x[i], y[i], nb_pt)
23+
return x_new, y_new
24+
25+
26+
@njit(cache=True)
27+
def uniform_resample_polys(x, y, nb_pt):
28+
nb = x.shape[0]
29+
x_new = empty((nb, nb_pt), dtype=x.dtype)
30+
y_new = empty((nb, nb_pt), dtype=y.dtype)
31+
for i in range(nb):
32+
x_new[i], y_new[i] = uniform_resample(x[i], y[i], fixed_size=nb_pt)
33+
return x_new, y_new
34+
35+
36+
def update_line(num):
37+
nb = 50 - num - 20
38+
x_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb)
39+
for i, (x_, y_) in enumerate(zip(x_v, y_v)):
40+
lines_v[i].set_data(x_, y_)
41+
x_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb)
42+
for i, (x_, y_) in enumerate(zip(x_u, y_u)):
43+
lines_u[i].set_data(x_, y_)
44+
scores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0
45+
scores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0
46+
for i, (s_v, s_u) in enumerate(zip(scores_v, scores_u)):
47+
texts[i].set_text(f"Score uniform {s_u:.1f} %\nScore visvalingam {s_v:.1f} %")
48+
title.set_text(f"{nb} points by contour in place of 50")
49+
return (title, *lines_u, *lines_v, *texts)
50+
51+
52+
# %%
53+
# Load detection files
54+
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
55+
a = a.extract_with_mask((abs(a.lat) < 66) * (abs(a.radius_e) > 80e3))
56+
57+
nb_pt = 10
58+
x_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)
59+
x_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)
60+
scores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0
61+
scores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0
62+
d_6 = scores_v - scores_u
63+
nb_pt = 18
64+
x_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)
65+
x_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)
66+
scores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0
67+
scores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0
68+
d_12 = scores_v - scores_u
69+
a = a.index(array((d_6.argmin(), d_6.argmax(), d_12.argmin(), d_12.argmax())))
70+
71+
72+
# %%
73+
fig = plt.figure()
74+
axs = [
75+
fig.add_subplot(221),
76+
fig.add_subplot(222),
77+
fig.add_subplot(223),
78+
fig.add_subplot(224),
79+
]
80+
lines_u, lines_v, texts, score_text = list(), list(), list(), list()
81+
for i, obs in enumerate(a):
82+
axs[i].set_aspect("equal")
83+
axs[i].grid()
84+
axs[i].set_xticklabels([]), axs[i].set_yticklabels([])
85+
axs[i].plot(
86+
obs["contour_lon_e"], obs["contour_lat_e"], "r", lw=6, label="Original contour"
87+
)
88+
lines_v.append(axs[i].plot([], [], color="limegreen", lw=4, label="visvalingam")[0])
89+
lines_u.append(
90+
axs[i].plot([], [], color="black", lw=2, label="uniform resampling")[0]
91+
)
92+
texts.append(axs[i].set_title("", fontsize=8))
93+
axs[0].legend(fontsize=8)
94+
title = fig.suptitle("")
95+
ani = animation.FuncAnimation(fig, update_line, 27)

‎notebooks/python_module/14_generic_tools/pet_visvalingam.ipynb‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {
7+
"collapsed": false
8+
},
9+
"outputs": [],
10+
"source": [
11+
"%matplotlib inline"
12+
]
13+
},
14+
{
15+
"cell_type": "markdown",
16+
"metadata": {},
17+
"source": [
18+
"\n# Visvalingam algorithm\n"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {
25+
"collapsed": false
26+
},
27+
"outputs": [],
28+
"source": [
29+
"import matplotlib.animation as animation\nfrom matplotlib import pyplot as plt\nfrom numba import njit\nfrom numpy import array, empty\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.generic import uniform_resample\nfrom py_eddy_tracker.observations.observation import EddiesObservations\nfrom py_eddy_tracker.poly import vertice_overlap, visvalingam\n\n\n@njit(cache=True)\ndef visvalingam_polys(x, y, nb_pt):\n nb = x.shape[0]\n x_new = empty((nb, nb_pt), dtype=x.dtype)\n y_new = empty((nb, nb_pt), dtype=y.dtype)\n for i in range(nb):\n x_new[i], y_new[i] = visvalingam(x[i], y[i], nb_pt)\n return x_new, y_new\n\n\n@njit(cache=True)\ndef uniform_resample_polys(x, y, nb_pt):\n nb = x.shape[0]\n x_new = empty((nb, nb_pt), dtype=x.dtype)\n y_new = empty((nb, nb_pt), dtype=y.dtype)\n for i in range(nb):\n x_new[i], y_new[i] = uniform_resample(x[i], y[i], fixed_size=nb_pt)\n return x_new, y_new\n\n\ndef update_line(num):\n nb = 50 - num - 20\n x_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb)\n for i, (x_, y_) in enumerate(zip(x_v, y_v)):\n lines_v[i].set_data(x_, y_)\n x_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb)\n for i, (x_, y_) in enumerate(zip(x_u, y_u)):\n lines_u[i].set_data(x_, y_)\n scores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0\n scores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0\n for i, (s_v, s_u) in enumerate(zip(scores_v, scores_u)):\n texts[i].set_text(f\"Score uniform {s_u:.1f} %\\nScore visvalingam {s_v:.1f} %\")\n title.set_text(f\"{nb} points by contour in place of 50\")\n return (title, *lines_u, *lines_v, *texts)"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"Load detection files\n\n"
37+
]
38+
},
39+
{
40+
"cell_type": "code",
41+
"execution_count": null,
42+
"metadata": {
43+
"collapsed": false
44+
},
45+
"outputs": [],
46+
"source": [
47+
"a = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20190223.nc\"))\na = a.extract_with_mask((abs(a.lat) < 66) * (abs(a.radius_e) > 80e3))\n\nnb_pt = 10\nx_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)\nx_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)\nscores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0\nscores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0\nd_6 = scores_v - scores_u\nnb_pt = 18\nx_v, y_v = visvalingam_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)\nx_u, y_u = uniform_resample_polys(a.contour_lon_e, a.contour_lat_e, nb_pt)\nscores_v = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_v, y_v) * 100.0\nscores_u = vertice_overlap(a.contour_lon_e, a.contour_lat_e, x_u, y_u) * 100.0\nd_12 = scores_v - scores_u\na = a.index(array((d_6.argmin(), d_6.argmax(), d_12.argmin(), d_12.argmax())))"
48+
]
49+
},
50+
{
51+
"cell_type": "code",
52+
"execution_count": null,
53+
"metadata": {
54+
"collapsed": false
55+
},
56+
"outputs": [],
57+
"source": [
58+
"fig = plt.figure()\naxs = [\n fig.add_subplot(221),\n fig.add_subplot(222),\n fig.add_subplot(223),\n fig.add_subplot(224),\n]\nlines_u, lines_v, texts, score_text = list(), list(), list(), list()\nfor i, obs in enumerate(a):\n axs[i].set_aspect(\"equal\")\n axs[i].grid()\n axs[i].set_xticklabels([]), axs[i].set_yticklabels([])\n axs[i].plot(\n obs[\"contour_lon_e\"], obs[\"contour_lat_e\"], \"r\", lw=6, label=\"Original contour\"\n )\n lines_v.append(axs[i].plot([], [], color=\"limegreen\", lw=4, label=\"visvalingam\")[0])\n lines_u.append(\n axs[i].plot([], [], color=\"black\", lw=2, label=\"uniform resampling\")[0]\n )\n texts.append(axs[i].set_title(\"\", fontsize=8))\naxs[0].legend(fontsize=8)\ntitle = fig.suptitle(\"\")\nani = animation.FuncAnimation(fig, update_line, 27)"
59+
]
60+
}
61+
],
62+
"metadata": {
63+
"kernelspec": {
64+
"display_name": "Python 3",
65+
"language": "python",
66+
"name": "python3"
67+
},
68+
"language_info": {
69+
"codemirror_mode": {
70+
"name": "ipython",
71+
"version": 3
72+
},
73+
"file_extension": ".py",
74+
"mimetype": "text/x-python",
75+
"name": "python",
76+
"nbconvert_exporter": "python",
77+
"pygments_lexer": "ipython3",
78+
"version": "3.7.7"
79+
}
80+
},
81+
"nbformat": 4,
82+
"nbformat_minor": 0
83+
}

‎src/py_eddy_tracker/appli/grid.py‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,19 @@ def eddy_id(args=None):
7979
parser.add_argument("--height_unit", default=None, help="Force height unit")
8080
parser.add_argument("--speed_unit", default=None, help="Force speed unit")
8181
parser.add_argument("--unregular", action="store_true", help="if grid is unregular")
82+
parser.add_argument(
83+
"--sampling",
84+
default=50,
85+
type=int,
86+
help="Array size used to build contour, first and last point will be the same",
87+
)
88+
parser.add_argument(
89+
"--sampling_method",
90+
default="visvalingam",
91+
type=str,
92+
choices=("visvalingam", "uniform"),
93+
help="Method to resample contour",
94+
)
8295
help = "Output will be wrote in zarr"
8396
parser.add_argument("--zarr", action="store_true", help=help)
8497
help = "Indexs to select grid : --indexs time=2, will select third step along time dimensions"
@@ -109,6 +122,8 @@ def eddy_id(args=None):
109122
cut_wavelength=args.cut_wavelength,
110123
filter_order=args.filter_order,
111124
indexs=args.indexs,
125+
sampling=args.sampling,
126+
sampling_method=args.sampling_method,
112127
**kwargs,
113128
)
114129
out_name = date.strftime("%(path)s/%(sign_type)s_%Y%m%d.nc")

‎src/py_eddy_tracker/dataset/grid.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@
6868
get_pixel_in_regular,
6969
poly_area,
7070
poly_contain_poly,
71+
visvalingam,
7172
winding_number_poly,
7273
)
7374

@@ -602,6 +603,7 @@ def eddy_identification(
602603
step=0.005,
603604
shape_error=55,
604605
sampling=50,
606+
sampling_method="visvalingam",
605607
pixel_limit=None,
606608
precision=None,
607609
force_height_unit=None,
@@ -618,6 +620,7 @@ def eddy_identification(
618620
:param float,int step: Height between two layers in m
619621
:param float,int shape_error: Maximal error allowed for outter contour in %
620622
:param int sampling: Number of points to store contours and speed profile
623+
:param str sampling_method: Method to resample 'uniform' or 'visvalingam'
621624
:param (int,int),None pixel_limit:
622625
Min and max number of pixels inside the inner and the outer contour to be considered as an eddy
623626
:param float,None precision: Truncate values at the defined precision in m
@@ -693,6 +696,7 @@ def eddy_identification(
693696
self.contours = Contours(x, y, data, levels, wrap_x=self.is_circular())
694697

695698
out_sampling = dict(fixed_size=sampling)
699+
resample = visvalingam if sampling_method == "visvalingam" else uniform_resample
696700
track_extra_variables = [
697701
"height_max_speed_contour",
698702
"height_external_contour",
@@ -832,10 +836,10 @@ def eddy_identification(
832836
obs.amplitude[:] = amp.amplitude
833837
obs.speed_average[:] = max_average_speed
834838
obs.num_point_e[:] = contour.lon.shape[0]
835-
xy_e = uniform_resample(contour.lon, contour.lat, **out_sampling)
839+
xy_e = resample(contour.lon, contour.lat, **out_sampling)
836840
obs.contour_lon_e[:], obs.contour_lat_e[:] = xy_e
837841
obs.num_point_s[:] = speed_contour.lon.shape[0]
838-
xy_s = uniform_resample(
842+
xy_s = resample(
839843
speed_contour.lon, speed_contour.lat, **out_sampling
840844
)
841845
obs.contour_lon_s[:], obs.contour_lat_s[:] = xy_s

‎src/py_eddy_tracker/poly.py‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -671,12 +671,12 @@ def tri_area2(x, y, i0, i1, i2):
671671

672672

673673
@njit(cache=True)
674-
def visvalingam(x, y, nb_pt=18):
674+
def visvalingam(x, y, fixed_size=18):
675675
"""Polygon simplification with visvalingam algorithm
676676
677677
:param array x:
678678
:param array y:
679-
:param int nb_pt: array size of out
679+
:param int fixed_size: array size of out
680680
:return: New (x, y) array
681681
:rtype: array,array
682682
"""
@@ -696,7 +696,7 @@ def visvalingam(x, y, nb_pt=18):
696696
i0 = i1
697697
i1 = i
698698
# we continue until we are equal to nb_pt
699-
while len(h) >= nb_pt:
699+
while len(h) >= fixed_size:
700700
# We pop lower area
701701
_, (i0, i1, i2) = heapq.heappop(h)
702702
# We check if triangle is valid(i0 or i2 not removed)
@@ -714,13 +714,14 @@ def visvalingam(x, y, nb_pt=18):
714714
# in this case we replace two point
715715
i0, i2 = i_p, i_n
716716
heapq.heappush(h, (tri_area2(x, y, i0, i1, i2), (i0, i1, i2)))
717-
x_new, y_new = empty(nb_pt, dtype=x.dtype), empty(nb_pt, dtype=y.dtype)
717+
x_new, y_new = empty(fixed_size, dtype=x.dtype), empty(fixed_size, dtype=y.dtype)
718718
j = 0
719719
for i, i_n in enumerate(i_next):
720720
if i_n == -1:
721721
x_new[j] = x[i]
722722
y_new[j] = y[i]
723723
j += 1
724-
x_new[j] = x_new[0]
725-
y_new[j] = y_new[0]
724+
# we copy first value to fill array end
725+
x_new[j:] = x_new[0]
726+
y_new[j:] = y_new[0]
726727
return x_new, y_new

0 commit comments

Comments
 (0)