Skip to content
Snippets Groups Projects
calc.py 4.04 KiB
Newer Older
kgao's avatar
kgao committed
import math

import numpy as np

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

NaN = float("nan")


def is_nan(a):
    return a != a
kgao's avatar
kgao committed


def knots_to_mps(knots):
    return knots * 0.51444


David Hoese's avatar
David Hoese committed
def dewpoint(temp_c, relhum):
    """Convert air temperature and relative humidity to dewpoint.
kgao's avatar
kgao committed

David Hoese's avatar
David Hoese committed
    Algorithm from Tom Whittaker.

    :param temp_c: temperature in celsius
kgao's avatar
kgao committed
    :param relhum: relative humidity as a percentage
    """
David Hoese's avatar
David Hoese committed
    if temp_c is None or relhum is None:
kgao's avatar
kgao committed
        return NaN

    gasconst = 461.5
    latheat = 2500800.0

David Hoese's avatar
David Hoese committed
    dp = 1.0 / (1.0 / (273.15 + temp_c) - gasconst * np.log((0.0 + relhum) / 100) / (latheat - temp_c * 2397.5))
kgao's avatar
kgao committed

    if pd is not None and isinstance(dp, pd.Series):
David Hoese's avatar
David Hoese committed
        return pd.concat([dp - 273.15, temp_c], axis=1).min(axis=1)
    return np.min(dp - 273.15, temp_c)
kgao's avatar
kgao committed


David Hoese's avatar
David Hoese committed
def relhum(air_temp_k, dewpoint_temp_k):
    """Calculate relative humidity from air temperature and dewpoint temperature.
kgao's avatar
kgao committed

David Hoese's avatar
David Hoese committed
    :param air_temp_k: air temperature in Kelvin
    :param dewpoint_temp_k: dewpoint temp in Kelvin
kgao's avatar
kgao committed
    """
David Hoese's avatar
David Hoese committed
    if air_temp_k is None or dewpoint_temp_k is None:
kgao's avatar
kgao committed
        return NaN

    gas_constant = 461.5
    latheat = 2500800.0

    # Only one section of the equation
David Hoese's avatar
David Hoese committed
    latpart = latheat - (air_temp_k - 273.15) * 2397.5
    return 100 * math.e ** ((latpart / air_temp_k - latpart / dewpoint_temp_k) / gas_constant)
kgao's avatar
kgao committed


David Hoese's avatar
David Hoese committed
def potentialtemp(air_temp_k, pressure_mb):
    """Algorithm from David Hoese to calculate potential temperature.
kgao's avatar
kgao committed

David Hoese's avatar
David Hoese committed
    :param air_temp_k: air temperature in Kelvin
    :param pressure_mb: air pressure in millibars
kgao's avatar
kgao committed
    """
David Hoese's avatar
David Hoese committed
    if air_temp_k is None or pressure_mb is None:
kgao's avatar
kgao committed
        return NaN

David Hoese's avatar
David Hoese committed
    return air_temp_k * (pressure_mb.max() / pressure_mb) ** 0.286
kgao's avatar
kgao committed


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 = 0.190284
    c1 = 0.0065 * pow(1013.25, n) / 288.0
    c2 = alt / pow((p - 0.3), n)
    ff = pow(1.0 + c1 * c2, 1.0 / n)
    return (p - 0.3) * ff * 29.92 / 1013.25
kgao's avatar
kgao committed


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'
    """
David Hoese's avatar
David Hoese committed
    if not (val >= 0 and val < 360):  # noqa: PLR2004
        msg = f"'{val}' out of range"
        raise ValueError(msg)
    cardinal_dirs = ("NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW")
kgao's avatar
kgao committed

David Hoese's avatar
David Hoese committed
    if (val >= 348.75 and val <= 360) or val >= 0 and val < 11.25:  # noqa: PLR2004
kgao's avatar
kgao committed

    # 1/2 degree increment between the directions
David Hoese's avatar
David Hoese committed
    for card_dir in cardinal_dirs:
kgao's avatar
kgao committed
        if val >= i and val < (i + 22.5):
David Hoese's avatar
David Hoese committed
            return card_dir
kgao's avatar
kgao committed
        i += 22.5
kgao's avatar
kgao committed


def wind_vector_components(windspd, winddir):
David Hoese's avatar
David Hoese committed
    """Decompose polar wind direction and speed into the horizontal and vertical vector components and speed vector.
kgao's avatar
kgao committed

    Inputs can be scalar or arrays.
    """
    dir_rad = np.deg2rad(winddir)
    spd_arr = np.array(windspd)
David Hoese's avatar
David Hoese committed
    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
kgao's avatar
kgao committed


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

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


def mean_wind_vector(windspd, winddir):
David Hoese's avatar
David Hoese committed
    v_e, v_n, v_spd = wind_vector_components(windspd, winddir)
    avg_dir = wind_vector_degrees(np.mean(v_e), np.mean(v_n))
kgao's avatar
kgao committed

David Hoese's avatar
David Hoese committed
    return avg_dir, np.mean(v_spd)