Define functions fig_hist, fig_distr_funct and plot_all in module

distribution_function. This allows to plot comparisons of various
combinations of runs.

Update plot_snapshot.py to work with version 2 of pyshp.

Analysis/distribution_function.py

@@ -3,16 +3,17 @@
 import shapefile
 import numpy as np
 import matplotlib.pyplot as plt
+import itertools
 
-def plot_distr_funct(X, label = None, ax = None):
-    """Sort and plot distribution function."""
+def plot_distr_funct(x, label = None, ax = None):
+    """Sort and plot distribution function. x should be a Numpy array."""
     
-    X.sort()
-    nx = np.size(X)
+    x.sort()
+    nx = np.size(x)
     if ax is None: ax = plt.gca()
-    ax.plot(X, (1 + np.arange(nx)) / nx, label = label)
-    print("minimum value:", X[0])
-    print("maximum value:", X[-1])
+    ax.plot(x, (1 + np.arange(nx)) / nx, label = label)
+    print("minimum value:", x[0])
+    print("maximum value:", x[-1])
     print("number of values:", nx, "\n")
 
 def read(end_filename = ""):
@@ -51,90 +52,84 @@ def read(end_filename = ""):
             "rad_speed": np.array(rad_speed), "amp_outer": np.array(amp_outer),
             "amp_speed": np.array(amp_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()
-    d = read(args.end_filename)
-
-    xlabel = "speed, in m s-1"
+def fig_hist(xlabel, x, label):
+    """Complete figure with histogram. x and label can also be sequences."""
     
     plt.figure()
-    plt.hist(d["speed"], bins = "auto")
+    plt.hist(x, bins = "auto", histtype = "step", label = label)
     plt.xlabel(xlabel)
     plt.ylabel("number of occurrences")
+    if label is not None: plt.legend()
 
-    plt.figure()
-    print("speed, when positive:")
-    plot_distr_funct(d["speed"][d["speed"] >= 0], label = "positive")
-    print("- speed, when 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(d["rad_outer"], bins = "auto")
-    plt.xlabel(xlabel)
-    plt.ylabel("number of occurrences")
+def fig_distr_funct(xlabel, x, label = None):
+    """Complete figure with distribution function. x must be a sequence of
+    Numpy arrays and label None or a sequence of strings.
 
-    plt.figure()
-    print(xlabel + ":")
-    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(d["rad_speed"], bins = "auto")
-    plt.xlabel(xlabel)
-    plt.ylabel("number of occurrences")
-
     plt.figure()
     print(xlabel + ":")
-    plot_distr_funct(d["rad_speed"])
+
+    label_None = label is None
+    if label_None: label = itertools.repeat(None)
+        
+    for xi, li in zip(x, label):
+        plot_distr_funct(xi, label = li)
+        
     plt.xlabel(xlabel)
     plt.ylabel("distribution function")
-
-    xlabel = "amplitude of outermost contour, in m"
+    if not label_None: plt.legend()
     
-    plt.figure()
-    plt.hist(d["amp_outer"], bins = "auto")
-    plt.xlabel(xlabel)
-    plt.ylabel("number of occurrences")
+def plot_all(dict_list, label = None):
+    """dict_list: list of dictionaries. label: list of labels, one label
+    for each dictionary.
 
-    plt.figure()
-    print("amplitude of outermost contour, anticyclones:")
-    plot_distr_funct(d["amp_outer"][d["amp_outer"] >= 0],
-                     label = "anticyclones")
-    print("amplitude of outermost contour, cyclones:")
-    plot_distr_funct(- d["amp_outer"][d["amp_outer"] < 0], label = "cyclones")
-    plt.xlabel(xlabel)
-    plt.ylabel("distribution function")
-    plt.legend()
+    """
     
-    xlabel = "amplitude of max-speed contour, in m"
+    fig_hist("speed, in m s-1", [d["speed"] for d in dict_list], label)
+    fig_distr_funct("speed, when positive, in m s-1",
+                    [d["speed"][d["speed"] >= 0] for d in dict_list], label)
+    fig_distr_funct("speed, when negative, in m s-1",
+                    [- d["speed"][d["speed"] < 0] for d in dict_list], label)
     
-    plt.figure()
-    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(d["amp_speed"][d["amp_speed"] >= 0],
-                     label = "anticyclones")
-    print("amplitude of max-speed contour, cyclones:")
-    plot_distr_funct(- d["amp_speed"][d["amp_speed"] < 0], label = "cyclones")
-    plt.xlabel(xlabel)
-    plt.ylabel("distribution function")
-    plt.legend()
+    fig_hist("equivalent radius of outermost contour, in km",
+             [d["rad_outer"] for d in dict_list], label)
+    fig_distr_funct("equivalent radius of outermost contour, in km",
+                    [d["rad_outer"] for d in dict_list], label)
+
+    fig_hist("equivalent radius of max-speed contour, in km",
+             [d["rad_speed"] for d in dict_list], label)
+    fig_distr_funct("equivalent radius of max-speed contour, in km",
+                    [d["rad_speed"] for d in dict_list], label)
+
+    fig_hist("amplitude of outermost contour, in m",
+             [d["amp_outer"] for d in dict_list], label)
+    fig_distr_funct("amplitude of outermost contour, anticyclones, in m",
+                    [d["amp_outer"][d["amp_outer"] >= 0] for d in dict_list],
+                    label)
+    fig_distr_funct("amplitude of outermost contour, cyclones, in m",
+                    [- d["amp_outer"][d["amp_outer"] < 0] for d in dict_list],
+                    label)
+
+    fig_hist("amplitude of max-speed contour, in m",
+             [d["amp_speed"] for d in dict_list], label)
+    fig_distr_funct("amplitude of max-speed contour, anticyclones, in m",
+                    [d["amp_speed"][d["amp_speed"] >= 0] for d in dict_list],
+                    label)
+    fig_distr_funct("amplitude of max-speed contour, cyclones, in m",
+                    [- d["amp_speed"][d["amp_speed"] < 0] for d in dict_list],
+                    label)
+    
+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()
+    d = read(args.end_filename)
+    plot_all([d])
     
     plt.show()

Analysis/plot_snapshot.py

@@ -48,8 +48,7 @@ try:
 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)
+    m_s_iterShapes = itertools.repeat(None)
 else:
     m_s_iterShapes = reader_m_s.iterShapes()
 
@@ -129,12 +128,12 @@ for shape_rec_extr, shape_outer, shape_m_s \
                     # Invalid outermost contour
                     lines[0].set_marker("s")
                     lines[0].set_fillstyle("none")
-                elif shape_m_s.shapeType == shapefile.NULL:
+                elif shape_m_s == None or shape_m_s.shapeType == shapefile.NULL:
                     lines[0].set_marker("x")
                 else:
                     lines[0].set_marker("o")
 
-            if shape_m_s.shapeType != shapefile.NULL:
+            if shape_m_s != None and shape_m_s.shapeType != shapefile.NULL:
                 points = np.array(shape_m_s.points)
                 
                 if shape_rec_extr.record.cyclone == 0: