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(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)) if pd is not None and isinstance(dp, pd.Series): return pd.concat([dp - 273.15, tempC], axis=1).min(axis=1) return np.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 is None or dewpointTempK is 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 is None or pressureMB is None: return NaN pT = airTempK * (pressureMB.max() / pressureMB) ** 0.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 = 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' """ 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 return None 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)