distribution_function.py

#!/usr/bin/env python3

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

def plot_distr_funct(X, xlabel, 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):
    """Return speed, radius of outermost contour and radius of maximum speed
    contour from three dbf files, 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 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 np.array(speed), np.array(r_outer), np.array(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 = read(args.end_filename)

    X = speed[speed >= 0]
    xlabel = "speed, when positive, in m s-1"
    
    plt.figure()
    plt.hist(X, bins = "auto")
    plt.xlabel(xlabel)
    plt.ylabel("number of occurrences")

    plt.figure()
    print(xlabel + ":")
    plot_distr_funct(X, "speed, in m s-1")
    ax_distr_speed = plt.gca()

    X = - speed[speed < 0]
    xlabel = "- speed, when negative, in m s-1"
    plt.figure()
    plt.hist(X, bins = "auto")
    plt.xlabel(xlabel)
    plt.ylabel("number of occurrences")

    print(xlabel + ":")
    plot_distr_funct(X, xlabel, ax = ax_distr_speed)
    
    xlabel = "equivalent radius of outermost contour, in km"

    plt.figure()
    plt.hist(r_outer, bins = "auto")
    plt.xlabel(xlabel)
    plt.ylabel("number of occurrences")

    plt.figure()
    print(xlabel + ":")
    plot_distr_funct(r_outer, xlabel)
    plt.xlabel(xlabel)
    plt.ylabel("distribution function")
    
    xlabel = "equivalent radius of max-speed contour, in km"

    plt.figure()
    plt.hist(r_speed, bins = "auto")
    plt.xlabel(xlabel)
    plt.ylabel("number of occurrences")

    plt.figure()
    print(xlabel + ":")
    plot_distr_funct(r_speed, xlabel)
    plt.xlabel(xlabel)
    plt.ylabel("distribution function")
    
    plt.show()