diff --git a/Analysis/distribution_function.py b/Analysis/distribution_function.py
new file mode 100755
index 0000000000000000000000000000000000000000..15b6fce505bcff32572df2ff431e18211e1e176e
--- /dev/null
+++ b/Analysis/distribution_function.py
@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+
+import shapefile
+import numpy as np
+import matplotlib.pyplot as plt
+
+def plot_distr_funct(X, xlabel):
+ plt.figure()
+ plt.hist(X, bins="auto")
+ plt.xlabel(xlabel)
+ plt.ylabel("number of occurrences")
+
+ X = np.array(X)
+ X.sort()
+ nx = np.size(X)
+
+ print(xlabel + ":")
+ print("minimum value:", X[0])
+ print("maximum value:", X[-1])
+ print("number of values:", nx, "\n")
+
+ plt.figure()
+ plt.plot(X, (1 + np.arange(nx)) / nx)
+ plt.xlabel(xlabel)
+ plt.ylabel("distribution function")
+
+def get_var(end_filename):
+ with shapefile.Reader("extremum" + end_filename) as extremum, \
+ shapefile.Reader("outermost_contour" + end_filename) as outermost_cont:
+ speed = []
+ r_outer = []
+
+ for rec_extr, rec_outer in zip(extremum.iterRecords(),
+ outermost_cont.iterRecords()):
+ if rec_extr.speed != 1e4: speed.append(rec_extr.speed)
+ if rec_extr.valid == 1: r_outer.append(rec_outer.r_eq_area)
+
+ with shapefile.Reader("max_speed_contour" + end_filename) \
+ as max_speed_contour:
+ r_speed = []
+
+ for rec in max_speed_contour.iterRecords():
+ if rec.r_eq_area != - 100: r_speed.append(rec.r_eq_area)
+
+ return speed, r_outer, r_speed
+
+if __name__ == "__main__":
+ import argparse
+
+ parser = argparse.ArgumentParser()
+ parser.add_argument("-e", "--end_filename",
+ help = "add character string after "
+ "extremum, outermost contour and max_speed_contour",
+ default = "")
+ args = parser.parse_args()
+
+ speed, r_outer, r_speed = get_var(args.end_filename)
+
+ plot_distr_funct(speed, "speed, in m s-1")
+ plt.axhline(0.5, linestyle = "--", color = "black")
+ plt.axvline(linestyle = "--", color = "black")
+
+ plot_distr_funct(r_outer, "equivalent radius of outermost contour, in km")
+ plot_distr_funct(r_speed, "equivalent radius of max-speed contour, in km")
+
+ plt.show()
diff --git a/Tests/examine_ssh_values.py b/Analysis/examine_ssh_values.py
similarity index 100%
rename from Tests/examine_ssh_values.py
rename to Analysis/examine_ssh_values.py
diff --git a/Analysis/filter.py b/Analysis/filter.py
new file mode 100755
index 0000000000000000000000000000000000000000..4ed0faa0a8bd07f44cf1546eb895b00f88b25aa2
--- /dev/null
+++ b/Analysis/filter.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+
+import shapefile
+
+with shapefile.Reader("extremum_1") as extremum, \
+ shapefile.Reader("outermost_contour_1") as outermost_cont, \
+ shapefile.Reader("max_speed_contour_1") as max_speed_cont, \
+ shapefile.Writer("extremum_1_filtered") as extremum_filt, \
+ shapefile.Writer("outermost_contour_1_filtered") as outermost_cont_filt, \
+ shapefile.Writer("max_speed_contour_1_filtered") as max_speed_cont_filt:
+ extremum_filt.fields = extremum.fields[1:]
+ outermost_cont_filt.fields = outermost_cont.fields[1:]
+ max_speed_cont_filt.fields = max_speed_cont.fields[1:]
+
+ for shape_rec_extr, shape_rec_outer, shape_rec_max \
+ in zip(extremum, outermost_cont, max_speed_cont):
+ if shape_rec_extr.record.speed == 1e4 \
+ and shape_rec_outer.record.radius4 >= 2:
+ extremum_filt.record(*shape_rec_extr.record)
+ extremum_filt.shape(shape_rec_extr.shape)
+
+ outermost_cont_filt.record(*shape_rec_outer.record)
+ outermost_cont_filt.shape(shape_rec_outer.shape)
+
+ max_speed_cont_filt.record(*shape_rec_max.record)
+ max_speed_cont_filt.shape(shape_rec_max.shape)
diff --git a/Tests/plot_snapshot.py b/Analysis/plot_snapshot.py
similarity index 87%
rename from Tests/plot_snapshot.py
rename to Analysis/plot_snapshot.py
index e2781e4baa5d7476952f592d2ab69efb499044de..ec5ded3502f8e896669f1c7933a3911e6c039ac5 100755
--- a/Tests/plot_snapshot.py
+++ b/Analysis/plot_snapshot.py
@@ -34,6 +34,9 @@ parser.add_argument("-w", "--window", help = "choose a limited plot window",
action = "store_true")
parser.add_argument("-l", "--light", help = "lighter plot",
action = "store_true")
+parser.add_argument("-e", "--end_filename", help = "add character string after "
+ "extremum, outermost contour and max_speed_contour",
+ default = "")
args = parser.parse_args()
with netCDF4.Dataset("h.nc") as f:
@@ -68,11 +71,11 @@ if args.grid:
ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
marker = "+", color = "gray", linestyle = "None")
-reader_extr = shapefile.Reader("extremum_1")
-reader_outer = shapefile.Reader("outermost_contour_1")
+reader_extr = shapefile.Reader("extremum" + args.end_filename)
+reader_outer = shapefile.Reader("outermost_contour" + args.end_filename)
try:
- reader_m_s = shapefile.Reader("max_speed_contour_1")
+ reader_m_s = shapefile.Reader("max_speed_contour" + args.end_filename)
except shapefile.ShapefileException:
print("Shapefile max_speed_contour_1 not found. "
"Max-speed contours will not be plotted.")
@@ -82,8 +85,7 @@ else:
m_s_iterShapes = reader_m_s.iterShapes()
for shape_rec_extr, shape_outer, shape_m_s \
- in zip(reader_extr.iterShapeRecords(), reader_outer.iterShapes(),
- m_s_iterShapes):
+ in zip(reader_extr, reader_outer.iterShapes(), m_s_iterShapes):
points = shape_rec_extr.shape.points[0]
if args.window:
@@ -92,7 +94,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
if not args.window or (points[0] <= urcrnrlon
and llcrnrlat <= points[1] <= urcrnrlat):
- if shape_rec_extr.record[4] == 0:
+ if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "red"
else:
@@ -101,7 +103,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
lines = ax.plot(points[0], points[1], markersize = 10, color = color,
fillstyle = "none", transform = src_crs)
- if shape_rec_extr.record[5] != 0:
+ if shape_rec_extr.record.valid != 0:
if args.light:
lines[0].set_marker("+")
else:
@@ -111,7 +113,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
lines[0].set_marker("s")
if not args.light:
- ax.annotate(str(shape_rec_extr.record[2]),
+ ax.annotate(str(shape_rec_extr.record.eddy_index),
projection.transform_point(points[0], points[1],
src_crs),
xytext = (3, 3), textcoords = "offset points")
@@ -122,7 +124,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
transform = src_crs)
if not args.light:
- if shape_rec_extr.record[5] == 0:
+ if shape_rec_extr.record.valid == 0:
# Invalid outermost contour
lines[0].set_marker("s")
lines[0].set_fillstyle("none")
@@ -134,7 +136,7 @@ for shape_rec_extr, shape_outer, shape_m_s \
if shape_m_s.shapeType != shapefile.NULL:
points = np.array(shape_m_s.points)
- if shape_rec_extr.record[4] == 0:
+ if shape_rec_extr.record.cyclone == 0:
# Anti-cyclone
color = "magenta"
else:
diff --git a/Tests/read_overlap.py b/Analysis/read_overlap.py
similarity index 100%
rename from Tests/read_overlap.py
rename to Analysis/read_overlap.py
diff --git a/Tests/stat.py b/Analysis/stat.py
similarity index 72%
rename from Tests/stat.py
rename to Analysis/stat.py
index fb2de32809b3fddb4b37518a21edbbfe6dc04f71..e32e4edb3a90b3a9dc1ae5c42a4148aa5eb8291f 100755
--- a/Tests/stat.py
+++ b/Analysis/stat.py
@@ -1,8 +1,15 @@
#!/usr/bin/env python3
+import argparse
import shapefile
import numpy as np
+parser = argparse.ArgumentParser()
+parser.add_argument("-e", "--end_filename", help = "add character string after "
+ "extremum, outermost contour and max_speed_contour",
+ default = "")
+args = parser.parse_args()
+
n_valid_cycl = 0
n_valid = 0
n_minima = 0
@@ -11,9 +18,11 @@ n_outer = 0
n_points_valid = 0
n_points = 0
n_missing_speed = np.zeros(3, dtype = int)
+n_valid_speed = 0
-with shapefile.Reader("extremum_1") as extremum, \
- shapefile.Reader("outermost_contour_1") as outermost_contour:
+with shapefile.Reader("extremum" + args.end_filename) as extremum, \
+ shapefile.Reader("outermost_contour" + args.end_filename) \
+ as outermost_contour:
n_extr = len(extremum)
for rec_extr, shape_rec_outer in zip(extremum.iterRecords(),
@@ -33,7 +42,7 @@ with shapefile.Reader("extremum_1") as extremum, \
n_radius4[0] += 1
elif shape_rec_outer.record.radius4 == 1:
n_radius4[1] += 1
- else:
+ elif shape_rec_outer.record.radius4 >= 2:
n_radius4[2] += 1
if rec_extr.speed == 1e4:
@@ -41,8 +50,10 @@ with shapefile.Reader("extremum_1") as extremum, \
n_missing_speed[0] += 1
elif shape_rec_outer.record.radius4 == 1:
n_missing_speed[1] += 1
- else:
+ elif shape_rec_outer.record.radius4 >= 2:
n_missing_speed[2] += 1
+ else:
+ n_valid_speed += 1
print("number of extrema of SSH:", n_extr)
print("number of minima of SSH:", n_minima)
@@ -71,16 +82,20 @@ print("average number of points of outermost contour per extremum (valid eddy "
"or not):", n_points / n_extr)
print()
+print("number of valid speed values:", n_valid_speed)
print("number of missing speed values with radius4 == 0:", n_missing_speed[0])
print("number of missing speed values with radius4 == 1:", n_missing_speed[1])
-print("number of missing speed values with radius4 >= 2:", n_missing_speed[2])
+print("number of missing speed values with radius4 >= 2 (non null shape in "
+ "max_speed_contour and missing speed in extremum.dbf):",
+ n_missing_speed[2])
print()
n_max_speed = 0
n_points_valid = 0
n_points = 0
-with shapefile.Reader("max_speed_contour_1") as max_speed_contour:
+with shapefile.Reader("max_speed_contour" + args.end_filename) \
+ as max_speed_contour:
for shape_rec in max_speed_contour:
l = len(shape_rec.shape.points)
n_points += l
@@ -89,9 +104,12 @@ with shapefile.Reader("max_speed_contour_1") as max_speed_contour:
n_points_valid += l
print("number of non-null maximum-speed contours (distinct from the "
- "outermost contour):",
+ "outermost contour, implies radius4 >= 2):",
n_max_speed)
-print("average number of points of non-null maximum-speed contours:",
- n_points_valid / n_max_speed)
+
+if n_max_speed != 0:
+ print("average number of points of non-null maximum-speed contours:",
+ n_points_valid / n_max_speed)
+
print("average number of points of maximum-speed contours, null or not, per "
"extremum:", n_points / n_extr)
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index ad064dbf839ca0f2bedb3bd25bcf706b9a19fe93..c3cc3e3960d5aca71f19ba9cdbe136529e2cb6c0 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -118,7 +118,10 @@ set\_max\_speed est NaN. Si radius4 $\ge$ 2 et speed dans extremum.dbf
vaut missing\_speed alors speed\_outerm dans set\_max\_speed est NaN,
e\%speed\_cont est non vide dans max\_speed\_contour.shp et
mean\_speed sur e\%speed\_cont a donné un NaN. En d'autres termes, les
-deux appels à mean\_speed ont tous les deux produit un NaN.
+deux appels à mean\_speed ont tous les deux produit un
+NaN. (e\%speed\_cont est non vide dans max\_speed\_contour.shp et
+speed dans extremum.dbf vaut missing\_speed) si et seulement si
+(radius4 $\ge$ 2 et speed dans extremum.dbf vaut missing\_speed).
Cas où les composantes out\_cont ou speed\_cont d'un tourbillon sont
vides, c'est-à-dire égales à null\_ssh\_contour(). On peut distinguer