calc.py

import math

import numpy as np

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

NaN = float("nan")


def is_nan(a):
    return a != a


def knots_to_mps(knots):
    return knots * 0.51444


def dewpoint(temp_c, relhum):
    """Convert air temperature and relative humidity to dewpoint.

    Algorithm from Tom Whittaker.

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

    gasconst = 461.5
    latheat = 2500800.0

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

    if pd is not None and isinstance(dp, pd.Series):
        return pd.concat([dp - 273.15, temp_c], axis=1).min(axis=1)
    return np.min(dp - 273.15, temp_c)


def relhum(air_temp_k, dewpoint_temp_k):
    """Calculate relative humidity from air temperature and dewpoint temperature.

    :param air_temp_k: air temperature in Kelvin
    :param dewpoint_temp_k: dewpoint temp in Kelvin
    """
    if air_temp_k is None or dewpoint_temp_k is None:
        return NaN

    gas_constant = 461.5
    latheat = 2500800.0

    # Only one section of the equation
    latpart = latheat - (air_temp_k - 273.15) * 2397.5
    return 100 * math.e ** ((latpart / air_temp_k - latpart / dewpoint_temp_k) / gas_constant)


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

    :param air_temp_k: air temperature in Kelvin
    :param pressure_mb: air pressure in millibars
    """
    if air_temp_k is None or pressure_mb is None:
        return NaN

    return air_temp_k * (pressure_mb.max() / pressure_mb) ** 0.286


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


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

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

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


def wind_vector_components(windspd, winddir):
    """Decompose polar wind direction and speed 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)