Create directory Analysis.

Analysis/distribution_function.py

+#!/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()

Analysis/filter.py

+#!/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)

Tests/plot_snapshot.py → 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:

Tests/stat.py → Analysis/stat.py

 #!/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)

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