distribution_function.py

#!/usr/bin/env python3

import shapefile
import numpy as np
import matplotlib.pyplot as plt

def plot_distr_funct(X, label = None, ax = None):
    """Sort and plot distribution function."""
    
    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])
    print("number of values:", nx, "\n")

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).

    """
    
    with shapefile.Reader("extremum" + end_filename) as extremum, \
         shapefile.Reader("outermost_contour" + end_filename) as outerm_cont, \
         shapefile.Reader("max_speed_contour" + end_filename) as max_speed_cont:
        speed = []
        rad_outer = []
        rad_speed = []
        amp_outer = []
        amp_speed = []

        for rec_extr, rec_outer, rec_max in zip(extremum.iterRecords(),
                                                outerm_cont.iterRecords(),
                                                max_speed_cont.iterRecords()):
            if rec_extr.speed != 1e4: speed.append(rec_extr.speed)

            if rec_extr.valid == 1:
                rad_outer.append(rec_outer.r_eq_area)
                amp_outer.append(rec_extr.ssh - rec_outer.ssh)
                
            if rec_max.r_eq_area != - 100:
                rad_speed.append(rec_max.r_eq_area)
                amp_speed.append(rec_extr.ssh - rec_max.ssh)

    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
    
    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"
    
    plt.figure()
    plt.hist(d["speed"], bins = "auto")
    plt.xlabel(xlabel)
    plt.ylabel("number of occurrences")

    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")

    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"])
    plt.xlabel(xlabel)
    plt.ylabel("distribution function")

    xlabel = "amplitude of outermost contour, in m"
    
    plt.figure()
    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(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"
    
    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()
    
    plt.show()