Add default value of argument end_filename of function def read in

script distribution_function.py. Return dictionary instead of tuple:
easier to use.

In script plot_snapshot.py, move definition of readers to the
beginning so there is no risk of entering the window for nothing if we
forgot the end_filename option. Also, enter the four values defining
the window in a single input: easier to copy and paste.

Analysis/distribution_function.py

@@ -15,9 +15,10 @@ def plot_distr_funct(X, label = None, ax = None):
     print("maximum value:", X[-1])
     print("number of values:", nx, "\n")
 
-def read(end_filename):
-    """Return speed, radius of outermost contour and radius of maximum speed
-    contour from three dbf files, as Numpy arrays.
+def read(end_filename = ""):
+    """Read the three dbf files and return speed, radius and amplitude of
+    outermost contour, radius and amplitude of maximum speed contour,
+    as Numpy arrays.
 
     Select valid speed values and valid outermost contour (with
     sufficient area).
@@ -46,8 +47,9 @@ def read(end_filename):
                 rad_speed.append(rec_max.r_eq_area)
                 amp_speed.append(rec_extr.ssh - rec_max.ssh)
 
-    return np.array(speed), np.array(rad_outer), np.array(rad_speed), \
-        np.array(amp_outer), np.array(amp_speed)
+    return {"speed": np.array(speed), "rad_outer": np.array(rad_outer),
+            "rad_speed": np.array(rad_speed), "amp_outer": np.array(amp_outer),
+            "amp_speed": np.array(amp_speed)}
 
 if __name__ == "__main__":
     import argparse
@@ -58,62 +60,62 @@ if __name__ == "__main__":
                         "extremum, outermost_contour and max_speed_contour",
                         default = "")
     args = parser.parse_args()
-
-    speed, rad_outer, rad_speed, amp_outer, amp_speed = read(args.end_filename)
+    d = read(args.end_filename)
 
     xlabel = "speed, in m s-1"
     
     plt.figure()
-    plt.hist(speed, bins = "auto")
+    plt.hist(d["speed"], bins = "auto")
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
 
     plt.figure()
     print("speed, when positive:")
-    plot_distr_funct(speed[speed >= 0], label = "positive")
+    plot_distr_funct(d["speed"][d["speed"] >= 0], label = "positive")
     print("- speed, when negative:")
-    plot_distr_funct(- speed[speed < 0], label = "negative")
+    plot_distr_funct(- d["speed"][d["speed"] < 0], label = "negative")
     plt.xlabel(xlabel)
     plt.legend()
     
     xlabel = "equivalent radius of outermost contour, in km"
 
     plt.figure()
-    plt.hist(rad_outer, bins = "auto")
+    plt.hist(d["rad_outer"], bins = "auto")
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
 
     plt.figure()
     print(xlabel + ":")
-    plot_distr_funct(rad_outer)
+    plot_distr_funct(d["rad_outer"])
     plt.xlabel(xlabel)
     plt.ylabel("distribution function")
     
     xlabel = "equivalent radius of max-speed contour, in km"
 
     plt.figure()
-    plt.hist(rad_speed, bins = "auto")
+    plt.hist(d["rad_speed"], bins = "auto")
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
 
     plt.figure()
     print(xlabel + ":")
-    plot_distr_funct(rad_speed)
+    plot_distr_funct(d["rad_speed"])
     plt.xlabel(xlabel)
     plt.ylabel("distribution function")
 
     xlabel = "amplitude of outermost contour, in m"
     
     plt.figure()
-    plt.hist(amp_outer, bins = "auto")
+    plt.hist(d["amp_outer"], bins = "auto")
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
 
     plt.figure()
     print("amplitude of outermost contour, anticyclones:")
-    plot_distr_funct(amp_outer[amp_outer >= 0], label = "anticyclones")
+    plot_distr_funct(d["amp_outer"][d["amp_outer"] >= 0],
+                     label = "anticyclones")
     print("amplitude of outermost contour, cyclones:")
-    plot_distr_funct(- amp_outer[amp_outer < 0], label = "cyclones")
+    plot_distr_funct(- d["amp_outer"][d["amp_outer"] < 0], label = "cyclones")
     plt.xlabel(xlabel)
     plt.ylabel("distribution function")
     plt.legend()
@@ -121,15 +123,16 @@ if __name__ == "__main__":
     xlabel = "amplitude of max-speed contour, in m"
     
     plt.figure()
-    plt.hist(amp_speed, bins = "auto")
+    plt.hist(d["amp_speed"], bins = "auto")
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
 
     plt.figure()
     print("amplitude of max-speed contour, anticyclones:")
-    plot_distr_funct(amp_speed[amp_speed >= 0], label = "anticyclones")
+    plot_distr_funct(d["amp_speed"][d["amp_speed"] >= 0],
+                     label = "anticyclones")
     print("amplitude of max-speed contour, cyclones:")
-    plot_distr_funct(- amp_speed[amp_speed < 0], label = "cyclones")
+    plot_distr_funct(- d["amp_speed"][d["amp_speed"] < 0], label = "cyclones")
     plt.xlabel(xlabel)
     plt.ylabel("distribution function")
     plt.legend()

Analysis/plot_snapshot.py

@@ -4,10 +4,11 @@
 
 Red for anti-cyclones, blue for cyclones. Empty circles for extrema
 with a valid outermost contour, empty squares for other
-extrema. Squares for a well-defined but invalid outermost
-contour. Crosses for a valid outermost contour but with no distinct
-max-speed contour. Filled circles for a valid outermost contour with a
-distinct max-speed contour.
+extrema. Squares on outermost contour for a well-defined but invalid
+outermost contour. Crosses on outermost contour for a valid outermost
+contour but with no distinct max-speed contour. Filled circles on
+outermost contour and max-speed contour for a valid outermost contour
+with a distinct max-speed contour.
 
 This scripts takes about 30 s of CPU for a 90° by 90° window, about 4
 mn for a 180° by 180° window.
@@ -39,15 +40,27 @@ parser.add_argument("-e", "--end_filename", help = "add character string after "
                     default = "")
 args = parser.parse_args()
 
+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" + args.end_filename)
+except shapefile.ShapefileException:
+    print("Shapefile max_speed_contour_1 not found. "
+          "Max-speed contours will not be plotted.")
+    w = shapefile.Writer(shapeType = shapefile.NULL)
+    m_s_iterShapes = itertools.repeat(w)
+else:
+    m_s_iterShapes = reader_m_s.iterShapes()
+
 with netCDF4.Dataset("h.nc") as f:
     longitude = f.variables["lon"][:]
     latitude = f.variables["lat"][:]
     
 if args.window:
-    llcrnrlon = float(input("llcrnrlon = ? "))
-    llcrnrlat = float(input("llcrnrlat = ? "))
-    urcrnrlon = float(input("urcrnrlon = ? "))
-    urcrnrlat = float(input("urcrnrlat = ? "))
+    l = input("llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = ? ")\
+                                                            .split(sep = ",")
+    llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = [float(v) for v in l]
 
     if urcrnrlon - llcrnrlon > 360:
         sys.exit("bad values of urcrnrlon and llcrnrlon")
@@ -64,6 +77,7 @@ else:
 plt.figure()
 src_crs = ccrs.PlateCarree()
 projection = ccrs.PlateCarree()
+##projection = ccrs.NorthPolarStereo()
 ax = plt.axes(projection = projection)
 
 if args.grid:
@@ -71,19 +85,6 @@ 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" + args.end_filename)
-reader_outer = shapefile.Reader("outermost_contour" + args.end_filename)
-
-try:
-    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.")
-    w = shapefile.Writer(shapeType = shapefile.NULL)
-    m_s_iterShapes = itertools.repeat(w)
-else:
-    m_s_iterShapes = reader_m_s.iterShapes()
-
 for shape_rec_extr, shape_outer, shape_m_s \
     in zip(reader_extr, reader_outer.iterShapes(), m_s_iterShapes):
     points = shape_rec_extr.shape.points[0]
@@ -103,7 +104,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.valid != 0:
+        if shape_rec_extr.record.valid == 1:
             if args.light:
                 lines[0].set_marker("+")
             else:

Documentation_texfol/documentation.tex

@@ -20,6 +20,8 @@
 
 \usepackage{algorithm}
 
+\newcommand{\ud}{\mathrm{d}}
+
 \renewcommand{\algorithmicfor}{\textbf{pour}}
 \renewcommand{\algorithmicend}{\textbf{fin}}
 \renewcommand{\algorithmicdo}{\textbf{faire}}
@@ -1277,6 +1279,55 @@ weight à u et v, c'est-à-dire qu'il faudrait conserver u, v à
 plusieurs dates en dehors de get\_snpashot. Cela remettrait en
 question tout l'algorithme.
 
+\subsection{spherical\_polyline\_area}
+
+La procédure de calcul des contours donne une suite de sommets. Les
+courbes reliant ces sommets ne sont en fait pas imposés
+intrinsèquement par cette procédure. Par exemple, il n'est ni plus ni
+moins juste de supposer que ces courbes sont des portions de grands
+cercles ou de supposer qu'elles donnent des segments en projection
+plate carrée. Néanmoins, dans le reste du programme, pour décider si
+un contour entoure un extremum, ou pour calculer l'intersection de
+deux tourbillons, nous supposons que les contours sont des polygones
+en $(\lambda, \phi)$ (c'est-à-dire que les courbes reliant les sommets
+donnent des segments en projection plate carrée).
+
+Un domaine $S$ quelconque à la surface de la Terre a pour aire :
+\begin{equation*}
+  a^2  \iint_S \cos \phi\ \ud \phi\ \ud \lambda
+  = a^2  \iint_S \ud \lambda\ \ud(\sin \phi)
+\end{equation*}
+Nous avons une procédure qui calcule la surface $\iint \ud x\ \ud y$
+d'un polygone en $(x, y)$. Nous calculons la surface du polygone en
+$(\lambda, \sin \phi)$ passant par les sommets du
+contour. C'est-à-dire que nous supposons que les courbes reliant les
+somments donnent des segments en projection conforme cylindrique.
+
+Rémi fait un autre calcul : la surface du polygone en
+$(\lambda \cos \phi, \phi)$ passant par les sommets du contour. Le
+changement de variable :
+\begin{align*}
+  & x = \lambda \cos \phi \\
+  & y = \phi
+\end{align*}
+a pour jacobien $\cos \phi$ donc on a bien :
+\begin{equation*}
+  \iint \cos \phi\ \ud \phi\ \ud \lambda = \iint \ud x\ \ud y
+\end{equation*}
+Les deux surfaces obtenues, celle du polygone en
+$(\lambda, \sin \phi)$ et celle du polygone en
+$(\lambda \cos \phi, \phi)$ sont très proches si les sommets sont très
+proches. C'est le cas puisque les sommets sont séparés au maximum de
+\np{0.25}° en longitude et latitude. On peut aussi le voir en
+considérant la formule de calcul de surface du polygone :
+\begin{align*}
+  S & = \frac{1}{2} \sum_i x_i (y_{i + 1} - y_{i - 1}) \\
+    & =
+      \frac{1}{2} \sum_i \lambda_i \cos \phi_i (\phi_{i + 1} - \phi_{i - 1}) \\
+    & \approx \frac{1}{2} \sum_i \lambda_i (\sin \phi_{i + 1} - \sin \phi_{i - 1})
+\end{align*}
+si les valeurs successives des $\phi_i$ sont proches.
+
 \subsection{successive\_overlap}
 
 j, k, scalaires entiers, données. flow : donnée-résultat. Cherche les

spherical_polygon_area.f

@@ -6,8 +6,8 @@ contains
 
   pure real function spherical_polygon_area(p)
 
-    ! Assuming p is a polygon in longitude, latitude, in rad, compute
-    ! the area of a polygon in longitude, sin(latitude) with the same
+    ! Assuming p is a polygon in longitude, latitude, compute the area
+    ! of a polygon in longitude, sin(latitude) with the same
     ! vertices. Result in m2.
 
     ! Libraries:
@@ -15,7 +15,7 @@ contains
     
     use spherical_polyline_area_m, only: spherical_polyline_area
 
-    type(polygon), intent(in):: p
+    type(polygon), intent(in):: p ! in rad
 
     ! Local:
 

spherical_polyline_area.f

@@ -6,15 +6,15 @@ contains
 
   pure real function spherical_polyline_area(p)
 
-    ! Assuming p is a polyline in longitude, latitude, in rad, compute
-    ! the area (positive) of a polyline in longitude, sin(latitude)
-    ! with the same vertices. Result in m2.
+    ! Assuming p is a polyline in longitude, latitude, compute the
+    ! area (positive) of a polyline in longitude, sin(latitude) with
+    ! the same vertices. Result in m2.
 
     ! Libraries:
     use contour_531, only: polyline
     use geometry, only: polygon_area_2d
 
-    type(polyline), intent(in):: p ! should be closed
+    type(polyline), intent(in):: p ! Should be closed. In rad.
 
     ! Local: