Finding stats for the available Metocean data. (Requires data file fields.py, not included in this gist.)

met.py

#!/usr/bin/env python
# encoding: utf-8
"""
met.py -- Finding stats for regions

2015-7-1: Created by ckunte.
2015-9-11: More dataset options added.

NOTE:
    Two regional sets are available: sarawak, sabah

    Data (ref: MRD) in datasets are in the following 
    specific units:
     -- Hs is in m.
     -- Tp is in s.
     -- d (Average field water depth) is in m.
    
"""
from operator import add
from math import *
from itertools import *

"""
fields.py format: Below is an example of the data 
format of one area with various fields (fields can 
be any in number, and ellipsis (...) is used to 
indicate that there are/can be more such line 
entries):

sar100y = (
{'field': 'AR1', 'hs': 5.7, 'tp': 11.7, 'd': 94.0},
{'field': 'K3', 'hs': 4.6, 'tp': 10.5, 'd': 58.0},
...
{'field': 'ZB9', 'hs': 5.3, 'tp': 11.3, 'd': 48.0}
)

"""
# Import datasets from fields.py
from fields import *

# Constants
g = 9.81    # Acceleration due to gravity (m/s^2)
Y = 0.01005 # Seawater density (MN/m^3)

def main():
    print "Available Metocean datasets: sar100y, sar1m, sab100y, sab1m"
    dataset = input("Enter dataset to use: ")

    # Make filtered lists of Hs, Tp, and d from the chosen dataset
    hs_all = map(lambda x: x['hs'], filter(lambda x: 'hs' in x, dataset))
    tp_all = map(lambda x: x['tp'], filter(lambda x: 'tp' in x, dataset))
    d_all = map(lambda x: x['d'], filter(lambda x: 'd' in x, dataset))

    print "Number of tp entries:", len(tp_all)
    print "Number of hs entries:", len(hs_all)
    print "Number of water depths:", len(d_all)
    # Note: length mismatch between tp, hs, and d will cause TypeError.

    # Wave length (Deepwater)
    if len(tp_all) > 0:
        Ld_all = map(lambda x: g * x**2 / (2 * pi), tp_all)

    # Wave length (Intermediate)
    if len(Ld_all) > 0:
        Li_all = map(lambda x, y: y * sqrt(tanh(2 * pi * x / y)), d_all, Ld_all)

    # Wave length (Shallow)
    if len(tp_all) > 0:
        Ls_all = map(lambda x, y: x * ((g * y)**0.5), tp_all, d_all)

    # WAVE LENGTHS
    # Populate wavelengths in to a list: L_all
    L_all = []

    # Call four separate lists using zip(), compare and 
    # cherry pick one that matches the condition.

    for i, j, k, l in zip(d_all, Ld_all, Li_all, Ls_all):
        # Select appropriate wave length type (deep, 
        # intermediate, or shallow), and value:
        if ((i / j) >= 0.5):
            L = j
        if (0.05 < (i / k) < 0.5):
            L = k
        if ((i / l) <= 0.05):
            L = l
        L_all.append(L)
    #print L_all

    # WAVE NUMBERS
    k_all = []

    for i in L_all:
        k = (2.0 * pi / i)
        k_all.append(k)
    #print k_all

    # WAVE CELERITIES
    c_all = []

    for i, j in zip(L_all, tp_all):
        c = (i / j)
        c_all.append(c)
    #print c_all

    # WAVE ENERGIES
    Energy_all = []

    for i, j in zip(hs_all, L_all):
        Energy = (Y * g * i**2 * j / 8.0)
        Energy_all.append(Energy)
    #print Energy_all

    # WAVE STEEPNESSES
    s_all = []

    for i, j in zip(hs_all, L_all):
        s = i / j
        s_all.append(s)
    #print s_all

    # MAX. WAVE HEIGHT IN A SEASTATE (WEIGEL, 1949)
    h_mhs_all = []
    # AVG. HIGHEST WAVE IN 1,000 WAVES (RAYLEIGH ASSUMPTION)
    h_ahw_all = []
    # MOST PROBABLE HEIGHEST WAVE IN 1,000 WAVES (RAYLEIGH ASSUMPTION)
    h_mphw_all = []
    # STANDARD DEVIATION OF FREE SURFACE
    h_sdfs_all = []
    # MODE WAVE HEIGHT
    h_modh_all = []
    # MEDIAN WAVE HEIGHT
    h_medh_all = []
    # MEAN WAVE HEIGHT
    h_meanh_all = []
    # WAVE HEIGHT ROOT MEAN SQUARE
    h_rms_all = []
    ## POWER WEIGHTED MEAN WAVE HEIGHTS (for Fatigue)
    #h_wmwh1_all = []
    #h_wmwh2_all = []
    #h_wmwh3_all = []
    for i in hs_all:
        h_mhs = i * 1.8
        h_mhs_all.append(h_mhs)
        h_ahw = i * 1.93
        h_ahw_all.append(h_ahw)
        h_mphw = i * 1.86
        h_mphw_all.append(h_mphw)
        h_sdfs = i * 0.25
        h_sdfs_all.append(h_sdfs)
        h_modh = i * 0.499
        h_modh_all.append(h_modh)
        h_medh = i * 0.588
        h_medh_all.append(h_medh)
        h_meanh = i * 0.626
        h_meanh_all.append(h_meanh)
        h_rms = i / sqrt(2.0)
        h_rms_all.append(h_rms)
        #h_wmwh1 = i * 0.776
        #h_wmwh1_all.append(h_wmwh1)
        #h_wmwh2 = i * 0.869
        #h_wmwh2_all.append(h_wmwh2)
        #h_wmwh3 = i * 0.952
        #h_wmwh3_all.append(h_wmwh3)

    params = [hs_all, tp_all, d_all, L_all, k_all, c_all, Energy_all, s_all, 
              h_mhs_all, h_ahw_all, h_mphw_all, h_sdfs_all, 
              h_modh_all, h_medh_all, h_meanh_all, h_rms_all]
    
    print "Legend:"
    print " ( 0 ): Hs -- Significant wave height"
    print " ( 1 ): Tp -- Peak wave period"
    print " ( 2 ): d -- Water depth"
    print " ( 3 ): L -- Wave length"
    print " ( 4 ): k -- Wave number"
    print " ( 5 ): c -- Wave celerity"
    print " ( 6 ): Energy -- Wave energy"
    print " ( 7 ): s -- Wave steepness"
    print " ( 8 ): H_mhs -- Max. height in a seastate (Weigel, 1949)"
    print " ( 9 ): H_ahw -- Avg. highest wave in 1,000 waves (Rayleigh assumption)"
    print "( 10 ): H_mphw -- Most prob. highest wave in 1,000 waves (Rayleigh assumption)"
    print "( 11 ): H_sdfs -- Standard deviation of free surface"
    print "( 12 ): H_modh -- Mode height"
    print "( 13 ): H_medh -- Median height"
    print "( 14 ): h_meanh -- Mean height"
    print "( 15 ): H_rms -- Root mean square"
    print ""
    print "Results for are as follow:"

    for i in params:
        if len(i) > 0:
            mean = reduce(add, i) / len(i)
            vfn = map(lambda x: (x - mean)**2, i)
            variance = reduce(add, vfn) / len(i)
            # Print results (to the second decimal)
            print ""
            print "Maximum (",params.index(i),"):", round(max(i), 4)
            print "Mean    (",params.index(i),"):", round(mean, 4)
            print "Minimum (",params.index(i),"):", round(min(i), 4)
            print "Std.dev.(",params.index(i),"):", round(sqrt(variance), 4)
            print "Variance(",params.index(i),"):", round(variance, 4)
            pass

if __name__ == '__main__':
    main()