calc.py

import math

import numpy as np

try:
    from pandas import Series
except ImportError:
    # expected to use for isinstance
    Series = np.ndarray

NaN = float('nan')
is_nan = lambda a: a != a


def knots_to_mps(knots):
    return knots * 0.51444


def dewpoint(tempC, relhum):
    """
    Algorithm from Tom Whittaker tempC is the temperature in degrees Celsius,
    relhum is the relative humidity as a percentage.

    :param tempC: temperature in celsius
    :param relhum: relative humidity as a percentage
    """
    if tempC is None or relhum is None:
        return NaN

    gasconst = 461.5
    latheat = 2500800.0

    dp = 1.0 / (1.0 / (273.15 + tempC) - gasconst * np.log((0.0 + relhum) / 100) /
                (latheat - tempC * 2397.5))

    return min(dp - 273.15, tempC)


def relhum(airTempK, dewpointTempK):
    """
    Algorithm derived by David Hoese from the above
    dewpoint(tempC, relhum) function, both parameters are in Kelvin units.

    :param airTempK: air temperature in Kelvin
    :param dewpointTempK: dewpoint temp in Kelvin
    """
    if airTempK == None or dewpointTempK == None:
        return NaN

    gas_constant = 461.5
    latheat = 2500800.0

    # Only one section of the equation
    latpart = (latheat - (airTempK - 273.15) * 2397.5)
    relativehum = 100 * math.e ** ((latpart / airTempK - latpart / dewpointTempK) / gas_constant)

    return relativehum


def potentialtemp(airTempK, pressureMB):
    """
    Algorithm from David Hoese to calculate potential temperature.

    :param airTempK: air temperature in Kelvin
    :param pressureMB: air pressure in millibars
    """
    if airTempK == None or pressureMB == None:
        return NaN

    pT = airTempK * (pressureMB.max() / pressureMB) ** .286

    return pT


def altimeter(p, alt):
    """Compute altimeter from pressure and altitude.

    Converted from code provided by TomW.

    :param p: pressure in hPa.
    :param alt: altitude of the measurement in meters.

    :returns: altimeter in inHg
    """
    n = .190284
    c1 = .0065 * pow(1013.25, n) / 288.
    c2 = alt / pow((p - .3), n)
    ff = pow(1. + c1 * c2, 1. / n)
    return ((p - .3) * ff * 29.92 / 1013.25)


def dir2txt(val):
    """Convert degrees [0, 360) to a textual representation.

    :param val: decimal degrees

    >>> dir2txt(0)
    'N'
    >>> dir2txt(90)
    'E'
    >>> dir2txt(180)
    'S'
    >>> dir2txt(270)
    'W'
    >>> dir2txt(359)
    'N'
    """
    assert val >= 0 and val < 360, "'%s' out of range" % val
    dirs = ("NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW")

    if ((val >= 348.75 and val <= 360) or val >= 0 and val < 11.25): return "N"

    # 1/2 degree increment between the directions
    i = 11.25;
    for dir in dirs:
        if val >= i and val < (i + 22.5):
            return dir
        i += 22.5


def wind_vector_components(windspd, winddir):
    """Decompose scalar or list/array polar wind direction and speed data
    into the horizontal and vertical vector components and speed vector.

    Inputs can be scalar or arrays.
    """
    dir_rad = np.deg2rad(winddir)
    spd_arr = np.array(windspd)
    V_e = spd_arr * np.sin(dir_rad)
    V_n = spd_arr * np.cos(dir_rad)
    U_spd = np.sqrt(pow(V_e, 2) + pow(V_n, 2))
    return V_e, V_n, U_spd


def wind_vector_degrees(vector_east, vector_north):
    """Re-compose horizontal (east/west) and vertical (north/south) vector
    components into wind direction in degrees.

    Inputs can be scalar or arrays.
    """
    rads = np.arctan2(vector_east, vector_north)
    winddir = np.rad2deg(rads)
    if isinstance(winddir, (np.ndarray, Series)):
        winddir[np.less(winddir, 0)] += 360
    elif winddir < 0:
        winddir += 360
    return winddir % 360


def mean_wind_vector(windspd, winddir):
    V_e, V_n, V_spd = wind_vector_components(windspd, winddir)
    avg_dir = wind_vector_degrees(np.mean(V_e), np.mean(V_n))

    return avg_dir, np.mean(V_spd)