From 821647a6957751ea2b91cf293f3e144855863c14 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Wed, 24 Apr 2024 16:24:49 +0100 Subject: [PATCH] Implement array submodule (#6) * Implement array submodule --------- Co-authored-by: Olivier Iffrig --- docs/conf.py | 7 + docs/skip_api_rules.py | 28 + earthkit/meteo/__init__.py | 4 +- earthkit/meteo/constants/__init__.py | 9 + earthkit/meteo/constants/constants.py | 11 +- earthkit/meteo/extreme/__init__.py | 17 + earthkit/meteo/extreme/array/__init__.py | 16 + earthkit/meteo/extreme/array/cpf.py | 99 + earthkit/meteo/extreme/array/efi.py | 148 ++ earthkit/meteo/extreme/array/sot.py | 173 ++ earthkit/meteo/extreme/cpf.py | 102 +- earthkit/meteo/extreme/efi.py | 151 +- earthkit/meteo/extreme/sot.py | 176 +- earthkit/meteo/geo/__init__.py | 1 - earthkit/meteo/geo/geo.py | 173 -- earthkit/meteo/score/__init__.py | 17 + .../meteo/score/array/__init__.py | 11 +- earthkit/meteo/score/array/crps.py | 60 + earthkit/meteo/score/crps.py | 59 +- earthkit/meteo/solar/__init__.py | 237 +-- earthkit/meteo/solar/array/__init__.py | 14 + earthkit/meteo/solar/array/solar.py | 241 +++ earthkit/meteo/solar/solar.py | 34 + earthkit/meteo/stats/__init__.py | 17 + earthkit/meteo/stats/array/__init__.py | 15 + earthkit/meteo/stats/array/numpy_extended.py | 55 + earthkit/meteo/stats/array/quantiles.py | 77 + earthkit/meteo/stats/numpy_extended.py | 54 +- earthkit/meteo/stats/quantiles.py | 76 +- earthkit/meteo/thermo/__init__.py | 18 + earthkit/meteo/thermo/array/__init__.py | 14 + earthkit/meteo/thermo/array/thermo.py | 1791 ++++++++++++++++ earthkit/meteo/thermo/thermo.py | 1794 +---------------- earthkit/meteo/wind/__init__.py | 17 + earthkit/meteo/wind/array/__init__.py | 14 + earthkit/meteo/wind/array/wind.py | 295 +++ earthkit/meteo/wind/wind.py | 305 +-- tests/{ => extreme}/test_extreme.py | 54 +- tests/{ => score}/test_score.py | 3 +- tests/solar/test_solar.py | 27 + tests/{ => stats}/test_stats.py | 3 +- tests/test_geo.py | 171 -- tests/test_version.py | 3 +- tests/thermo/test_thermo.py | 19 + .../test_thermo_array.py} | 177 +- tests/{ => wind}/test_wind.py | 4 +- 46 files changed, 3548 insertions(+), 3243 deletions(-) create mode 100644 docs/skip_api_rules.py create mode 100644 earthkit/meteo/extreme/array/__init__.py create mode 100644 earthkit/meteo/extreme/array/cpf.py create mode 100644 earthkit/meteo/extreme/array/efi.py create mode 100644 earthkit/meteo/extreme/array/sot.py delete mode 100644 earthkit/meteo/geo/__init__.py delete mode 100644 earthkit/meteo/geo/geo.py rename tests/test_solar.py => earthkit/meteo/score/array/__init__.py (75%) create mode 100644 earthkit/meteo/score/array/crps.py create mode 100644 earthkit/meteo/solar/array/__init__.py create mode 100644 earthkit/meteo/solar/array/solar.py create mode 100644 earthkit/meteo/solar/solar.py create mode 100644 earthkit/meteo/stats/array/__init__.py create mode 100644 earthkit/meteo/stats/array/numpy_extended.py create mode 100644 earthkit/meteo/stats/array/quantiles.py create mode 100644 earthkit/meteo/thermo/array/__init__.py create mode 100644 earthkit/meteo/thermo/array/thermo.py create mode 100644 earthkit/meteo/wind/array/__init__.py create mode 100644 earthkit/meteo/wind/array/wind.py rename tests/{ => extreme}/test_extreme.py (94%) rename tests/{ => score}/test_score.py (99%) create mode 100644 tests/solar/test_solar.py rename tests/{ => stats}/test_stats.py (98%) delete mode 100644 tests/test_geo.py create mode 100644 tests/thermo/test_thermo.py rename tests/{test_thermo.py => thermo/test_thermo_array.py} (74%) rename tests/{ => wind}/test_wind.py (99%) diff --git a/docs/conf.py b/docs/conf.py index 05c02a9..cae2e00 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -55,6 +55,7 @@ autoapi_root = "_api" autoapi_member_order = "alphabetical" autoapi_add_toctree_entry = False +autoapi_own_page_level = "function" # napoleon configuration napoleon_google_docstring = False @@ -86,3 +87,9 @@ html_css_files = ["style.css"] html_logo = "_static/earthkit-meteo.png" + + +def setup(app): + from skip_api_rules import _skip_api_items + + app.connect("autoapi-skip-member", _skip_api_items) diff --git a/docs/skip_api_rules.py b/docs/skip_api_rules.py new file mode 100644 index 0000000..e5a9159 --- /dev/null +++ b/docs/skip_api_rules.py @@ -0,0 +1,28 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + + +# define skip rules for autoapi +def _skip_api_items(app, what, name, obj, skip, options): + # print(f"{what=} {name=}") + + if ( + what == "module" + and ".array" not in name + and name not in ["earthkit.meteo.solar", "earthkit.meteo.solar.array"] + ): + skip = True + elif what == "package" and ".array" not in name and len(name.split(".")) > 2: + skip = True + elif what == "function" and ".array" not in name: + skip = True + + # if not skip: + # print(f"{what} {name}") + return skip diff --git a/earthkit/meteo/__init__.py b/earthkit/meteo/__init__.py index c521135..d64217c 100644 --- a/earthkit/meteo/__init__.py +++ b/earthkit/meteo/__init__.py @@ -1,11 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. - +# try: # NOTE: the `version.py` file must not be present in the git repository diff --git a/earthkit/meteo/constants/__init__.py b/earthkit/meteo/constants/__init__.py index 5b1f107..afcccbb 100644 --- a/earthkit/meteo/constants/__init__.py +++ b/earthkit/meteo/constants/__init__.py @@ -1 +1,10 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + from .constants import * # noqa diff --git a/earthkit/meteo/constants/constants.py b/earthkit/meteo/constants/constants.py index 087d849..adc1ba1 100644 --- a/earthkit/meteo/constants/constants.py +++ b/earthkit/meteo/constants/constants.py @@ -1,10 +1,11 @@ -# (C) Copyright 2021- ECMWF. - +# (C) Copyright 2021 ECMWF. +# # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# import numpy as np diff --git a/earthkit/meteo/extreme/__init__.py b/earthkit/meteo/extreme/__init__.py index c42ccca..d3ee0b6 100644 --- a/earthkit/meteo/extreme/__init__.py +++ b/earthkit/meteo/extreme/__init__.py @@ -1,3 +1,20 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Extreme index functions. + +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. +""" + from .cpf import * # noqa from .efi import * # noqa from .sot import * # noqa diff --git a/earthkit/meteo/extreme/array/__init__.py b/earthkit/meteo/extreme/array/__init__.py new file mode 100644 index 0000000..dcf2242 --- /dev/null +++ b/earthkit/meteo/extreme/array/__init__.py @@ -0,0 +1,16 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Extreme index functions operating on numpy arrays. +""" + +from .cpf import * # noqa +from .efi import * # noqa +from .sot import * # noqa diff --git a/earthkit/meteo/extreme/array/cpf.py b/earthkit/meteo/extreme/array/cpf.py new file mode 100644 index 0000000..16596de --- /dev/null +++ b/earthkit/meteo/extreme/array/cpf.py @@ -0,0 +1,99 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + + +def cpf(clim, ens, sort_clim=True, sort_ens=True): + """Compute Crossing Point Forecast (CPF) + + WARNING: this code is experimental, use at your own risk! + + Parameters + ---------- + clim: numpy array (nclim, npoints) + Per-point climatology + ens: numpy array (nens, npoints) + Ensemble forecast + sort_clim: bool + If True, sort the climatology first + sort_ens: bool + If True, sort the ensemble first + + Returns + ------- + numpy array (npoints) + CPF values + """ + nclim, npoints = clim.shape + nens, npoints_ens = ens.shape + assert npoints == npoints_ens + + cpf = np.ones(npoints, dtype=np.float32) + mask = np.zeros(npoints, dtype=np.bool_) + + if sort_clim: + clim = np.sort(clim, axis=0) + if sort_ens: + ens = np.sort(ens, axis=0) + + for icl in range(1, nclim - 1): + # quantile level of climatology + tau_c = icl / (nclim - 1.0) + for iq in range(nens): + # quantile level of forecast + tau_f = (iq + 1.0) / (nens + 1.0) + if tau_f >= tau_c: + # quantile values of forecast and climatology + qv_f = ens[iq, :] + qv_c = clim[icl, :] + + # lowest climate quantile: interpolate between 2 consecutive quantiles + if iq < 2: + # quantile value and quantile level of climatology at previous + qv_c_2 = clim[icl - 1, :] + tau_c_2 = (icl - 1) / (nclim - 1) + + # condition of crossing situtaion: + idx = (qv_f < qv_c) & (qv_c_2 < qv_c) + + # intersection between two lines + tau_i = ( + tau_c * (qv_c_2[idx] - qv_f[idx]) + + tau_c_2 * (qv_f[idx] - qv_c[idx]) + ) / (qv_c_2[idx] - qv_c[idx]) + + # populate matrix, no values below 0 + cpf[idx] = np.maximum(tau_i, 0) + mask[idx] = True + + # check crossing cases + idx = (qv_f < qv_c) & (~mask) + cpf[idx] = tau_f + mask[idx] = True + + # largest climate quantile: interpolate + if iq == nens - 1: + qv_c_2 = clim[nclim - 1, :] + tau_c_2 = 1.0 + + idx = (qv_f > qv_c) & (qv_c_2 > qv_c) & (~mask) + + tau_i = ( + tau_c * (qv_c_2[idx] - qv_f[idx]) + + tau_c_2 * (qv_f[idx] - qv_c[idx]) + ) / (qv_c_2[idx] - qv_c[idx]) + + # populate matrix, no values above 1 + cpf[idx] = np.minimum(tau_i, 1) + + # speed up process + break + + return cpf diff --git a/earthkit/meteo/extreme/array/efi.py b/earthkit/meteo/extreme/array/efi.py new file mode 100644 index 0000000..455222b --- /dev/null +++ b/earthkit/meteo/extreme/array/efi.py @@ -0,0 +1,148 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + +# import numba +# from numba import float64, float32 + + +def efi(clim, ens, eps=-0.1): + """Compute Extreme Forecast Index (EFI) + + Parameters + ---------- + clim: numpy array (nclim, npoints) + Sorted per-point climatology + ens: numpy array (nens, npoints) + Ensemble forecast + eps: (float) + Epsilon factor for zero values + + Returns + ------- + numpy array (npoints) + EFI values + """ + # locate missing values + missing_mask = np.logical_or( + np.sum(np.isnan(clim), axis=0), np.sum(np.isnan(ens), axis=0) + ) + + # Compute fraction of the forecast below climatology + nclim, npoints = clim.shape + nens, npoints_ens = ens.shape + assert npoints == npoints_ens + frac = np.zeros_like(clim) + ################################## + for icl in range(nclim): + frac[icl, :] = np.sum(ens[:, :] <= clim[icl, np.newaxis, :], axis=0) + ################################## + frac /= nens + + # Compute formula coefficients + p = np.linspace(0.0, 1.0, nclim) + dp = 1 / (nclim - 1) + dFdp = np.diff(frac, axis=0) / dp + + acosdiff = np.diff(np.arccos(np.sqrt(p))) + proddiff = np.diff(np.sqrt(p * (1.0 - p))) + + acoef = (1.0 - 2.0 * p[:-1]) * acosdiff + proddiff + + # compute EFI from coefficients + efi = np.zeros(npoints) + ################################## + if eps > 0: + efimax = np.zeros(npoints) + for icl in range(nclim - 1): + mask = clim[icl + 1, :] > eps + dEFI = np.where( + mask, + (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] + + acoef[icl] * dFdp[icl, :] + - proddiff[icl], + 0.0, + ) + defimax = np.where(mask, -acosdiff[icl] - proddiff[icl], 0.0) + efi += dEFI + efimax += defimax + efimax = np.fmax(efimax, eps) + efi /= efimax + else: + for icl in range(nclim - 1): + dEFI = ( + (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] + + acoef[icl] * dFdp[icl, :] + - proddiff[icl] + ) + efi += dEFI + efi *= 2.0 / np.pi + ################################## + + # apply missing values + efi[missing_mask] = np.nan + + return efi + + +# @numba.jit(float64[:](float64[:,:], float64[:,:]), fastmath=False, nopython=True, nogil=True, cache=True) +# @numba.jit(nopython=True) +# def efi_numba(clim, ens): +# """Compute EFI + +# Parameters +# ---------- +# clim: numpy array (nclim, npoints) +# Sorted per-point climatology +# ens: numpy array (nens, npoints) +# Ensemble forecast + +# Returns +# ------- +# numpy array (npoints) +# EFI values +# """ + +# # Compute fraction of the forecast below climatology +# nclim, npoints = clim.shape +# nens, npoints_ens = ens.shape +# assert npoints == npoints_ens +# frac = np.zeros_like(clim) +# ################################## +# for ifo in numba.prange(nens): +# for icl in range(nclim): +# for i in range(npoints): +# if ens[ifo, i] <= clim[icl, i]: +# frac[icl, i] += 1 +# ################################## +# frac /= nens + +# # Compute formula coefficients +# p = np.linspace(0., 1., nclim) +# dp = 1 / (nclim - 1) #np.diff(p) + +# acosdiff = np.diff(np.arccos(np.sqrt(p))) +# proddiff = np.diff(np.sqrt(p * (1. - p))) + +# acoef = (1. - 2. * p[:-1]) * acosdiff + proddiff + +# # TODO: handle epsilon +# efi = np.zeros(npoints) +# ################################## +# for icl in numba.prange(nclim-1): +# for i in range(npoints): +# dFdp = (frac[icl+1, i] - frac[icl, i]) / dp +# # XXX: why proddiff here?! +# dEFI = (2. * frac[icl, i] - 1.) * acosdiff[icl] + acoef[icl] * dFdp - proddiff[icl] +# efi[i] += dEFI +# efi *= 2. / np.pi +# ################################## + +# return efi diff --git a/earthkit/meteo/extreme/array/sot.py b/earthkit/meteo/extreme/array/sot.py new file mode 100644 index 0000000..f556afa --- /dev/null +++ b/earthkit/meteo/extreme/array/sot.py @@ -0,0 +1,173 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + + +def sot_func(qc_tail, qc, qf, eps=-1e-4, lower_bound=-10, upper_bound=10): + """Compute basic Shift of Tails (SOT) + using already computed percentiles + + Parameters + ---------- + qc_tail: numpy array (npoints) + Tail percentile value (99% or 1%) + Model climatology + qc: numpy array (npoints) + Upper or lower percentile (at 90% or 10%) + Model climatology + qf: numpy array (npoints) + Upper or lower percentile (at 90% or 10%) + Ensemble forecast + eps: (float) + Epsilon factor for zero values + missing: (float) + missing points values where denominator is zero + + Returns + ------- + numpy array (npoints) + SOT values + """ + # avoid divided by zero warning + err = np.seterr(divide="ignore", invalid="ignore") + + min_den = np.fmax(eps, 0) + sot = np.where( + np.fabs(qc_tail - qc) > min_den, (qf - qc_tail) / (qc_tail - qc), np.nan + ) + + # revert to original error state + np.seterr(**err) + + mask_missing = np.isnan(sot) + + # upper and lower bounds + mask2 = np.logical_and(np.logical_not(mask_missing), sot < lower_bound) + sot[mask2] = lower_bound + mask3 = np.logical_and(np.logical_not(mask_missing), sot > upper_bound) + sot[mask3] = upper_bound + + return sot + + +def sot(clim, ens, perc, eps=-1e4): + """Compute Shift of Tails (SOT) + from climatology percentiles (sorted) + and ensemble forecast (not sorted) + + Parameters + ---------- + clim: numpy array (nclim, npoints) + Model climatology (percentiles) + ens: numpy array (nens, npoints) + Ensemble forecast + perc: int + Percentile value (typically 10 or 90) + eps: (float) + Epsilon factor for zero values + + Returns + ------- + numpy array (npoints) + SOT values + """ + if not (isinstance(perc, int) or isinstance(perc, np.int64)) or ( + perc < 2 or perc > 98 + ): + raise Exception( + "Percentile value should be and Integer between 2 and 98, is {}".format( + perc + ) + ) + + if clim.shape[0] != 101: + raise Exception( + "Climatology array should contain 101 percentiles, it has {} values".format( + clim.shape + ) + ) + + qc = clim[perc] + # if eps>0, set to zero everything below eps + if eps > 0: + ens = np.where(ens < eps, 0.0, ens) + qc = np.where(qc < eps, 0.0, qc) + + qf = np.percentile(ens, q=perc, axis=0) + if perc > 50: + qc_tail = clim[99] + elif perc < 50: + qc_tail = clim[1] + else: + raise Exception( + "Percentile value to be computed cannot be 50 for sot, has to be in the upper or lower half" + ) + + sot = sot_func(qc_tail, qc, qf, eps=eps) + + return sot + + +def sot_unsorted(clim, ens, perc, eps=-1e4): + """Compute Shift of Tails (SOT) + from climatology percentiles (sorted) + and ensemble forecast (not sorted) + + Parameters + ---------- + clim: numpy array (nclim, npoints) + Model climatology (percentiles) + ens: numpy array (nens, npoints) + Ensemble forecast + perc: int + Percentile value (typically 10 or 90) + eps: (float) + Epsilon factor for zero values + + Returns + ------- + numpy array (npoints) + SOT values + """ + if not (isinstance(perc, int) or isinstance(perc, np.int64)) or ( + perc < 2 or perc > 98 + ): + raise Exception( + "Percentile value should be and Integer between 2 and 98, is {}".format( + perc + ) + ) + + if clim.shape[0] != 101: + raise Exception( + "Climatology array should contain 101 percentiles, it has {} values".format( + clim.shape + ) + ) + + if eps > 0: + ens = np.where(ens < eps, 0.0, ens) + clim = np.where(clim < eps, 0.0, clim) + + qf = np.percentile(ens, q=perc, axis=0) + qc = np.percentile(clim, q=perc, axis=0) + if perc > 50: + perc_tail = 99 + elif perc < 50: + perc_tail = 1 + else: + raise Exception( + "Percentile value to be computed cannot be 50 for sot, has to be in the upper or lower half" + ) + qc_tail = np.percentile(clim, q=perc_tail, axis=0) + + sot = sot_func(qc_tail, qc, qf, eps=eps) + + return sot diff --git a/earthkit/meteo/extreme/cpf.py b/earthkit/meteo/extreme/cpf.py index 8a4e177..5b9e88a 100644 --- a/earthkit/meteo/extreme/cpf.py +++ b/earthkit/meteo/extreme/cpf.py @@ -1,98 +1,14 @@ -# (C) Copyright 2021- ECMWF. - +# (C) Copyright 2021 ECMWF. +# # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. - -import numpy as np - - -def cpf(clim, ens, sort_clim=True, sort_ens=True): - """Compute Crossing Point Forecast (CPF) - - WARNING: this code is experimental, use at your own risk! - - Parameters - ---------- - clim: numpy array (nclim, npoints) - Per-point climatology - ens: numpy array (nens, npoints) - Ensemble forecast - sort_clim: bool - If True, sort the climatology first - sort_ens: bool - If True, sort the ensemble first - - Returns - ------- - numpy array (npoints) - CPF values - """ - nclim, npoints = clim.shape - nens, npoints_ens = ens.shape - assert npoints == npoints_ens - - cpf = np.ones(npoints, dtype=np.float32) - mask = np.zeros(npoints, dtype=np.bool_) - - if sort_clim: - clim = np.sort(clim, axis=0) - if sort_ens: - ens = np.sort(ens, axis=0) - - for icl in range(1, nclim - 1): - # quantile level of climatology - tau_c = icl / (nclim - 1.0) - for iq in range(nens): - # quantile level of forecast - tau_f = (iq + 1.0) / (nens + 1.0) - if tau_f >= tau_c: - # quantile values of forecast and climatology - qv_f = ens[iq, :] - qv_c = clim[icl, :] - - # lowest climate quantile: interpolate between 2 consecutive quantiles - if iq < 2: - # quantile value and quantile level of climatology at previous - qv_c_2 = clim[icl - 1, :] - tau_c_2 = (icl - 1) / (nclim - 1) - - # condition of crossing situtaion: - idx = (qv_f < qv_c) & (qv_c_2 < qv_c) - - # intersection between two lines - tau_i = ( - tau_c * (qv_c_2[idx] - qv_f[idx]) - + tau_c_2 * (qv_f[idx] - qv_c[idx]) - ) / (qv_c_2[idx] - qv_c[idx]) - - # populate matrix, no values below 0 - cpf[idx] = np.maximum(tau_i, 0) - mask[idx] = True - - # check crossing cases - idx = (qv_f < qv_c) & (~mask) - cpf[idx] = tau_f - mask[idx] = True - - # largest climate quantile: interpolate - if iq == nens - 1: - qv_c_2 = clim[nclim - 1, :] - tau_c_2 = 1.0 - - idx = (qv_f > qv_c) & (qv_c_2 > qv_c) & (~mask) - - tau_i = ( - tau_c * (qv_c_2[idx] - qv_f[idx]) - + tau_c_2 * (qv_f[idx] - qv_c[idx]) - ) / (qv_c_2[idx] - qv_c[idx]) - - # populate matrix, no values above 1 - cpf[idx] = np.minimum(tau_i, 1) +from . import array # noqa - # speed up process - break - return cpf +def cpf(*args, **kwargs): + return array.cpf(*args, **kwargs) diff --git a/earthkit/meteo/extreme/efi.py b/earthkit/meteo/extreme/efi.py index 960dd01..cb312b1 100644 --- a/earthkit/meteo/extreme/efi.py +++ b/earthkit/meteo/extreme/efi.py @@ -1,147 +1,14 @@ -# (C) Copyright 2021- ECMWF. - +# (C) Copyright 2021 ECMWF. +# # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. - -import numpy as np - -# import numba -# from numba import float64, float32 - - -def efi(clim, ens, eps=-0.1): - """Compute Extreme Forecast Index (EFI) - - Parameters - ---------- - clim: numpy array (nclim, npoints) - Sorted per-point climatology - ens: numpy array (nens, npoints) - Ensemble forecast - eps: (float) - Epsilon factor for zero values - - Returns - ------- - numpy array (npoints) - EFI values - """ - # locate missing values - missing_mask = np.logical_or( - np.sum(np.isnan(clim), axis=0), np.sum(np.isnan(ens), axis=0) - ) - - # Compute fraction of the forecast below climatology - nclim, npoints = clim.shape - nens, npoints_ens = ens.shape - assert npoints == npoints_ens - frac = np.zeros_like(clim) - ################################## - for icl in range(nclim): - frac[icl, :] = np.sum(ens[:, :] <= clim[icl, np.newaxis, :], axis=0) - ################################## - frac /= nens - - # Compute formula coefficients - p = np.linspace(0.0, 1.0, nclim) - dp = 1 / (nclim - 1) - dFdp = np.diff(frac, axis=0) / dp - - acosdiff = np.diff(np.arccos(np.sqrt(p))) - proddiff = np.diff(np.sqrt(p * (1.0 - p))) - - acoef = (1.0 - 2.0 * p[:-1]) * acosdiff + proddiff - - # compute EFI from coefficients - efi = np.zeros(npoints) - ################################## - if eps > 0: - efimax = np.zeros(npoints) - for icl in range(nclim - 1): - mask = clim[icl + 1, :] > eps - dEFI = np.where( - mask, - (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] - + acoef[icl] * dFdp[icl, :] - - proddiff[icl], - 0.0, - ) - defimax = np.where(mask, -acosdiff[icl] - proddiff[icl], 0.0) - efi += dEFI - efimax += defimax - efimax = np.fmax(efimax, eps) - efi /= efimax - else: - for icl in range(nclim - 1): - dEFI = ( - (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] - + acoef[icl] * dFdp[icl, :] - - proddiff[icl] - ) - efi += dEFI - efi *= 2.0 / np.pi - ################################## - - # apply missing values - efi[missing_mask] = np.nan - - return efi - - -# @numba.jit(float64[:](float64[:,:], float64[:,:]), fastmath=False, nopython=True, nogil=True, cache=True) -# @numba.jit(nopython=True) -# def efi_numba(clim, ens): -# """Compute EFI - -# Parameters -# ---------- -# clim: numpy array (nclim, npoints) -# Sorted per-point climatology -# ens: numpy array (nens, npoints) -# Ensemble forecast - -# Returns -# ------- -# numpy array (npoints) -# EFI values -# """ - -# # Compute fraction of the forecast below climatology -# nclim, npoints = clim.shape -# nens, npoints_ens = ens.shape -# assert npoints == npoints_ens -# frac = np.zeros_like(clim) -# ################################## -# for ifo in numba.prange(nens): -# for icl in range(nclim): -# for i in range(npoints): -# if ens[ifo, i] <= clim[icl, i]: -# frac[icl, i] += 1 -# ################################## -# frac /= nens - -# # Compute formula coefficients -# p = np.linspace(0., 1., nclim) -# dp = 1 / (nclim - 1) #np.diff(p) - -# acosdiff = np.diff(np.arccos(np.sqrt(p))) -# proddiff = np.diff(np.sqrt(p * (1. - p))) - -# acoef = (1. - 2. * p[:-1]) * acosdiff + proddiff +from . import array # noqa -# # TODO: handle epsilon -# efi = np.zeros(npoints) -# ################################## -# for icl in numba.prange(nclim-1): -# for i in range(npoints): -# dFdp = (frac[icl+1, i] - frac[icl, i]) / dp -# # XXX: why proddiff here?! -# dEFI = (2. * frac[icl, i] - 1.) * acosdiff[icl] + acoef[icl] * dFdp - proddiff[icl] -# efi[i] += dEFI -# efi *= 2. / np.pi -# ################################## -# return efi +def efi(*args, **kwargs): + return array.efi(*args, **kwargs) diff --git a/earthkit/meteo/extreme/sot.py b/earthkit/meteo/extreme/sot.py index 3751b50..7546d4a 100644 --- a/earthkit/meteo/extreme/sot.py +++ b/earthkit/meteo/extreme/sot.py @@ -1,172 +1,18 @@ -# (C) Copyright 2021- ECMWF. - +# (C) Copyright 2021 ECMWF. +# # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. - -import numpy as np - - -def sot_func(qc_tail, qc, qf, eps=-1e-4, lower_bound=-10, upper_bound=10): - """Compute basic Shift of Tails (SOT) - using already computed percentiles - - Parameters - ---------- - qc_tail: numpy array (npoints) - Tail percentile value (99% or 1%) - Model climatology - qc: numpy array (npoints) - Upper or lower percentile (at 90% or 10%) - Model climatology - qf: numpy array (npoints) - Upper or lower percentile (at 90% or 10%) - Ensemble forecast - eps: (float) - Epsilon factor for zero values - missing: (float) - missing points values where denominator is zero - - Returns - ------- - numpy array (npoints) - SOT values - """ - # avoid divided by zero warning - err = np.seterr(divide="ignore", invalid="ignore") - - min_den = np.fmax(eps, 0) - sot = np.where( - np.fabs(qc_tail - qc) > min_den, (qf - qc_tail) / (qc_tail - qc), np.nan - ) - - # revert to original error state - np.seterr(**err) - - mask_missing = np.isnan(sot) - - # upper and lower bounds - mask2 = np.logical_and(np.logical_not(mask_missing), sot < lower_bound) - sot[mask2] = lower_bound - mask3 = np.logical_and(np.logical_not(mask_missing), sot > upper_bound) - sot[mask3] = upper_bound - - return sot - - -def sot(clim, ens, perc, eps=-1e4): - """Compute Shift of Tails (SOT) - from climatology percentiles (sorted) - and ensemble forecast (not sorted) - - Parameters - ---------- - clim: numpy array (nclim, npoints) - Model climatology (percentiles) - ens: numpy array (nens, npoints) - Ensemble forecast - perc: int - Percentile value (typically 10 or 90) - eps: (float) - Epsilon factor for zero values - - Returns - ------- - numpy array (npoints) - SOT values - """ - if not (isinstance(perc, int) or isinstance(perc, np.int64)) or ( - perc < 2 or perc > 98 - ): - raise Exception( - "Percentile value should be and Integer between 2 and 98, is {}".format( - perc - ) - ) - - if clim.shape[0] != 101: - raise Exception( - "Climatology array should contain 101 percentiles, it has {} values".format( - clim.shape - ) - ) - - qc = clim[perc] - # if eps>0, set to zero everything below eps - if eps > 0: - ens = np.where(ens < eps, 0.0, ens) - qc = np.where(qc < eps, 0.0, qc) - - qf = np.percentile(ens, q=perc, axis=0) - if perc > 50: - qc_tail = clim[99] - elif perc < 50: - qc_tail = clim[1] - else: - raise Exception( - "Percentile value to be computed cannot be 50 for sot, has to be in the upper or lower half" - ) - - sot = sot_func(qc_tail, qc, qf, eps=eps) - - return sot - - -def sot_unsorted(clim, ens, perc, eps=-1e4): - """Compute Shift of Tails (SOT) - from climatology percentiles (sorted) - and ensemble forecast (not sorted) - - Parameters - ---------- - clim: numpy array (nclim, npoints) - Model climatology (percentiles) - ens: numpy array (nens, npoints) - Ensemble forecast - perc: int - Percentile value (typically 10 or 90) - eps: (float) - Epsilon factor for zero values - - Returns - ------- - numpy array (npoints) - SOT values - """ - if not (isinstance(perc, int) or isinstance(perc, np.int64)) or ( - perc < 2 or perc > 98 - ): - raise Exception( - "Percentile value should be and Integer between 2 and 98, is {}".format( - perc - ) - ) - - if clim.shape[0] != 101: - raise Exception( - "Climatology array should contain 101 percentiles, it has {} values".format( - clim.shape - ) - ) +from . import array # noqa - if eps > 0: - ens = np.where(ens < eps, 0.0, ens) - clim = np.where(clim < eps, 0.0, clim) - qf = np.percentile(ens, q=perc, axis=0) - qc = np.percentile(clim, q=perc, axis=0) - if perc > 50: - perc_tail = 99 - elif perc < 50: - perc_tail = 1 - else: - raise Exception( - "Percentile value to be computed cannot be 50 for sot, has to be in the upper or lower half" - ) - qc_tail = np.percentile(clim, q=perc_tail, axis=0) +def sot(*args, **kwargs): + return array.sot(*args, **kwargs) - sot = sot_func(qc_tail, qc, qf, eps=eps) - return sot +def sot_unsorted(*args, **kwargs): + return array.sot_unsorted(*args, **kwargs) diff --git a/earthkit/meteo/geo/__init__.py b/earthkit/meteo/geo/__init__.py deleted file mode 100644 index bf69c35..0000000 --- a/earthkit/meteo/geo/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .geo import * # noqa diff --git a/earthkit/meteo/geo/geo.py b/earthkit/meteo/geo/geo.py deleted file mode 100644 index ace6187..0000000 --- a/earthkit/meteo/geo/geo.py +++ /dev/null @@ -1,173 +0,0 @@ -# (C) Copyright 2021- ECMWF. - -# This software is licensed under the terms of the Apache Licence Version 2.0 -# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. - - -import numpy as np - -from earthkit.meteo import constants - - -def distance(lat1, lon1, lat2, lon2): - r"""Computes the spherical distance between two (sets of) points on Earth. - - Parameters - ---------- - lat1: number or ndarray - Latitudes of first set of points (degrees) - lon1: number or ndarray - Longitudes of first set of points (degrees). Must have the same shape as ``lat1``. - lat2: number or ndarray - Latitudes of second set of points (degrees) - lon2: number or ndarray - Longitudes of second set of points (degrees). Must have the same shape as ``lat2``. - - Returns - ------- - number or ndarray - Spherical distance on the surface in Earth (m) - - - The computation is based on the following formula: - - .. math:: - - d = R_{e}\; arccos(sin\phi_{1} sin\phi_{2} + cos\phi_{1} cos\phi_{2} cos\Delta\lambda) - - where :math:`\phi` and :math:`\lambda` stand for the latitude and longitude, - and :math:`R_{e}` is the radius of the Earth - (see :data:`earthkit.meteo.constants.R_earth`). - - As for the input points, the following restrictions apply: - - * at least one of them is a single point - * or both sets of points have the same shape - - """ - lat1 = np.asarray(lat1) - lon1 = np.asarray(lon1) - lat2 = np.asarray(lat2) - lon2 = np.asarray(lon2) - - if lat1.shape != lon1.shape: - raise ValueError( - f"distance: lat1 and lon1 must have the same shape! {lat1.shape} != {lon1.shape}" - ) - - if lat2.shape != lon2.shape: - raise ValueError( - f"distance: lat2 and lon2 must have the same shape! {lat2.shape} != {lon2.shape}" - ) - - if lat1.size > 1 and lat2.size > 1 and lat2.shape != lat1.shape: - raise ValueError( - f"distance: incompatible shapes! lat1={lat1.shape}, lat2={lat2.shape}" - ) - - lat1_rad = lat1 * constants.radian - lat2_rad = lat2 * constants.radian - - x = np.sin(lat1_rad) * np.sin(lat2_rad) + np.cos(lat1_rad) * np.cos( - lat2_rad - ) * np.cos(constants.radian * (lon2 - lon1)) - np.clip(x, -1, 1, out=x) - return np.arccos(x) * constants.R_earth - - -def bearing(lat_ref, lon_ref, lat, lon): - r"""Computes the bearing with respect to a given reference location. - - Parameters - ---------- - lat_ref: number or ndarray - Latitudes of reference points (degrees) - lon_ref: number or ndarray - Longitudes of reference points (degrees). Must have the same shape as ``lat_ref``. - lat: number or ndarray - Latitudes (degrees) - lon: number or ndarray - Longitudes (degrees). Must have the same shape as ``lat``. - - Returns - ------- - number or ndarray - Bearing with respect to the reference points (degree). - - - The **bearing** is the angle between the Northward meridian going through the reference point - and the great circle connecting the reference point and the other point. It is measured in - degrees clockwise from North. If a point is located on the same latitude as the reference - point the bearing is regarded constant: it is either 90° (East) or 270° (West). If the point - is co-located with the reference point the bearing is set to np.nan. - - As for the reference and other points, the following restrictions apply: - - * at least one of them is a single point - * or both the reference and other points have the same shape - - """ - lat_ref = np.asarray(lat_ref) - lon_ref = np.asarray(lon_ref) - lat = np.asarray(lat) - lon = np.asarray(lon) - - if lat_ref.shape != lon_ref.shape: - raise ValueError( - f"bearing: lat_ref and lon_ref must have the same shape! {lat_ref.shape} != {lon_ref.shape}" - ) - - if lat.shape != lon.shape: - raise ValueError( - f"bearing: lat and lon must have the same shape! {lat.shape} != {lon.shape}" - ) - - if lat_ref.size > 1 and lat.size > 1 and lat_ref.shape != lat.shape: - raise ValueError( - f"bearing: incompatible shapes! lat_ref={lat_ref.shape}, lat={lat.shape}" - ) - - eps = 1e-9 - if lat_ref.shape: - br = np.full(lat_ref.shape, np.nan) - else: - br = np.full(lat.shape, np.nan) - - # computes longitude difference - d_lon = lon.copy() - d_lon[lon > 180.0] -= 360.0 - d_lon = (d_lon - lon_ref) * constants.radian - - lat_ref_rad = lat_ref * constants.radian - lat_rad = lat * constants.radian - - # the bearing is constant on the same latitude - mask = np.fabs(lat - lat_ref) < eps - mask_1 = mask & (np.fabs(d_lon) >= eps) - br[mask_1 & (d_lon <= 0.0)] = 270.0 - br[mask_1 & (d_lon > 0.0)] = 90.0 - - # deals with different latitudes - mask = ~mask - if mask.shape == lat_ref.shape: - lat_ref_rad = lat_ref_rad[mask] - if mask.shape == lat.shape: - lat_rad = lat_rad[mask] - - # bearing in degrees, x axis points to East, anti-clockwise - br_m = np.arctan2( - np.cos(lat_ref_rad) * np.sin(lat_rad) - - np.sin(lat_ref_rad) * np.cos(lat_rad) * np.cos(d_lon[mask]), - np.sin(d_lon[mask]) * np.cos(lat_ref_rad), - ) - - # transforms to the required coordinate system: x axis points to North, clockwise - br_m = (np.pi / 2.0) - br_m - br_m[br_m < 0] += 2.0 * np.pi - br_m *= constants.degree - - br[mask] = br_m - return br diff --git a/earthkit/meteo/score/__init__.py b/earthkit/meteo/score/__init__.py index ac48af7..f11e64d 100644 --- a/earthkit/meteo/score/__init__.py +++ b/earthkit/meteo/score/__init__.py @@ -1 +1,18 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Verification functions. + +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. +""" + from .crps import * # noqa diff --git a/tests/test_solar.py b/earthkit/meteo/score/array/__init__.py similarity index 75% rename from tests/test_solar.py rename to earthkit/meteo/score/array/__init__.py index 5643a00..ba1acf9 100644 --- a/tests/test_solar.py +++ b/earthkit/meteo/score/array/__init__.py @@ -1,13 +1,14 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# -# import earthkit.meteo - +""" +Verification functions operating on numpy arrays. +""" -def test_solar() -> None: - pass +from .crps import * # noqa diff --git a/earthkit/meteo/score/array/crps.py b/earthkit/meteo/score/array/crps.py new file mode 100644 index 0000000..9595094 --- /dev/null +++ b/earthkit/meteo/score/array/crps.py @@ -0,0 +1,60 @@ +# +# # (C) Copyright 2021 ECMWF. +# # +# # This software is licensed under the terms of the Apache Licence Version 2.0 +# # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# # In applying this licence, ECMWF does not waive the privileges and immunities +# # granted to it by virtue of its status as an intergovernmental organisation +# # nor does it submit to any jurisdiction. +# ## +import numpy as np + + +def crps(x, y): + """Computes Continuous Ranked Probability Score (CRPS). + + Parameters + ---------- + x: numpy array (n_ens, n_points) + Ensemble forecast + y: numpy array (n_points) + Observation/analysis + + Returns + ------- + numpy array (n_points) + CRPS values + + + The method is described in [Hersbach2000]_. + """ + # first sort ensemble + x.sort(axis=0) + + # construct alpha and beta, size nens+1 + n_ens = x.shape[0] + shape = (n_ens + 1,) + x.shape[1:] + alpha = np.zeros(shape) + beta = np.zeros(shape) + + # x[i+1]-x[i] and x[i]-y[i] arrays + diffxy = x - y.reshape(1, *(y.shape)) + diffxx = x[1:] - x[:-1] # x[i+1]-x[i], size ens-1 + + # if i == 0 + alpha[0] = 0 + beta[0] = np.fmax(diffxy[0], 0) # x(0)-y + # if i == n_ens + alpha[-1] = np.fmax(-diffxy[-1], 0) # y-x(n) + beta[-1] = 0 + # else + alpha[1:-1] = np.fmin( + diffxx, np.fmax(-diffxy[:-1], 0) + ) # x(i+1)-x(i) or y-x(i) or 0 + beta[1:-1] = np.fmin(diffxx, np.fmax(diffxy[1:], 0)) # 0 or x(i+1)-y or x(i+1)-x(i) + + # compute crps + p_exp = (np.arange(n_ens + 1) / float(n_ens)).reshape(n_ens + 1, *([1] * y.ndim)) + crps = np.sum(alpha * (p_exp**2) + beta * ((1 - p_exp) ** 2), axis=0) + + return crps diff --git a/earthkit/meteo/score/crps.py b/earthkit/meteo/score/crps.py index c6174a1..4423b84 100644 --- a/earthkit/meteo/score/crps.py +++ b/earthkit/meteo/score/crps.py @@ -1,51 +1,14 @@ -import numpy as np +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# +from . import array -def crps(x, y): - """Computes Continuous Ranked Probability Score (CRPS). - Parameters - ---------- - x: numpy array (n_ens, n_points) - Ensemble forecast - y: numpy array (n_points) - Observation/analysis - - Returns - ------- - numpy array (n_points) - CRPS values - - - The method is described in [Hersbach2000]_. - """ - # first sort ensemble - x.sort(axis=0) - - # construct alpha and beta, size nens+1 - n_ens = x.shape[0] - shape = (n_ens + 1,) + x.shape[1:] - alpha = np.zeros(shape) - beta = np.zeros(shape) - - # x[i+1]-x[i] and x[i]-y[i] arrays - diffxy = x - y.reshape(1, *(y.shape)) - diffxx = x[1:] - x[:-1] # x[i+1]-x[i], size ens-1 - - # if i == 0 - alpha[0] = 0 - beta[0] = np.fmax(diffxy[0], 0) # x(0)-y - # if i == n_ens - alpha[-1] = np.fmax(-diffxy[-1], 0) # y-x(n) - beta[-1] = 0 - # else - alpha[1:-1] = np.fmin( - diffxx, np.fmax(-diffxy[:-1], 0) - ) # x(i+1)-x(i) or y-x(i) or 0 - beta[1:-1] = np.fmin(diffxx, np.fmax(diffxy[1:], 0)) # 0 or x(i+1)-y or x(i+1)-x(i) - - # compute crps - p_exp = (np.arange(n_ens + 1) / float(n_ens)).reshape(n_ens + 1, *([1] * y.ndim)) - crps = np.sum(alpha * (p_exp**2) + beta * ((1 - p_exp) ** 2), axis=0) - - return crps +def crps(*args, **kwargs): + return array.crps(*args, **kwargs) diff --git a/earthkit/meteo/solar/__init__.py b/earthkit/meteo/solar/__init__.py index 212bb87..9bad33c 100644 --- a/earthkit/meteo/solar/__init__.py +++ b/earthkit/meteo/solar/__init__.py @@ -1,240 +1,19 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# -import datetime - -import numpy as np - -DAYS_PER_YEAR = 365.25 """ -WARNING: This code is incomplete and not tested. DO NOT USE. -""" - - -def julian_day(date): - delta = date - datetime.datetime(date.year, 1, 1) - return delta.days + delta.seconds / 86400.0 - - -def solar_declination_angle(date): - angle = julian_day(date) / DAYS_PER_YEAR * np.pi * 2 - - # declination in [degrees] - declination = ( - 0.396372 - - 22.91327 * np.cos(angle) - + 4.025430 * np.sin(angle) - - 0.387205 * np.cos(2 * angle) - + 0.051967 * np.sin(2 * angle) - - 0.154527 * np.cos(3 * angle) - + 0.084798 * np.sin(3 * angle) - ) - # time correction in [ h.degrees ] - time_correction = ( - 0.004297 - + 0.107029 * np.cos(angle) - - 1.837877 * np.sin(angle) - - 0.837378 * np.cos(2 * angle) - - 2.340475 * np.sin(2 * angle) - ) - return declination, time_correction - - -def cos_solar_zenith_angle(date, latitudes, longitudes): - """Cosine of solar zenith angle. - - Parameters - ---------- - date: datetime.datetime - Date - lat: float array - Latitude [degrees] - lon: float array - Longitude [degrees] - - Returns - ------- - float array - Cosine of the solar zenith angle (all values, including negatives) - [Hogan_and_Hirahara2015]_. See also: - http://answers.google.com/answers/threadview/id/782886.html - - """ - # declination angle + time correction for solar angle - declination, time_correction = solar_declination_angle(date) - - # solar_declination_angle returns degrees - declination = np.deg2rad(declination) - - latitudes = np.deg2rad(latitudes) - - sindec_sinlat = np.sin(declination) * np.sin(latitudes) - cosdec_coslat = np.cos(declination) * np.cos(latitudes) - - # solar hour angle [h.deg] - solar_angle = np.deg2rad((date.hour - 12) * 15 + longitudes + time_correction) - zenith_angle = sindec_sinlat + cosdec_coslat * np.cos(solar_angle) - - # Clip negative values - return np.clip(zenith_angle, 0, None) - +Solar computation functions. -def _integrate( - func, - begin_date, - end_date, - latitudes, - longitudes, - *, - intervals_per_hour=1, - integration_order=3, -): - # Gauss-Integration coefficients - if integration_order == 3: # default, good speed and accuracy (3 points) - E = np.array([-np.sqrt(3.0 / 5.0), 0.0, np.sqrt(3.0 / 5.0)]) - W = np.array([5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0]) - elif integration_order == 1: # fastest, worse accuracy (1 point) - E = np.array([0.0]) - W = np.array([2.0]) - elif integration_order == 2: # faster, less accurate (2 points) - E = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)]) - W = np.array([1.0, 1.0]) - elif integration_order == 4: # slower, more accurate (4 points) - E = np.array( - [ - -np.sqrt(3.0 / 7.0 + 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), - -np.sqrt(3.0 / 7.0 - 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), - np.sqrt(3.0 / 7.0 - 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), - np.sqrt(3.0 / 7.0 + 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), - ] - ) - W = np.array( - [ - (18 - np.sqrt(30)) / 36, - (18 + np.sqrt(30)) / 36, - (18 + np.sqrt(30)) / 36, - (18 - np.sqrt(30)) / 36, - ] - ) - else: - raise ValueError("Invalid integration order %d", integration_order) - - assert intervals_per_hour > 0 - assert end_date > begin_date - - date = begin_date - interval_size_hours = (end_date - begin_date).total_seconds() / 3600.0 - - nsplits = int(interval_size_hours * intervals_per_hour + 0.5) - - assert nsplits > 0 - - time_steps = np.linspace(0, interval_size_hours, num=nsplits + 1) - - integral = np.zeros_like(latitudes) - for s in range(len(time_steps) - 1): - ti = time_steps[s] - tf = time_steps[s + 1] - - deltat = tf - ti - jacob = deltat / 2.0 - - w = jacob * W - w /= interval_size_hours # average of integral - t = jacob * E - t += (tf + ti) / 2.0 - - for n in range(len(w)): - integral += w[n] * func( - date + datetime.timedelta(hours=t[n]), - latitudes, - longitudes, - ) - - return integral - - -def cos_solar_zenith_angle_integrated( - begin_date, - end_date, - latitudes, - longitudes, - *, - intervals_per_hour=1, - integration_order=3, -): - """Average of solar zenith angle based on numerical integration. - - Parameters - ---------- - begin_date: datetime.datetime - end_date: datetime.datetime - lat: int darray - Latitude [degrees]. - lon: int darray - Longitude [degrees]. - tbegin: int - Offset in hours from forecast time to begin of time interval for integration. - tend: int - Offset in hours from forecast time to end of time interval for integration. - intervals_per_hour: int - Number of time integrations per hour. - integration order: int - Order of gauss integration, valid = (1, 2, 3, 4) - - Returns - ------- - float array - Average of cosine of the solar zenith angle during interval [degrees]. Based on - numerical integration using the 3 point - `Gauss integration `_ rule. - [Hogan_and_Hirahara2015]_, [Biricombe2022]_ - - """ - return _integrate( - cos_solar_zenith_angle, - begin_date, - end_date, - latitudes, - longitudes, - intervals_per_hour=intervals_per_hour, - integration_order=integration_order, - ) - - -def incoming_solar_radiation(date): - # To be replaced with improved formula - (a, b) = (165120.0, 4892416.0) - angle = julian_day(date) / DAYS_PER_YEAR * np.pi * 2 - return np.cos(angle) * a + b - - -def toa_incident_solar_radiation( - begin_date, - end_date, - latitudes, - longitudes, - *, - intervals_per_hour=1, - integration_order=3, -): - def func(date, latitudes, longitudes): - isr = incoming_solar_radiation(date) - csza = cos_solar_zenith_angle(date, latitudes, longitudes) - return isr * csza +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets* +""" - return _integrate( - func, - begin_date, - end_date, - latitudes, - longitudes, - intervals_per_hour=intervals_per_hour, - integration_order=integration_order, - ) +from .solar import * # noqa diff --git a/earthkit/meteo/solar/array/__init__.py b/earthkit/meteo/solar/array/__init__.py new file mode 100644 index 0000000..da2e210 --- /dev/null +++ b/earthkit/meteo/solar/array/__init__.py @@ -0,0 +1,14 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Solar computation functions operating on numpy arrays. +""" + +from .solar import * # noqa diff --git a/earthkit/meteo/solar/array/solar.py b/earthkit/meteo/solar/array/solar.py new file mode 100644 index 0000000..21848b4 --- /dev/null +++ b/earthkit/meteo/solar/array/solar.py @@ -0,0 +1,241 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import datetime + +import numpy as np + +DAYS_PER_YEAR = 365.25 + +""" +WARNING: This code is incomplete and not tested. DO NOT USE. +""" + + +def julian_day(date): + delta = date - datetime.datetime(date.year, 1, 1) + return delta.days + delta.seconds / 86400.0 + + +def solar_declination_angle(date): + angle = julian_day(date) / DAYS_PER_YEAR * np.pi * 2 + + # declination in [degrees] + declination = ( + 0.396372 + - 22.91327 * np.cos(angle) + + 4.025430 * np.sin(angle) + - 0.387205 * np.cos(2 * angle) + + 0.051967 * np.sin(2 * angle) + - 0.154527 * np.cos(3 * angle) + + 0.084798 * np.sin(3 * angle) + ) + # time correction in [ h.degrees ] + time_correction = ( + 0.004297 + + 0.107029 * np.cos(angle) + - 1.837877 * np.sin(angle) + - 0.837378 * np.cos(2 * angle) + - 2.340475 * np.sin(2 * angle) + ) + return declination, time_correction + + +def cos_solar_zenith_angle(date, latitudes, longitudes): + """Cosine of solar zenith angle. + + Parameters + ---------- + date: datetime.datetime + Date + lat: float array + Latitude [degrees] + lon: float array + Longitude [degrees] + + Returns + ------- + float array + Cosine of the solar zenith angle (all values, including negatives) + [Hogan_and_Hirahara2015]_. See also: + http://answers.google.com/answers/threadview/id/782886.html + + """ + # declination angle + time correction for solar angle + declination, time_correction = solar_declination_angle(date) + + # solar_declination_angle returns degrees + declination = np.deg2rad(declination) + + latitudes = np.deg2rad(latitudes) + + sindec_sinlat = np.sin(declination) * np.sin(latitudes) + cosdec_coslat = np.cos(declination) * np.cos(latitudes) + + # solar hour angle [h.deg] + solar_angle = np.deg2rad((date.hour - 12) * 15 + longitudes + time_correction) + zenith_angle = sindec_sinlat + cosdec_coslat * np.cos(solar_angle) + + # Clip negative values + return np.clip(zenith_angle, 0, None) + + +def _integrate( + func, + begin_date, + end_date, + latitudes, + longitudes, + *, + intervals_per_hour=1, + integration_order=3, +): + # Gauss-Integration coefficients + if integration_order == 3: # default, good speed and accuracy (3 points) + E = np.array([-np.sqrt(3.0 / 5.0), 0.0, np.sqrt(3.0 / 5.0)]) + W = np.array([5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0]) + elif integration_order == 1: # fastest, worse accuracy (1 point) + E = np.array([0.0]) + W = np.array([2.0]) + elif integration_order == 2: # faster, less accurate (2 points) + E = np.array([-1.0 / np.sqrt(3.0), 1.0 / np.sqrt(3.0)]) + W = np.array([1.0, 1.0]) + elif integration_order == 4: # slower, more accurate (4 points) + E = np.array( + [ + -np.sqrt(3.0 / 7.0 + 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), + -np.sqrt(3.0 / 7.0 - 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), + np.sqrt(3.0 / 7.0 - 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), + np.sqrt(3.0 / 7.0 + 2.0 / 7.0 * np.sqrt(6.0 / 5.0)), + ] + ) + W = np.array( + [ + (18 - np.sqrt(30)) / 36, + (18 + np.sqrt(30)) / 36, + (18 + np.sqrt(30)) / 36, + (18 - np.sqrt(30)) / 36, + ] + ) + else: + raise ValueError("Invalid integration order %d", integration_order) + + assert intervals_per_hour > 0 + assert end_date > begin_date + + date = begin_date + interval_size_hours = (end_date - begin_date).total_seconds() / 3600.0 + + nsplits = int(interval_size_hours * intervals_per_hour + 0.5) + + assert nsplits > 0 + + time_steps = np.linspace(0, interval_size_hours, num=nsplits + 1) + + integral = np.zeros_like(latitudes) + for s in range(len(time_steps) - 1): + ti = time_steps[s] + tf = time_steps[s + 1] + + deltat = tf - ti + jacob = deltat / 2.0 + + w = jacob * W + w /= interval_size_hours # average of integral + t = jacob * E + t += (tf + ti) / 2.0 + + for n in range(len(w)): + integral += w[n] * func( + date + datetime.timedelta(hours=t[n]), + latitudes, + longitudes, + ) + + return integral + + +def cos_solar_zenith_angle_integrated( + begin_date, + end_date, + latitudes, + longitudes, + *, + intervals_per_hour=1, + integration_order=3, +): + """Average of solar zenith angle based on numerical integration. + + Parameters + ---------- + begin_date: datetime.datetime + end_date: datetime.datetime + lat: int darray + Latitude [degrees]. + lon: int darray + Longitude [degrees]. + tbegin: int + Offset in hours from forecast time to begin of time interval for integration. + tend: int + Offset in hours from forecast time to end of time interval for integration. + intervals_per_hour: int + Number of time integrations per hour. + integration order: int + Order of gauss integration, valid = (1, 2, 3, 4) + + Returns + ------- + float array + Average of cosine of the solar zenith angle during interval [degrees]. Based on + numerical integration using the 3 point + `Gauss integration `_ rule. + [Hogan_and_Hirahara2015]_, [Biricombe2022]_ + + """ + return _integrate( + cos_solar_zenith_angle, + begin_date, + end_date, + latitudes, + longitudes, + intervals_per_hour=intervals_per_hour, + integration_order=integration_order, + ) + + +def incoming_solar_radiation(date): + # To be replaced with improved formula + (a, b) = (165120.0, 4892416.0) + angle = julian_day(date) / DAYS_PER_YEAR * np.pi * 2 + return np.cos(angle) * a + b + + +def toa_incident_solar_radiation( + begin_date, + end_date, + latitudes, + longitudes, + *, + intervals_per_hour=1, + integration_order=3, +): + def func(date, latitudes, longitudes): + isr = incoming_solar_radiation(date) + csza = cos_solar_zenith_angle(date, latitudes, longitudes) + return isr * csza + + return _integrate( + func, + begin_date, + end_date, + latitudes, + longitudes, + intervals_per_hour=intervals_per_hour, + integration_order=integration_order, + ) diff --git a/earthkit/meteo/solar/solar.py b/earthkit/meteo/solar/solar.py new file mode 100644 index 0000000..77c754b --- /dev/null +++ b/earthkit/meteo/solar/solar.py @@ -0,0 +1,34 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +from . import array + + +def julian_day(*args, **kwargs): + return array.julian_day(*args, **kwargs) + + +def solar_declination_angle(*args, **kwargs): + return array.solar_declination_angle(*args, **kwargs) + + +def cos_solar_zenith_angle(*args, **kwargs): + return array.cos_solar_zenith_angle(*args, **kwargs) + + +def cos_solar_zenith_angle_integrated(*args, **kwargs): + return array.cos_solar_zenith_angle_integrated(*args, **kwargs) + + +def incoming_solar_radiation(*args, **kwargs): + return array.incoming_solar_radiation(*args, **kwargs) + + +def toa_incident_solar_radiation(*args, **kwargs): + return array.toa_incident_solar_radiation(*args, **kwargs) diff --git a/earthkit/meteo/stats/__init__.py b/earthkit/meteo/stats/__init__.py index 2482ebe..af8730a 100644 --- a/earthkit/meteo/stats/__init__.py +++ b/earthkit/meteo/stats/__init__.py @@ -1,2 +1,19 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Statistical functions. + +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. +""" + from .numpy_extended import * # noqa from .quantiles import * # noqa diff --git a/earthkit/meteo/stats/array/__init__.py b/earthkit/meteo/stats/array/__init__.py new file mode 100644 index 0000000..e43761a --- /dev/null +++ b/earthkit/meteo/stats/array/__init__.py @@ -0,0 +1,15 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Statistical functions operating on numpy arrays. +""" + +from .numpy_extended import * # noqa +from .quantiles import * # noqa diff --git a/earthkit/meteo/stats/array/numpy_extended.py b/earthkit/meteo/stats/array/numpy_extended.py new file mode 100644 index 0000000..7f78884 --- /dev/null +++ b/earthkit/meteo/stats/array/numpy_extended.py @@ -0,0 +1,55 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + + +def nanaverage(data, weights=None, **kwargs): + """A merge of the functionality of np.nanmean and np.average. + + Parameters + ---------- + data : numpy array + weights: Weights to apply to the data for averaging. + Weights will be normalised and must correspond to the + shape of the numpy data array and axis/axes that is/are + averaged over. + axis: axis/axes to compute the nanaverage over. + kwargs: any other np.nansum kwargs + + Returns + ------- + numpy array + mean of data (along axis) where nan-values are ignored + and weights applied if provided. + """ + if weights is not None: + # set weights to nan where data is nan: + this_weights = np.ones(data.shape) * weights + this_weights[np.isnan(data)] = np.nan + # Weights must be scaled to the sum of valid + # weights for each relevant axis: + this_denom = np.nansum(this_weights, **kwargs) + # If averaging over an axis then we must add dummy + # dimension[s] to the denominator to make compatible + # with the weights. + if kwargs.get("axis", None) is not None: + reshape = list(this_weights.shape) + reshape[kwargs.get("axis")] = 1 + this_denom = this_denom.reshape(reshape) + + # Scale weights to mean of valid weights: + this_weights = this_weights / this_denom + # Apply weights to data: + nanaverage = np.nansum(data * this_weights, **kwargs) + else: + # If no weights, then nanmean will suffice + nanaverage = np.nanmean(data, **kwargs) + + return nanaverage diff --git a/earthkit/meteo/stats/array/quantiles.py b/earthkit/meteo/stats/array/quantiles.py new file mode 100644 index 0000000..33e1cd6 --- /dev/null +++ b/earthkit/meteo/stats/array/quantiles.py @@ -0,0 +1,77 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +from typing import Iterable, List, Union + +import numpy as np + + +def iter_quantiles( + arr: np.ndarray, + which: Union[int, List[float]] = 100, + axis: int = 0, + method: str = "sort", +) -> Iterable[np.ndarray]: + """Iterate over the quantiles of a large array + + Parameters + ---------- + arr: numpy array + Data array + which: int or list of floats + List of quantiles to compute, e.g. `[0., 0.25, 0.5, 0.75, 1.]`, or + number of evenly-spaced intervals (e.g. 100 for percentiles). + axis: int + Axis along which to compute the quantiles + method: 'sort', 'numpy_bulk', 'numpy' + Method of computing the quantiles: + * sort: sort `arr` in place, then interpolates the quantiles one by one + * numpy_bulk: compute all the quantiles at once using `numpy.quantile` + * numpy: compute the quantiles one by one using `numpy.quantile` + + Returns + ------- + Iterable[numpy array] + Quantiles, in increasing order if `which` is an `int`, otherwise in the order specified + """ + if method not in ("sort", "numpy_bulk", "numpy"): + raise ValueError( + f"Invalid method {method!r}, expected 'sort', 'numpy_bulk', or 'numpy'" + ) + + if isinstance(which, int): + n = which + qs = np.linspace(0.0, 1.0, n + 1) + else: + qs = np.asarray(which) + + if method == "numpy_bulk": + quantiles = np.quantile(arr, qs, axis=axis) + yield from quantiles + return + + if method == "sort": + arr = np.asarray(arr) + arr.sort(axis=axis) + + for q in qs: + if method == "numpy": + yield np.quantile(arr, q, axis=axis) + + elif method == "sort": + m = arr.shape[axis] + f = (m - 1) * q + j = int(f) + x = f - j + quantile = arr.take(j, axis=axis) + quantile *= 1 - x + tmp = arr.take(min(j + 1, m - 1), axis=axis) + tmp *= x + quantile += tmp + yield quantile diff --git a/earthkit/meteo/stats/numpy_extended.py b/earthkit/meteo/stats/numpy_extended.py index 119e858..3a8fdec 100644 --- a/earthkit/meteo/stats/numpy_extended.py +++ b/earthkit/meteo/stats/numpy_extended.py @@ -1,46 +1,14 @@ -import numpy as np +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# +from . import array -def nanaverage(data, weights=None, **kwargs): - """A merge of the functionality of np.nanmean and np.average. - Parameters - ---------- - data : numpy array - weights: Weights to apply to the data for averaging. - Weights will be normalised and must correspond to the - shape of the numpy data array and axis/axes that is/are - averaged over. - axis: axis/axes to compute the nanaverage over. - kwargs: any other np.nansum kwargs - - Returns - ------- - numpy array - mean of data (along axis) where nan-values are ignored - and weights applied if provided. - """ - if weights is not None: - # set weights to nan where data is nan: - this_weights = np.ones(data.shape) * weights - this_weights[np.isnan(data)] = np.nan - # Weights must be scaled to the sum of valid - # weights for each relevant axis: - this_denom = np.nansum(this_weights, **kwargs) - # If averaging over an axis then we must add dummy - # dimension[s] to the denominator to make compatible - # with the weights. - if kwargs.get("axis", None) is not None: - reshape = list(this_weights.shape) - reshape[kwargs.get("axis")] = 1 - this_denom = this_denom.reshape(reshape) - - # Scale weights to mean of valid weights: - this_weights = this_weights / this_denom - # Apply weights to data: - nanaverage = np.nansum(data * this_weights, **kwargs) - else: - # If no weights, then nanmean will suffice - nanaverage = np.nanmean(data, **kwargs) - - return nanaverage +def nanaverage(*args, **kwargs): + return array.nanaverage(*args, **kwargs) diff --git a/earthkit/meteo/stats/quantiles.py b/earthkit/meteo/stats/quantiles.py index f31caab..4f607d0 100644 --- a/earthkit/meteo/stats/quantiles.py +++ b/earthkit/meteo/stats/quantiles.py @@ -1,68 +1,14 @@ -from typing import Iterable, List, Union +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# -import numpy as np +from . import array -def iter_quantiles( - arr: np.ndarray, - which: Union[int, List[float]] = 100, - axis: int = 0, - method: str = "sort", -) -> Iterable[np.ndarray]: - """Iterate over the quantiles of a large array - - Parameters - ---------- - arr: numpy array - Data array - which: int or list of floats - List of quantiles to compute, e.g. `[0., 0.25, 0.5, 0.75, 1.]`, or - number of evenly-spaced intervals (e.g. 100 for percentiles). - axis: int - Axis along which to compute the quantiles - method: 'sort', 'numpy_bulk', 'numpy' - Method of computing the quantiles: - * sort: sort `arr` in place, then interpolates the quantiles one by one - * numpy_bulk: compute all the quantiles at once using `numpy.quantile` - * numpy: compute the quantiles one by one using `numpy.quantile` - - Returns - ------- - Iterable[numpy array] - Quantiles, in increasing order if `which` is an `int`, otherwise in the order specified - """ - if method not in ("sort", "numpy_bulk", "numpy"): - raise ValueError( - f"Invalid method {method!r}, expected 'sort', 'numpy_bulk', or 'numpy'" - ) - - if isinstance(which, int): - n = which - qs = np.linspace(0.0, 1.0, n + 1) - else: - qs = np.asarray(which) - - if method == "numpy_bulk": - quantiles = np.quantile(arr, qs, axis=axis) - yield from quantiles - return - - if method == "sort": - arr = np.asarray(arr) - arr.sort(axis=axis) - - for q in qs: - if method == "numpy": - yield np.quantile(arr, q, axis=axis) - - elif method == "sort": - m = arr.shape[axis] - f = (m - 1) * q - j = int(f) - x = f - j - quantile = arr.take(j, axis=axis) - quantile *= 1 - x - tmp = arr.take(min(j + 1, m - 1), axis=axis) - tmp *= x - quantile += tmp - yield quantile +def iter_quantiles(*args, **kwargs): + return array.iter_quantiles(*args, **kwargs) diff --git a/earthkit/meteo/thermo/__init__.py b/earthkit/meteo/thermo/__init__.py index 8842a88..270a840 100644 --- a/earthkit/meteo/thermo/__init__.py +++ b/earthkit/meteo/thermo/__init__.py @@ -1 +1,19 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Thermodynamic functions. + +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. +""" + + from .thermo import * # noqa diff --git a/earthkit/meteo/thermo/array/__init__.py b/earthkit/meteo/thermo/array/__init__.py new file mode 100644 index 0000000..5507b1b --- /dev/null +++ b/earthkit/meteo/thermo/array/__init__.py @@ -0,0 +1,14 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Thermodynamic functions operating on numpy arrays. +""" + +from .thermo import * # noqa diff --git a/earthkit/meteo/thermo/array/thermo.py b/earthkit/meteo/thermo/array/thermo.py new file mode 100644 index 0000000..f0c2b7a --- /dev/null +++ b/earthkit/meteo/thermo/array/thermo.py @@ -0,0 +1,1791 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + + +import numpy as np + +from earthkit.meteo import constants + + +def celsius_to_kelvin(t): + """Converts temperature values from Celsius to Kelvin. + + Parameters + ---------- + t : number or ndarray + Temperature in Celsius units + + Returns + ------- + number or ndarray + Temperature in Kelvin units + + """ + return t + constants.T0 + + +def kelvin_to_celsius(t): + """Converts temperature values from Kelvin to Celsius. + + Parameters + ---------- + t : number or ndarray + Temperature in Kelvin units + + Returns + ------- + number or ndarray + Temperature in Celsius units + + """ + return t - constants.T0 + + +def specific_humidity_from_mixing_ratio(w): + r"""Computes the specific humidity from mixing ratio. + + Parameters + ---------- + w : number or ndarray + Mixing ratio (kg/kg) + + Returns + ------- + number or ndarray + Specific humidity (kg/kg) + + + The result is the specific humidity in kg/kg units. The computation is based on + the following definition [Wallace2006]_: + + .. math:: + + q = \frac {w}{1+w} + + """ + return w / (1 + w) + + +def mixing_ratio_from_specific_humidity(q): + r"""Computes the mixing ratio from specific humidity. + + Parameters + ---------- + q : number or ndarray + Specific humidity (kg/kg) + + Returns + ------- + number or ndarray + Mixing ratio (kg/kg) + + + The result is the mixing ratio in kg/kg units. The computation is based on + the following definition [Wallace2006]_: + + .. math:: + + w = \frac {q}{1-q} + + """ + return q / (1 - q) + + +def vapour_pressure_from_specific_humidity(q, p): + r"""Computes the vapour pressure from specific humidity. + + Parameters + ---------- + q: number or ndarray + Specific humidity (kg/kg) + p: number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Vapour pressure (Pa) + + + The computation is based on the following formula [Wallace2006]_: + + .. math:: + + e = \frac{p\;q}{\epsilon\; (1 + q(\frac{1}{\epsilon} -1 ))} + + with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + c = constants.epsilon * (1 / constants.epsilon - 1) + return (p * q) / (constants.epsilon + c * q) + + +def vapour_pressure_from_mixing_ratio(w, p): + r"""Computes the vapour pressure from mixing ratio. + + Parameters + ---------- + w: number or ndarray + Mixing ratio (kg/kg) + p: number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Vapour pressure (Pa) + + + The computation is based on the following formula: + + .. math:: + + e = \frac{p\;w}{\epsilon + w} + + with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + return (p * w) / (constants.epsilon + w) + + +def specific_humidity_from_vapour_pressure(e, p, eps=1e-4): + r"""Computes the specific humidity from vapour pressure. + + Parameters + ---------- + e: number or ndarray + Vapour pressure (Pa) + p: number or ndarray + Pressure (Pa) + eps: number + Where p - e < ``eps`` np.nan is returned. + + Returns + ------- + number or ndarray + Specific humidity (kg/kg) + + + The computation is based on the following formula: + + .. math:: + + q = \frac{\epsilon e}{p + e(\epsilon-1)} + + with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + if eps <= 0: + raise ValueError( + f"specific_humidity_from_vapour_pressure(): eps={eps} must be > 0" + ) + + v = np.asarray(p + (constants.epsilon - 1) * e) + v[np.asarray(p - e) < eps] = np.nan + return constants.epsilon * e / v + + +def mixing_ratio_from_vapour_pressure(e, p, eps=1e-4): + r"""Computes the mixing ratio from vapour pressure. + + Parameters + ---------- + e: number or ndarray + Vapour pressure (Pa) + p: number or ndarray + Pressure (Pa) + eps: number + Where p - e < ``eps`` np.nan is returned. + + Returns + ------- + number or ndarray + Mixing ratio (kg/kg). + + + The computation is based on the following formula: + + .. math:: + + w = \frac{\epsilon e}{p - e} + + with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + if eps <= 0: + raise ValueError(f"mixing_ratio_from_vapour_pressure(): eps={eps} must be > 0") + + v = np.asarray(p - e) + v[v < eps] = np.nan + return constants.epsilon * e / v + + +class _EsComp: + c1 = 611.21 + c3w = 17.502 + c4w = 32.19 + c3i = 22.587 + c4i = -0.7 + t0 = 273.16 + ti = t0 - 23 + PHASES = ["mixed", "water", "ice"] + + def _check_phase(self, phase): + if phase not in _EsComp.PHASES: + raise ValueError( + f"saturation_vapour_pressure(): invalid phase={phase}! Allowed values = {_EsComp.PHASES}" + ) + return True + + def compute_es(self, t, phase): + self._check_phase(phase) + if phase == "mixed": + return self._es_mixed(t) + elif phase == "water": + return self._es_water(t) + elif phase == "ice": + return self._es_ice(t) + + def compute_slope(self, t, phase): + self._check_phase(phase) + if phase == "mixed": + return self._es_mixed_slope(t) + elif phase == "water": + return self._es_water_slope(t) + elif phase == "ice": + return self._es_ice_slope(t) + + def t_from_es(self, es): + v = np.log(es / self.c1) + return (v * self.c4w - self.c3w * self.t0) / (v - self.c3w) + + def _es_water(self, t): + return self.c1 * np.exp(self.c3w * (t - self.t0) / (t - self.c4w)) + + def _es_ice(self, t): + return self.c1 * np.exp(self.c3i * (t - self.t0) / (t - self.c4i)) + + def _es_mixed(self, t): + # Fraction of liquid water (=alpha): + # t <= ti => alpha=0 + # t > ti and t < t0 => alpha=(t-ti)/(t0-ti))^2 + # t >= t0 => alpha=1 + # + # svp is interpolated between the ice and water phases: + # svp = alpha * es_water + (1.0 - alpha) * es_ice + + t = np.asarray(t) + svp = np.zeros(t.shape) + + # ice range + i_mask = t <= self.ti + svp[i_mask] = self._es_ice(t[i_mask]) + + # water range + w_mask = t >= self.t0 + svp[w_mask] = self._es_water(t[w_mask]) + + # mixed range + m_mask = ~(i_mask | w_mask) + alpha = np.square(t[m_mask] - self.ti) / np.square(self.t0 - self.ti) + svp[m_mask] = alpha * self._es_water(t[m_mask]) + (1.0 - alpha) * self._es_ice( + t[m_mask] + ) + return svp + + def _es_water_slope(self, t): + return ( + self._es_water(t) + * (self.c3w * (self.t0 - self.c4w)) + / np.square(t - self.c4w) + ) + + def _es_ice_slope(self, t): + return ( + self._es_ice(t) + * (self.c3i * (self.t0 - self.c4i)) + / np.square(t - self.c4i) + ) + + def _es_mixed_slope(self, t): + t = np.asarray(t) + d_svp = np.zeros(t.shape) + + # ice range + i_mask = t <= self.ti + d_svp[i_mask] = self._es_ice_slope(t[i_mask]) + + # water range + w_mask = t >= self.t0 + d_svp[w_mask] = self._es_water_slope(t[w_mask]) + + # mixed range + m_mask = ~(i_mask | w_mask) + alpha = np.square(t[m_mask] - self.ti) / np.square(self.t0 - self.ti) + d_alpha = (2.0 / np.square(self.t0 - self.ti)) * (t[m_mask] - self.ti) + t_m = t[m_mask] + d_svp[m_mask] = ( + d_alpha * self._es_water(t_m) + + alpha * self._es_water_slope(t_m) + - d_alpha * self._es_ice(t_m) + + (1.0 - alpha) * self._es_ice_slope(t_m) + ) + return d_svp + + +def saturation_vapour_pressure(t, phase="mixed"): + r"""Computes the saturation vapour pressure from temperature with respect to a phase. + + Parameters + ---------- + t: ndarray + Temperature (K) + phase: str, optional + Defines the phase with respect to the saturation vapour pressure is computed. + It is either “water”, “ice” or “mixed”. + + Returns + ------- + ndarray + Saturation vapour pressure (Pa) + + + The algorithm was taken from the IFS model [IFS-CY47R3-PhysicalProcesses]_ (see Chapter 12). + It uses the following formula when ``phase`` is "water" or "ice": + + .. math:: + + e_{sat} = a_{1}\;exp \left(a_{3}\frac{t-273.16}{t-a_{4}}\right) + + where the parameters are set as follows: + + * ``phase`` = "water": :math:`a_{1}` =611.21 Pa, :math:`a_{3}` =17.502 and :math:`a_{4}` =32.19 K + * ``phase`` = "ice": :math:`a_{1}` =611.21 Pa, :math:`a_{3}` =22.587 and :math:`a_{4}` =-0.7 K + + When ``phase`` is "mixed" the formula is based on the value of ``t``: + + * if :math:`t <= t_{i}`: the formula for ``phase`` = "ice" is used (:math:`t_{i} = 250.16 K`) + * if :math:`t >= t_{0}`: the formula for ``phase`` = "water" is used (:math:`t_{0} = 273.16 K`) + * for the range :math:`t_{i} < t < t_{0}` an interpolation is used between the "ice" and "water" phases: + + .. math:: + + \alpha(t) e_{wsat}(t) + (1 - \alpha(t)) e_{isat}(t) + + with :math:`\alpha(t) = (\frac{t-t_{i}}{t_{0}-t_{i}})^2`. + + """ + return _EsComp().compute_es(t, phase) + + +def saturation_mixing_ratio(t, p, phase="mixed"): + r"""Computes the saturation mixing ratio from temperature with respect to a phase. + + Parameters + ---------- + t: ndarray + Temperature (K) + p: ndarray + Pressure (Pa) + phase: str + Defines the phase with respect to the :func:`saturation_vapour_pressure` is computed. + It is either “water”, “ice” or “mixed”. + + Returns + ------- + ndarray + Saturation mixing ratio (kg/kg) + + + Equivalent to the following code: + + .. code-block:: python + + e = saturation_vapour_pressure(t, phase=phase) + return mixing_ratio_from_vapour_pressure(e, p) + + """ + e = saturation_vapour_pressure(t, phase=phase) + return mixing_ratio_from_vapour_pressure(e, p) + + +def saturation_specific_humidity(t, p, phase="mixed"): + r"""Computes the saturation specific humidity from temperature with respect to a phase. + + Parameters + ---------- + t: ndarray + Temperature (K) + p: ndarray + Pressure (Pa) + phase: str, optional + Defines the phase with respect to the :func:`saturation_vapour_pressure` is computed. + It is either “water”, “ice” or “mixed”. + + Returns + ------- + ndarray + Saturation specific humidity (kg/kg) + + + Equivalent to the following code: + + .. code-block:: python + + e = saturation_vapour_pressure(t, phase=phase) + return specific_humidity_from_vapour_pressure(e, p) + + """ + e = saturation_vapour_pressure(t, phase=phase) + return specific_humidity_from_vapour_pressure(e, p) + + +def saturation_vapour_pressure_slope(t, phase="mixed"): + r"""Computes the slope of saturation vapour pressure with respect to temperature. + + Parameters + ---------- + t: ndarray + Temperature (K) + phase: str, optional + Defines the phase with respect to the computation will be performed. + It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` + for details. + + Returns + ------- + ndarray + Slope of saturation vapour pressure (Pa/K) + + """ + return _EsComp().compute_slope(t, phase) + + +def saturation_mixing_ratio_slope( + t, p, es=None, es_slope=None, phase="mixed", eps=1e-4 +): + r"""Computes the slope of saturation mixing ratio with respect to temperature. + + Parameters + ---------- + t: ndarray + Temperature (K) + p: ndarray + Pressure (Pa) + es: ndarray or None, optional + :func:`saturation_vapour_pressure` pre-computed for the given ``phase`` (Pa) + es_slope: ndarray or None, optional + :func:`saturation_vapour_pressure_slope` pre-computed for the given ``phase`` (Pa/K) + phase: str, optional + Defines the phase with respect to the computation will be performed. + It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` + for details. + eps: number + Where p - es < ``eps`` np.nan is returned. + + Returns + ------- + ndarray + Slope of saturation mixing ratio (:math:`kg kg^{-1} K^{-1}`) + + + The computation is based on the following formula: + + .. math:: + + \frac{\partial w_{s}}{\partial t} = \frac{\epsilon\; p}{(p-e_{s})^{2}} \frac{d e_{s}}{d t} + + where + + * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + * :math:`e_{s}` is the :func:`saturation_vapour_pressure` for the given ``phase`` + + """ + if eps <= 0: + raise ValueError(f"saturation_mixing_ratio_slope(): eps={eps} must be > 0") + if es is None: + es = saturation_vapour_pressure(t, phase=phase) + if es_slope is None: + es_slope = saturation_vapour_pressure_slope(t, phase=phase) + + v = np.asarray(p - es) + v[v < eps] = np.nan + return constants.epsilon * es_slope * p / np.square(v) + + +def saturation_specific_humidity_slope( + t, p, es=None, es_slope=None, phase="mixed", eps=1e-4 +): + r"""Computes the slope of saturation specific humidity with respect to temperature. + + Parameters + ---------- + t: ndarray + Temperature (K) + p: ndarray + Pressure (Pa) + es: ndarray or None, optional + :func:`saturation_vapour_pressure` pre-computed for the given ``phase`` (Pa) + es_slope: ndarray or None, optional + :func:`saturation_vapour_pressure_slope` pre-computed for the given ``phase`` (Pa/K) + phase: str, optional + Defines the phase with respect to the computation will be performed. + It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` + for details. + eps: number + Where p - es < ``eps`` np.nan is returned. + + Returns + ------- + ndarray + Slope of saturation specific humidity (:math:`kg kg^{-1} K^{-1}`) + + + The computation is based on the following formula: + + .. math:: + + \frac{\partial q_{s}}{\partial t} = + \frac{\epsilon\; p}{(p+e_{s}(\epsilon - 1))^{2}} \frac{d e_{s}}{d t} + + where + + * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + * :math:`e_{s}` is the :func:`saturation_vapour_pressure` for the given ``phase`` + + """ + if eps <= 0: + raise ValueError(f"saturation_specific_humidity_slope(): eps={eps} must be > 0") + if es is None: + es = saturation_vapour_pressure(t, phase=phase) + if es_slope is None: + es_slope = saturation_vapour_pressure_slope(t, phase=phase) + + v = np.asarray(np.square(p + es * (constants.epsilon - 1.0))) + v[np.asarray(p - es) < eps] = np.nan + return constants.epsilon * es_slope * p / v + + +def temperature_from_saturation_vapour_pressure(es): + r"""Computes the temperature from saturation vapour pressure. + + Parameters + ---------- + es: ndarray + :func:`saturation_vapour_pressure` (Pa) + + Returns + ------- + ndarray + Temperature (K) + + + The computation is always based on the "water" phase of + the :func:`saturation_vapour_pressure` formulation irrespective of the + phase ``es`` was computed to. + + """ + return _EsComp().t_from_es(es) + + +def relative_humidity_from_dewpoint(t, td): + r"""Computes the relative humidity from dewpoint temperature. + + Parameters + ---------- + t: ndarray + Temperature (K) + td: ndarray + Dewpoint (K) + + Returns + ------- + ndarray + Relative humidity (%) + + + The computation is based on the following formula: + + .. math:: + + r = 100 \frac {e_{wsat}(td)}{e_{wsat}(t)} + + where :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water. + + """ + e = saturation_vapour_pressure(td, phase="water") + es = saturation_vapour_pressure(t, phase="water") + return 100.0 * e / es + + +def relative_humidity_from_specific_humidity(t, q, p): + r"""Computes the relative humidity from specific humidity. + + Parameters + ---------- + t: ndarray + Temperature (K) + q: ndarray + Specific humidity (kg/kg) + p: ndarray + Pressure (Pa) + + Returns + ------- + ndarray + Relative humidity (%) + + + The computation is based on the following formula: + + .. math:: + + r = 100 \frac {e(q, p)}{e_{msat}(t)} + + where: + + * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) + * :math:`e_{msat}` is the :func:`saturation_vapour_pressure` based on the "mixed" phase + + """ + svp = saturation_vapour_pressure(t) + e = vapour_pressure_from_specific_humidity(q, p) + return 100.0 * e / svp + + +def specific_humidity_from_dewpoint(td, p): + r"""Computes the specific humidity from dewpoint. + + Parameters + ---------- + td: ndarray + Dewpoint (K) + p: ndarray + Pressure (Pa) + + Returns + ------- + ndarray + Specific humidity (kg/kg) + + + The computation starts with determining the vapour pressure: + + .. math:: + + e(q, p) = e_{wsat}(td) + + where: + + * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) + * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water + * :math:`q` is the specific humidity + + Then `q` is computed from :math:`e` using :func:`specific_humidity_from_vapour_pressure`. + + """ + svp = saturation_vapour_pressure(td, phase="water") + return specific_humidity_from_vapour_pressure(svp, p) + + +def mixing_ratio_from_dewpoint(td, p): + r"""Computes the mixing ratio from dewpoint. + + Parameters + ---------- + td: ndarray + Dewpoint (K) + p: ndarray + Pressure (Pa) + + Returns + ------- + ndarray + Specific humidity (kg/kg) + + + The computation starts with determining the vapour pressure: + + .. math:: + + e(w, p) = e_{wsat}(td) + + where: + + * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_mixing_ratio`) + * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water + * :math:`w` is the mixing ratio + + Then `w` is computed from :math:`e` using :func:`mixing_ratio_from_vapour_pressure`. + + """ + svp = saturation_vapour_pressure(td, phase="water") + return mixing_ratio_from_vapour_pressure(svp, p) + + +def specific_humidity_from_relative_humidity(t, r, p): + r"""Computes the specific humidity from relative_humidity. + + Parameters + ---------- + t: ndarray + Temperature (K) + r: ndarray + Relative humidity(%) + p: ndarray + Pressure (Pa) + + Returns + ------- + ndarray + Specific humidity (kg/kg) units + + + The computation starts with determining the the vapour pressure: + + .. math:: + + e(q, p) = r\; \frac{e_{msat}(t)}{100} + + where: + + * :math:`e` is the vapour pressure (see :func:`vapour_pressure`) + * :math:`e_{msat}` is the :func:`saturation_vapour_pressure` based on the "mixed" phase + * :math:`q` is the specific humidity + + Then :math:`q` is computed from :math:`e` using :func:`specific_humidity_from_vapour_pressure`. + + """ + e = r * saturation_vapour_pressure(t) / 100.0 + return specific_humidity_from_vapour_pressure(e, p) + + +def dewpoint_from_relative_humidity(t, r): + r"""Computes the dewpoint temperature from relative humidity. + + Parameters + ---------- + t: ndarray + Temperature (K) + r: ndarray + Relative humidity (%) + + Returns + ------- + ndarray + Dewpoint (K) + + + The computation starts with determining the the saturation vapour pressure over + water at the dewpoint temperature: + + .. math:: + + e_{wsat}(td) = \frac{r\; e_{wsat}(t)}{100} + + where: + + * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water + * :math:`td` is the dewpoint. + + Then :math:`td` is computed from :math:`e_{wsat}(td)` by inverting the + equations used in :func:`saturation_vapour_pressure`. + + """ + es = saturation_vapour_pressure(t, phase="water") * r / 100.0 + return temperature_from_saturation_vapour_pressure(es) + + +def dewpoint_from_specific_humidity(q, p): + r"""Computes the dewpoint temperature from specific humidity. + + Parameters + ---------- + q: ndarray + Specific humidity (kg/kg) + p: ndarray + Pressure (Pa) + + Returns + ------- + ndarray + Dewpoint (K) + + + The computation starts with determining the the saturation vapour pressure over + water at the dewpoint temperature: + + .. math:: + + e_{wsat}(td) = e(q, p) + + where: + + * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) + * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water + * :math:`td` is the dewpoint + + Then :math:`td` is computed from :math:`e_{wsat}(td)` by inverting the equations + used in :func:`saturation_vapour_pressure`. + + """ + return temperature_from_saturation_vapour_pressure( + vapour_pressure_from_specific_humidity(q, p) + ) + + +def virtual_temperature(t, q): + r"""Computes the virtual temperature from temperature and specific humidity. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + q: number or ndarray + Specific humidity (kg/kg) + + Returns + ------- + number or ndarray + Virtual temperature (K) + + + The computation is based on the following formula [Wallace2006]_: + + .. math:: + + t_{v} = t (1 + \frac{1 - \epsilon}{\epsilon} q) + + with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + c1 = (1.0 - constants.epsilon) / constants.epsilon + return t * (1.0 + c1 * q) + + +def virtual_potential_temperature(t, q, p): + r"""Computes the virtual potential temperature from temperature and specific humidity. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + q: number or ndarray + Specific humidity (kg/kg) + p: number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Virtual potential temperature (K) + + + The computation is based on the following formula: + + .. math:: + + \Theta_{v} = \theta (1 + \frac{1 - \epsilon}{\epsilon} q) + + where: + + * :math:`\Theta` is the :func:`potential_temperature` + * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). + + """ + c1 = (1.0 - constants.epsilon) / constants.epsilon + return potential_temperature(t, p) * (1.0 + c1 * q) + + +def potential_temperature(t, p): + r"""Computes the potential temperature. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + p: number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Potential temperature (K) + + + The computation is based on the following formula [Wallace2006]_: + + .. math:: + + \theta = t (\frac{10^{5}}{p})^{\kappa} + + with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). + + """ + return t * np.power(constants.p0 / p, constants.kappa) + + +def temperature_from_potential_temperature(th, p): + r"""Computes the temperature from potential temperature. + + Parameters + ---------- + th: number or ndarray + Potential temperature (K) + p: number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Temperature (K) + + + The computation is based on the following formula: + + .. math:: + + t = \theta (\frac{p}{10^{5}})^{\kappa} + + with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). + + """ + return th * np.power(p / constants.p0, constants.kappa) + + +def pressure_on_dry_adiabat(t, t_def, p_def): + r"""Computes the pressure on a dry adiabat. + + Parameters + ---------- + t: number or ndarray + Temperature on the dry adiabat (K) + t_def: number or ndarray + Temperature defining the dry adiabat (K) + p_def: number or ndarray + Pressure defining the dry adiabat (Pa) + + Returns + ------- + number or ndarray + Pressure on the dry adiabat (Pa) + + + The computation is based on the following formula: + + .. math:: + + p = p_{def} (\frac{t}{t_{def}})^{\frac{1}{\kappa}} + + with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). + + """ + return p_def * np.power(t / t_def, 1 / constants.kappa) + + +def temperature_on_dry_adiabat(p, t_def, p_def): + r"""Computes the temperature on a dry adiabat. + + Parameters + ---------- + p: number or ndarray + Pressure on the dry adiabat (Pa) + t_def: number or ndarray + Temperature defining the dry adiabat (K) + p_def: number or ndarray + Pressure defining the dry adiabat (Pa) + + Returns + ------- + number or ndarray + Temperature on the dry adiabat (K) + + + The computation is based on the following formula: + + .. math:: + + t = t_{def} (\frac{p}{p_{def}})^{\kappa} + + with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). + + """ + return t_def * np.power(p / p_def, constants.kappa) + + +def lcl_temperature(t, td, method="davies"): + r"""Computes the Lifting Condenstaion Level (LCL) temperature from dewpoint. + + Parameters + ---------- + t: number or ndarray + Temperature at the start level (K) + td: number or ndarray + Dewpoint at the start level (K) + method: str, optional + The computation method: "davies" or "bolton". + + Returns + ------- + number or ndarray + Temperature of the LCL (K) + + + The actual computation is based on the ``method``: + + * "davies": the formula by [DaviesJones1983]_ is used (it is also used by the IFS model): + + .. math:: + + t_{LCL} = + td - (0.212 + 1.571\times 10^{-3} (td - t_{0}) - 4.36\times 10^{-4} (t - t_{0})) (t - td) + + where :math:`t_{0}` is the triple point of water (see :data:`earthkit.meteo.constants.T0`). + + * "bolton": the formula by [Bolton1980]_ is used: + + .. math:: + + t_{LCL} = 56.0 + \frac{1}{\frac{1}{td - 56} + \frac{log(\frac{t}{td})}{800}} + + """ + # Davies-Jones formula + if method == "davies": + t_lcl = td - ( + 0.212 + 1.571e-3 * (td - constants.T0) - 4.36e-4 * (t - constants.T0) + ) * (t - td) + return t_lcl + # Bolton formula + elif method == "bolton": + return 56.0 + 1 / (1 / (td - 56) + np.log(t / td) / 800) + else: + raise ValueError(f"lcl_temperature: invalid method={method} specified!") + + +def lcl(t, td, p, method="davies"): + r"""Computes the temperature and pressure of the Lifting Condenstaion Level (LCL) from dewpoint. + + Parameters + ---------- + t: number or ndarray + Temperature at the start level (K) + td: number or ndarray + Dewpoint at the start level (K) + p: number or ndarray + Pressure at the start level (Pa) + method: str + method: str, optional + The computation method: "davies" or "bolton". + + Returns + ------- + number or ndarray + Temperature of the LCL (K) + number or ndarray + Pressure of the LCL (Pa) + + + The LCL temperature is determined by :func:`lcl_temperature` with the given ``method`` + and the pressure is computed with :math:`t_{LCL}` using :func:`pressure_on_dry_adiabat`. + + """ + t_lcl = lcl_temperature(t, td, method=method) + p_lcl = pressure_on_dry_adiabat(t_lcl, t, p) + return t_lcl, p_lcl + + +class _ThermoState: + def __init__(self, t=None, td=None, q=None, w=None, p=None): + self.t = t + self.td = td + self.q = q + self.w = w + self.p = p + self.e = None + self.es = None + self.ws = None + self.qs = None + + +class _EptComp: + CM = {} + c_lambda = 1.0 / constants.kappa + + @staticmethod + def make(method): + return _EptComp.CM[method]() + + def is_mixing_ratio_based(self): + return True + + def compute_ept(self, t=None, td=None, q=None, p=None): + if td is None and q is None: + raise ValueError("ept: either td or q must have a valid value!") + + # compute dewpoint since all the methods require it + if td is None: + td = dewpoint_from_specific_humidity(q, p) + + ths = _ThermoState(t=t, td=td, q=q, p=p) + return self._ept(ths) + + def compute_ept_sat(self, t, p): + ths = _ThermoState(t=t, p=p) + return self._th_sat(ths) * np.exp(self._G_sat(ths)) + + def compute_wbpt(self, ept): + t0 = 273.16 + x = ept / t0 + a = [7.101574, -20.68208, 16.11182, 2.574631, -5.205688] + b = [1.0, -3.552497, 3.781782, -0.6899655, -0.5929340] + return ept - np.exp( + np.polynomial.polynomial.polyval(x, a) + / np.polynomial.polynomial.polyval(x, b) + ) + + def compute_t_on_ma_stipanuk(self, ept, p): + if isinstance(p, np.ndarray): + t = np.full(p.shape, constants.T0 - 20) + elif isinstance(ept, np.ndarray): + t = np.full(ept.shape, constants.T0 - 20) + else: + t = constants.T0 - 20 + max_iter = 12 + dt = 120.0 + for _ in range(max_iter): + ths = _ThermoState(t=t, p=p) + dt /= 2.0 + t += ( + np.sign(ept * np.exp(self._G_sat(ths, scale=-1.0)) - self._th_sat(ths)) + * dt + ) + # ths.t = t + # return ept - self._th_sat(ths) * np.exp(self._G_sat(ths)) + return t + + def compute_t_on_ma_davies(self, ept, p): + def _k1(pp): + """Function k1 in the article.""" + a = [-53.737, 137.81, -38.5] + return np.polynomial.polynomial.polyval(pp, a) + + def _k2(pp): + """Function k2 in the article.""" + a = [-0.384, 56.831, -4.392] + return np.polynomial.polynomial.polyval(pp, a) + + def _D(p): + """Function D in the article.""" + return 1.0 / (0.1859e-5 * p + 0.6512) + + max_iter = 1 + A = 2675 + t0 = 273.16 + + if not isinstance(p, np.ndarray): + p = np.full(ept.shape, p) + + tw = ept.copy() + pp = np.power(p / constants.p0, constants.kappa) + te = ept * pp + c_te = np.power(t0 / te, _EptComp.c_lambda) + + # initial guess + mask = c_te > _D(p) + if np.any(mask): + es = saturation_vapour_pressure(te[mask]) + ws = mixing_ratio_from_vapour_pressure(es, p[mask]) + d_es = saturation_vapour_pressure_slope(te[mask]) + tw[mask] = te[mask] - t0 - (A * ws) / (1 + A * ws * d_es / es) + + mask = (1 <= c_te) & (c_te <= _D(p)) + tw[mask] = _k1(pp[mask]) - _k2(pp[mask]) * c_te[mask] + + mask = (0.4 <= c_te) & (c_te < 1) + tw[mask] = (_k1(pp[mask]) - 1.21) - (_k2(pp[mask]) - 1.21) * c_te[mask] + + mask = c_te < 0.4 + tw[mask] = ( + (_k1(pp[mask]) - 2.66) + - (_k2(pp[mask]) - 1.21) * c_te[mask] + + 0.58 / c_te[mask] + ) + # tw has to be converted to Kelvin + tw = celsius_to_kelvin(tw) + + # np.place(tw, tw > 362.16, 362.1) + # print(f"ept={ept}") + # print(f"tw={tw}") + + for i in range(max_iter): + ths = _ThermoState(t=tw, p=p) + ths.c_tw = np.power(t0 / ths.t, _EptComp.c_lambda) + ths.es = saturation_vapour_pressure(ths.t) + if self.is_mixing_ratio_based(): + ths.ws = mixing_ratio_from_vapour_pressure(ths.es, ths.p) + else: + ths.qs = specific_humidity_from_vapour_pressure(ths.es, ths.p) + f_val = self._f(ths) + # print(f"{i}") + # print(f" p={p}") + # print(f" t={ths.t}") + # print(f" f_val={f_val}") + # print(f" f_val_cte={f_val - c_te}") + # print(f" d_lnf={self._d_lnf(ths)}") + # print(tw[mask].shape) + # print(f_val.shape) + tw -= (f_val - c_te) / (f_val * self._d_lnf(ths)) + # print(f"tw={tw}") + + # when ept is derived with the Bolton methods the iteration breaks down for extremeley + # hot and humid conditions (roughly about t > 50C and r > 95%) and results in negative + # tw values! + np.place(tw, tw <= 0, np.nan) + + # ths = _ThermoState(t=tw, p=p) + # return ept - self._th_sat(ths) * np.exp(self._G_sat(ths)) + return tw + + +class _EptCompIfs(_EptComp): + def __init__(self): + self.K0 = constants.Lv / constants.c_pd + + def is_mixing_ratio_based(self): + return False + + def _ept(self, ths): + th = potential_temperature(ths.t, ths.p) + t_lcl = lcl_temperature(ths.t, ths.td, method="davies") + if ths.q is None: + ths.q = specific_humidity_from_dewpoint(ths.td, ths.p) + return th * np.exp(self.K0 * ths.q / t_lcl) + + def _th_sat(self, ths): + return potential_temperature(ths.t, ths.p) + + def _G_sat(self, ths, scale=1.0): + qs = saturation_specific_humidity(ths.t, ths.p) + return (scale * self.K0) * qs / ths.t + + def _d_G_sat(self, ths): + if ths.qs is None: + ths.qs = saturation_specific_humidity(ths.t, ths.p) + return ( + -self.K0 * ths.qs / (ths.t**2) + + self.K0 * saturation_specific_humidity_slope(ths.t, ths.p) / ths.t + ) + + def _f(self, ths): + return ths.c_tw * np.exp(self._G_sat(ths, scale=-self.c_lambda)) + + def _d_lnf(self, ths): + return -self.c_lambda * (1 / ths.t + self._d_G_sat(ths)) + + +class _EptCompBolton35(_EptComp): + def __init__(self): + self.K0 = 2675.0 + self.K3 = 0.28 + + def _ept(self, ths): + t_lcl = lcl_temperature(ths.t, ths.td, method="bolton") + if ths.q is None: + w = mixing_ratio_from_dewpoint(ths.td, ths.p) + else: + w = mixing_ratio_from_specific_humidity(ths.q) + th = ths.t * np.power(constants.p0 / ths.p, constants.kappa * (1 - self.K3 * w)) + return th * np.exp(self.K0 * w / t_lcl) + + def _th_sat(self, ths): + if ths.ws is None: + ths.ws = saturation_mixing_ratio(ths.t, ths.p) + return ths.t * np.power( + constants.p0 / ths.p, constants.kappa * (1 - self.K3 * ths.ws) + ) + + def _G_sat(self, ths, scale=1.0): + if ths.ws is None: + ths.ws = saturation_mixing_ratio(ths.t, ths.p) + return (scale * self.K0) * ths.ws / ths.t + + def _d_G_sat(self, ths): + return ( + -self.K0 * ths.ws / np.square(ths.t) + + self.K0 * saturation_mixing_ratio_slope(ths.t, ths.p) / ths.t + ) + + def _f(self, ths): + # print(f" c_tw={ths.c_tw}") + # print(f" es_frac={ths.es / ths.p}") + # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}") + return ( + ths.c_tw + * np.power(ths.p / constants.p0, self.K3 * ths.ws) + * np.exp(self._G_sat(ths, scale=-self.c_lambda)) + ) + + def _d_lnf(self, ths): + return -self.c_lambda * ( + 1 / ths.t + + self.K3 + * np.log(ths.p / constants.p0) + * saturation_vapour_pressure_slope(ths.t) + + self._d_G_sat(ths) + ) + + +class _EptCompBolton39(_EptComp): + def __init__(self): + # Comment on the Bolton formulas. The constants used by Bolton to + # derive the formulas differ from the ones used by earthkit.meteo E.g. + # Bolton earthkit.meteo + # Rd 287.04 287.0597 + # cpd 1005.7 1004.79 + # kappa 0.2854 0.285691 + # epsilon 0.6220 0.621981 + + self.K0 = 3036.0 + self.K1 = 1.78 + self.K2 = 0.448 + self.K4 = 0.28 + + def _ept(self, ths): + t_lcl = lcl_temperature(ths.t, ths.td, method="bolton") + if ths.q is None: + w = mixing_ratio_from_dewpoint(ths.td, ths.p) + else: + w = mixing_ratio_from_specific_humidity(ths.q) + + e = vapour_pressure_from_mixing_ratio(w, ths.p) + th = potential_temperature(ths.t, ths.p - e) * np.power( + ths.t / t_lcl, self.K4 * w + ) + return th * np.exp((self.K0 / t_lcl - self.K1) * w * (1.0 + self.K2 * w)) + + def _th_sat(self, ths): + if ths.es is None: + ths.es = saturation_vapour_pressure(ths.t) + ths.es[ths.p - ths.es < 1e-4] = np.nan + return potential_temperature(ths.t, ths.p - ths.es) + + def _G_sat(self, ths, scale=1.0): + if ths.es is None: + ths.es = saturation_vapour_pressure(ths.t) + ths.es[ths.p - ths.es < 1e-4] = np.nan + # print(f" es={ths.es}") + ws = mixing_ratio_from_vapour_pressure(ths.es, ths.p) + # print(f" ws={ws}") + return ( + ((scale * self.K0) / ths.t - (scale * self.K1)) * ws * (1.0 + self.K2 * ws) + ) + + def _d_G_sat(self, ths): + # print(f" d_ws={saturation_mixing_ratio_slope(ths.t, ths.p)}") + return -self.K0 * (ths.ws + self.K2 * np.square(ths.ws)) / ( + np.square(ths.t) + ) + (self.K0 / ths.t - self.K1) * ( + 1 + (2 * self.K2) * ths.ws + ) * saturation_mixing_ratio_slope( + ths.t, ths.p + ) + + def _f(self, ths): + # print(f" c_tw={ths.c_tw}") + # print(f" es_frac={ths.es / ths.p}") + # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}") + return ( + ths.c_tw + * (1 - ths.es / ths.p) + * np.exp(self._G_sat(ths, scale=-self.c_lambda)) + ) + + def _d_lnf(self, ths): + return -self.c_lambda * ( + 1 / ths.t + + constants.kappa + * saturation_vapour_pressure_slope(ths.t) + / (ths.p - ths.es) + + self._d_G_sat(ths) + ) + + +_EptComp.CM = { + "ifs": _EptCompIfs, + "bolton35": _EptCompBolton35, + "bolton39": _EptCompBolton39, +} + + +def ept_from_dewpoint(t, td, p, method="ifs"): + r"""Computes the equivalent potential temperature from dewpoint. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + td: number or ndarray + Dewpoint (K) + p: number or ndarray + Pressure (Pa) + method: str, optional + Specifies the computation method. The possible values are: "ifs", "bolton35", "bolton39". + + Returns + ------- + number or ndarray + Equivalent potential temperature (K) + + + The actual computation is based on the value of ``method``: + + * "ifs": the formula from the IFS model [IFS-CY47R3-PhysicalProcesses]_ (Chapter 6.11) is used: + + .. math:: + + \Theta_{e} = \Theta\; exp(\frac{L_{v}\; q}{c_{pd}\; t_{LCL}}) + + * "bolton35": Eq (35) from [Bolton1980]_ is used: + + + .. math:: + + \Theta_{e} = \Theta (\frac{10^{5}}{p})^{\kappa 0.28 w} exp(\frac{2675 w}{t_{LCL}}) + + * "bolton39": Eq (39) from [Bolton1980]_ is used: + + .. math:: + + \Theta_{e} = + t (\frac{10^{5}}{p-e})^{\kappa} (\frac{t}{t_{LCL}})^{0.28 w} exp[(\frac{3036}{t_{LCL}} - + 1.78)w(1+0.448\; w)] + + where: + + * :math:`\Theta` is the :func:`potential_temperature` + * :math:`t_{LCL}` is the temperature at the Lifting Condestation Level computed + with :func:`lcl_temperature` using option: + + * method="davis" when ``method`` is "ifs" + * method="bolton" when ``method`` is "bolton35" or "bolton39" + * :math:`q` is the specific humidity computed with :func:`specific_humidity_from_dewpoint` + * :math:`w`: is the mixing ratio computed with :func:`mixing_ratio_from_dewpoint` + * :math:`e` is the vapour pressure computed with :func:`vapour_pressure_from_mixing_ratio` + * :math:`L_{v}`: is the latent heat of vaporisation + (see :data:`earthkit.meteo.constants.Lv`) + * :math:`c_{pd}` is the specific heat of dry air on constant pressure + (see :data:`earthkit.meteo.constants.c_pd`) + * :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`) + + """ + return _EptComp.make(method).compute_ept(t=t, td=td, p=p) + + +def ept_from_specific_humidity(t, q, p, method="ifs"): + r"""Computes the equivalent potential temperature from specific humidity. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + q: number or ndarray + Specific humidity (kg/kg) + p: number or ndarray + Pressure (Pa) + method: str, optional + Specifies the computation method. The possible values are: "ifs", + "bolton35", "bolton39. See :func:`ept_from_dewpoint` for details. + + Returns + ------- + number or ndarray + Equivalent potential temperature (K) + + + The computations are the same as in :func:`ept_from_dewpoint` + (the dewpoint is computed from q with :func:`dewpoint_from_specific_humidity`). + + """ + return _EptComp.make(method).compute_ept(t=t, q=q, p=p) + + +def saturation_ept(t, p, method="ifs"): + r"""Computes the saturation equivalent potential temperature. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + p: number or ndarray + Pressure (Pa) + method: str, optional + Specifies the computation method. The possible values are: "ifs", "bolton35", "bolton39". + + Returns + ------- + number or ndarray + Saturation equivalent potential temperature (K) + + + The actual computation is based on the ``method``: + + * "ifs": The formula is based on the equivalent potential temperature definition used + in the IFS model [IFS-CY47R3-PhysicalProcesses]_ (see Chapter 6.11) : + + .. math:: + + \Theta_{esat} = \Theta\; exp(\frac{L_{v}\; q_{sat}}{c_{pd}\; t}) + + * "bolton35": Eq (35) from [Bolton1980]_ is used: + + .. math:: + + \Theta_{e} = \Theta (\frac{10^{5}}{p})^{\kappa 0.28 w_{sat}}\; exp(\frac{2675\; w_{sat}}{t}) + + * "bolton39": Eq (39) from [Bolton1980]_ is used: + + .. math:: + + \Theta_{e} = + t (\frac{10^{5}}{p-e_{sat}})^{\kappa} exp[(\frac{3036}{t} - 1.78)w_{sat}(1+0.448\; w_{sat})] + + where: + + * :math:`\Theta` is the :func:`potential_temperature` + * :math:`e_{sat}` is the :func:`saturation_vapor_pressure` + * :math:`q_{sat}` is the :func:`saturation_specific_humidity` + * :math:`w_{sat}` is the :func:`saturation_mixing_ratio` + * :math:`L_{v}` is the specific latent heat of vaporization (see :data:`earthkit.meteo.constants.Lv`) + * :math:`c_{pd}` is the specific heat of dry air on constant pressure + (see :data:`earthkit.meteo.constants.c_pd`) + + """ + return _EptComp.make(method).compute_ept_sat(t, p) + + +def temperature_on_moist_adiabat(ept, p, ept_method="ifs", t_method="bisect"): + r"""Computes the temperature on a moist adiabat (pseudoadiabat) + + Parameters + ---------- + ept: number or ndarray + Equivalent potential temperature defining the moist adiabat (K) + p: number or ndarray + Pressure on the moist adiabat (Pa) + ept_method: str, optional + Specifies the computation method that was used to compute ``ept``. The possible + values are: "ifs", "bolton35", "bolton39". + (See :func:`ept_from_dewpoint` for details.) + t_method: str, optional + Specifies the iteration method along the moist adiabat to find the temperature + for the given ``p`` pressure. The possible values are as follows: + + * "bisect": a bisection method is used as defined in [Stipanuk1973]_ + * "newton": Newtons's method is used as defined by Eq (2.6) in [DaviesJones2008]_. + For extremely hot and humid conditions (``ept`` > 800 K) depending on + ``ept_method`` the computation might not be carried out + and np.nan will be returned. + + + Returns + ------- + number or ndarray + Temperature on the moist adiabat (K). For values where the computation cannot + be carried out np.nan is returned. + + """ + cm = _EptComp.make(ept_method) + if t_method == "bisect": + return cm.compute_t_on_ma_stipanuk(ept, p) + elif t_method == "newton": + return cm.compute_t_on_ma_davies(ept, p) + else: + raise ValueError( + f"temperature_on_moist_adiabat: invalid t_method={t_method} specified!" + ) + + +def wet_bulb_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="bisect"): + r"""Computes the pseudo adiabatic wet bulb temperature from dewpoint. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + td: number or ndarray + Dewpoint (K) + p: number or ndarray + Pressure (Pa) + ept_method: str, optional + Specifies the computation method for the equivalent potential temperature. + The possible values are: "ifs", "bolton35", "bolton39". + (See :func:`ept_from_dewpoint` for details.) + t_method: str, optional + Specifies the method to find the temperature along the moist adiabat defined + by the equivalent potential temperature. The possible values are as follows: + + * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used + * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used + + Returns + ------- + number or ndarray + Wet bulb temperature (K) + + + The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): + + * first the equivalent potential temperature is computed with the given + ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. + * then the wet bulb potential temperature is determined as the temperature at + pressure ``p`` on the moist adiabat with the given ``t_method``. + + """ + ept = ept_from_dewpoint(t, td, p, method=ept_method) + return temperature_on_moist_adiabat( + ept, p, ept_method=ept_method, t_method=t_method + ) + + +def wet_bulb_temperature_from_specific_humidity( + t, q, p, ept_method="ifs", t_method="bisect" +): + r"""Computes the pseudo adiabatic wet bulb temperature from specific humidity. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + q: number or ndarray + Specific humidity (kg/kg) + p: number or ndarray + Pressure (Pa) + ept_method: str, optional + Specifies the computation method for the equivalent potential temperature. + The possible values are: "ifs", "bolton35", "bolton39". + (See :func:`ept_from_dewpoint` for details.) + t_method: str, optional + Specifies the method to find the temperature along the moist adiabat + defined by the equivalent potential temperature. The possible values are + as follows: + + * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used + * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used + + Returns + ------- + number or ndarray + Wet bulb temperature (K) + + + The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): + + * first the equivalent potential temperature is computed with the given + ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. + * then the wet bulb potential temperature is determined as the temperature at + pressure ``p`` on the moist adiabat with the given ``t_method``. + + """ + ept = ept_from_specific_humidity(t, q, p, method=ept_method) + return temperature_on_moist_adiabat( + ept, p, ept_method=ept_method, t_method=t_method + ) + + +def wet_bulb_potential_temperature_from_dewpoint( + t, td, p, ept_method="ifs", t_method="direct" +): + r"""Computes the pseudo adiabatic wet bulb potential temperature from dewpoint. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + td: number or ndarray + Dewpoint (K) + p: number or ndarray + Pressure (Pa) + ept_method: str, optional + Specifies the computation method for the equivalent potential temperature. + The possible values are: "ifs", "bolton35", "bolton39". + (See :func:`ept_from_dewpoint` for details.) + t_method: str, optional + Specifies the method to find the temperature along the moist adiabat defined + by the equivalent potential temperature. The possible values are as follows: + + * "direct": the rational formula defined by Eq (3.8) in [DaviesJones2008]_ is used + * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used + * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used + + Returns + ------- + number or ndarray + Wet bulb potential temperature (K) + + + The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): + + * first the equivalent potential temperature is computed with the given + ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. + * then the wet bulb potential temperature is determined as the temperature at + pressure :math:`10^{5}` Pa on the moist adiabat with the given ``t_method``. + + """ + ept = ept_from_dewpoint(t, td, p, method=ept_method) + if t_method == "direct": + return _EptComp.make(ept_method).compute_wbpt(ept) + else: + return temperature_on_moist_adiabat( + ept, constants.p0, ept_method=ept_method, t_method=t_method + ) + + +def wet_bulb_potential_temperature_from_specific_humidity( + t, q, p, ept_method="ifs", t_method="direct" +): + r"""Computes the pseudo adiabatic wet bulb potential temperature from specific humidity. + + Parameters + ---------- + t: number or ndarray + Temperature (K) + q: number or ndarray + Specific humidity (kg/kg) + p: number or ndarray + Pressure (Pa) + ept_method: str, optional + Specifies the computation method for the equivalent potential temperature. + The possible values are: "ifs", "bolton35", "bolton39". + (See :func:`ept_from_dewpoint` for details.) + t_method: str, optional + Specifies the method to find the temperature along the moist adiabat + defined by the equivalent potential temperature. The possible values are as follows: + + * "direct": the rational formula defined by Eq (3.8) in [DaviesJones2008]_ is used + * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used + * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used + + Returns + ------- + number or ndarray + Wet bulb potential temperature (K) + + + The computations are the same as in + :func:`wet_bulb_potential_temperature_from_dewpoint` + (the dewpoint is computed from q with :func:`dewpoint_from_specific_humidity`). + + """ + ept = ept_from_specific_humidity(t, q, p, method=ept_method) + if t_method == "direct": + return _EptComp.make(ept_method).compute_wbpt(ept) + else: + return temperature_on_moist_adiabat( + ept, constants.p0, ept_method=ept_method, t_method=t_method + ) diff --git a/earthkit/meteo/thermo/thermo.py b/earthkit/meteo/thermo/thermo.py index 3270013..b6df320 100644 --- a/earthkit/meteo/thermo/thermo.py +++ b/earthkit/meteo/thermo/thermo.py @@ -1,1790 +1,162 @@ -# (C) Copyright 2021- ECMWF. - +# (C) Copyright 2021 ECMWF. +# # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. - - -import numpy as np - -from earthkit.meteo import constants - - -def celsius_to_kelvin(t): - """Converts temperature values from Celsius to Kelvin. - - Parameters - ---------- - t : number or ndarray - Temperature in Celsius units - - Returns - ------- - number or ndarray - Temperature in Kelvin units - - """ - return t + constants.T0 - - -def kelvin_to_celsius(t): - """Converts temperature values from Kelvin to Celsius. - - Parameters - ---------- - t : number or ndarray - Temperature in Kelvin units - - Returns - ------- - number or ndarray - Temperature in Celsius units - - """ - return t - constants.T0 - - -def specific_humidity_from_mixing_ratio(w): - r"""Computes the specific humidity from mixing ratio. - - Parameters - ---------- - w : number or ndarray - Mixing ratio (kg/kg) - - Returns - ------- - number or ndarray - Specific humidity (kg/kg) - - - The result is the specific humidity in kg/kg units. The computation is based on - the following definition [Wallace2006]_: - - .. math:: - - q = \frac {w}{1+w} - - """ - return w / (1 + w) - - -def mixing_ratio_from_specific_humidity(q): - r"""Computes the mixing ratio from specific humidity. - - Parameters - ---------- - q : number or ndarray - Specific humidity (kg/kg) - - Returns - ------- - number or ndarray - Mixing ratio (kg/kg) - - - The result is the mixing ratio in kg/kg units. The computation is based on - the following definition [Wallace2006]_: - - .. math:: - - w = \frac {q}{1-q} - - """ - return q / (1 - q) - - -def vapour_pressure_from_specific_humidity(q, p): - r"""Computes the vapour pressure from specific humidity. - - Parameters - ---------- - q: number or ndarray - Specific humidity (kg/kg) - p: number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Vapour pressure (Pa) - - - The computation is based on the following formula [Wallace2006]_: - - .. math:: - - e = \frac{p\;q}{\epsilon\; (1 + q(\frac{1}{\epsilon} -1 ))} - - with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - c = constants.epsilon * (1 / constants.epsilon - 1) - return (p * q) / (constants.epsilon + c * q) - - -def vapour_pressure_from_mixing_ratio(w, p): - r"""Computes the vapour pressure from mixing ratio. - - Parameters - ---------- - w: number or ndarray - Mixing ratio (kg/kg) - p: number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Vapour pressure (Pa) - - - The computation is based on the following formula: - - .. math:: - - e = \frac{p\;w}{\epsilon + w} - - with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - return (p * w) / (constants.epsilon + w) - - -def specific_humidity_from_vapour_pressure(e, p, eps=1e-4): - r"""Computes the specific humidity from vapour pressure. - - Parameters - ---------- - e: number or ndarray - Vapour pressure (Pa) - p: number or ndarray - Pressure (Pa) - eps: number - Where p - e < ``eps`` np.nan is returned. - - Returns - ------- - number or ndarray - Specific humidity (kg/kg) - - - The computation is based on the following formula: - - .. math:: - - q = \frac{\epsilon e}{p + e(\epsilon-1)} - - with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - if eps <= 0: - raise ValueError( - f"specific_humidity_from_vapour_pressure(): eps={eps} must be > 0" - ) - - v = np.asarray(p + (constants.epsilon - 1) * e) - v[np.asarray(p - e) < eps] = np.nan - return constants.epsilon * e / v - - -def mixing_ratio_from_vapour_pressure(e, p, eps=1e-4): - r"""Computes the mixing ratio from vapour pressure. - - Parameters - ---------- - e: number or ndarray - Vapour pressure (Pa) - p: number or ndarray - Pressure (Pa) - eps: number - Where p - e < ``eps`` np.nan is returned. - - Returns - ------- - number or ndarray - Mixing ratio (kg/kg). - - - The computation is based on the following formula: - - .. math:: - - w = \frac{\epsilon e}{p - e} - - with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - if eps <= 0: - raise ValueError(f"mixing_ratio_from_vapour_pressure(): eps={eps} must be > 0") - - v = np.asarray(p - e) - v[v < eps] = np.nan - return constants.epsilon * e / v - - -class _EsComp: - c1 = 611.21 - c3w = 17.502 - c4w = 32.19 - c3i = 22.587 - c4i = -0.7 - t0 = 273.16 - ti = t0 - 23 - PHASES = ["mixed", "water", "ice"] - - def _check_phase(self, phase): - if phase not in _EsComp.PHASES: - raise ValueError( - f"saturation_vapour_pressure(): invalid phase={phase}! Allowed values = {_EsComp.PHASES}" - ) - return True - - def compute_es(self, t, phase): - self._check_phase(phase) - if phase == "mixed": - return self._es_mixed(t) - elif phase == "water": - return self._es_water(t) - elif phase == "ice": - return self._es_ice(t) - - def compute_slope(self, t, phase): - self._check_phase(phase) - if phase == "mixed": - return self._es_mixed_slope(t) - elif phase == "water": - return self._es_water_slope(t) - elif phase == "ice": - return self._es_ice_slope(t) - - def t_from_es(self, es): - v = np.log(es / self.c1) - return (v * self.c4w - self.c3w * self.t0) / (v - self.c3w) - - def _es_water(self, t): - return self.c1 * np.exp(self.c3w * (t - self.t0) / (t - self.c4w)) - - def _es_ice(self, t): - return self.c1 * np.exp(self.c3i * (t - self.t0) / (t - self.c4i)) - - def _es_mixed(self, t): - # Fraction of liquid water (=alpha): - # t <= ti => alpha=0 - # t > ti and t < t0 => alpha=(t-ti)/(t0-ti))^2 - # t >= t0 => alpha=1 - # - # svp is interpolated between the ice and water phases: - # svp = alpha * es_water + (1.0 - alpha) * es_ice - - t = np.asarray(t) - svp = np.zeros(t.shape) - - # ice range - i_mask = t <= self.ti - svp[i_mask] = self._es_ice(t[i_mask]) - - # water range - w_mask = t >= self.t0 - svp[w_mask] = self._es_water(t[w_mask]) - - # mixed range - m_mask = ~(i_mask | w_mask) - alpha = np.square(t[m_mask] - self.ti) / np.square(self.t0 - self.ti) - svp[m_mask] = alpha * self._es_water(t[m_mask]) + (1.0 - alpha) * self._es_ice( - t[m_mask] - ) - return svp - - def _es_water_slope(self, t): - return ( - self._es_water(t) - * (self.c3w * (self.t0 - self.c4w)) - / np.square(t - self.c4w) - ) - - def _es_ice_slope(self, t): - return ( - self._es_ice(t) - * (self.c3i * (self.t0 - self.c4i)) - / np.square(t - self.c4i) - ) - - def _es_mixed_slope(self, t): - t = np.asarray(t) - d_svp = np.zeros(t.shape) - - # ice range - i_mask = t <= self.ti - d_svp[i_mask] = self._es_ice_slope(t[i_mask]) - - # water range - w_mask = t >= self.t0 - d_svp[w_mask] = self._es_water_slope(t[w_mask]) - - # mixed range - m_mask = ~(i_mask | w_mask) - alpha = np.square(t[m_mask] - self.ti) / np.square(self.t0 - self.ti) - d_alpha = (2.0 / np.square(self.t0 - self.ti)) * (t[m_mask] - self.ti) - t_m = t[m_mask] - d_svp[m_mask] = ( - d_alpha * self._es_water(t_m) - + alpha * self._es_water_slope(t_m) - - d_alpha * self._es_ice(t_m) - + (1.0 - alpha) * self._es_ice_slope(t_m) - ) - return d_svp - - -def saturation_vapour_pressure(t, phase="mixed"): - r"""Computes the saturation vapour pressure from temperature with respect to a phase. - - Parameters - ---------- - t: ndarray - Temperature (K) - phase: str, optional - Defines the phase with respect to the saturation vapour pressure is computed. - It is either “water”, “ice” or “mixed”. - - Returns - ------- - ndarray - Saturation vapour pressure (Pa) - - - The algorithm was taken from the IFS model [IFS-CY47R3-PhysicalProcesses]_ (see Chapter 12). - It uses the following formula when ``phase`` is "water" or "ice": - - .. math:: - - e_{sat} = a_{1}\;exp \left(a_{3}\frac{t-273.16}{t-a_{4}}\right) - - where the parameters are set as follows: - - * ``phase`` = "water": :math:`a_{1}` =611.21 Pa, :math:`a_{3}` =17.502 and :math:`a_{4}` =32.19 K - * ``phase`` = "ice": :math:`a_{1}` =611.21 Pa, :math:`a_{3}` =22.587 and :math:`a_{4}` =-0.7 K - - When ``phase`` is "mixed" the formula is based on the value of ``t``: - - * if :math:`t <= t_{i}`: the formula for ``phase`` = "ice" is used (:math:`t_{i} = 250.16 K`) - * if :math:`t >= t_{0}`: the formula for ``phase`` = "water" is used (:math:`t_{0} = 273.16 K`) - * for the range :math:`t_{i} < t < t_{0}` an interpolation is used between the "ice" and "water" phases: - - .. math:: - - \alpha(t) e_{wsat}(t) + (1 - \alpha(t)) e_{isat}(t) - - with :math:`\alpha(t) = (\frac{t-t_{i}}{t_{0}-t_{i}})^2`. - - """ - return _EsComp().compute_es(t, phase) - - -def saturation_mixing_ratio(t, p, phase="mixed"): - r"""Computes the saturation mixing ratio from temperature with respect to a phase. - - Parameters - ---------- - t: ndarray - Temperature (K) - p: ndarray - Pressure (Pa) - phase: str - Defines the phase with respect to the :func:`saturation_vapour_pressure` is computed. - It is either “water”, “ice” or “mixed”. - - Returns - ------- - ndarray - Saturation mixing ratio (kg/kg) - - - Equivalent to the following code: - - .. code-block:: python - - e = saturation_vapour_pressure(t, phase=phase) - return mixing_ratio_from_vapour_pressure(e, p) - - """ - e = saturation_vapour_pressure(t, phase=phase) - return mixing_ratio_from_vapour_pressure(e, p) - - -def saturation_specific_humidity(t, p, phase="mixed"): - r"""Computes the saturation specific humidity from temperature with respect to a phase. - - Parameters - ---------- - t: ndarray - Temperature (K) - p: ndarray - Pressure (Pa) - phase: str, optional - Defines the phase with respect to the :func:`saturation_vapour_pressure` is computed. - It is either “water”, “ice” or “mixed”. - - Returns - ------- - ndarray - Saturation specific humidity (kg/kg) - - - Equivalent to the following code: - - .. code-block:: python - - e = saturation_vapour_pressure(t, phase=phase) - return specific_humidity_from_vapour_pressure(e, p) - - """ - e = saturation_vapour_pressure(t, phase=phase) - return specific_humidity_from_vapour_pressure(e, p) - - -def saturation_vapour_pressure_slope(t, phase="mixed"): - r"""Computes the slope of saturation vapour pressure with respect to temperature. - - Parameters - ---------- - t: ndarray - Temperature (K) - phase: str, optional - Defines the phase with respect to the computation will be performed. - It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` - for details. - - Returns - ------- - ndarray - Slope of saturation vapour pressure (Pa/K) - - """ - return _EsComp().compute_slope(t, phase) - - -def saturation_mixing_ratio_slope( - t, p, es=None, es_slope=None, phase="mixed", eps=1e-4 -): - r"""Computes the slope of saturation mixing ratio with respect to temperature. - - Parameters - ---------- - t: ndarray - Temperature (K) - p: ndarray - Pressure (Pa) - es: ndarray or None, optional - :func:`saturation_vapour_pressure` pre-computed for the given ``phase`` (Pa) - es_slope: ndarray or None, optional - :func:`saturation_vapour_pressure_slope` pre-computed for the given ``phase`` (Pa/K) - phase: str, optional - Defines the phase with respect to the computation will be performed. - It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` - for details. - eps: number - Where p - es < ``eps`` np.nan is returned. - - Returns - ------- - ndarray - Slope of saturation mixing ratio (:math:`kg kg^{-1} K^{-1}`) - - - The computation is based on the following formula: - - .. math:: - - \frac{\partial w_{s}}{\partial t} = \frac{\epsilon\; p}{(p-e_{s})^{2}} \frac{d e_{s}}{d t} - - where - - * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - * :math:`e_{s}` is the :func:`saturation_vapour_pressure` for the given ``phase`` - - """ - if eps <= 0: - raise ValueError(f"saturation_mixing_ratio_slope(): eps={eps} must be > 0") - if es is None: - es = saturation_vapour_pressure(t, phase=phase) - if es_slope is None: - es_slope = saturation_vapour_pressure_slope(t, phase=phase) - - v = np.asarray(p - es) - v[v < eps] = np.nan - return constants.epsilon * es_slope * p / np.square(v) - - -def saturation_specific_humidity_slope( - t, p, es=None, es_slope=None, phase="mixed", eps=1e-4 -): - r"""Computes the slope of saturation specific humidity with respect to temperature. - - Parameters - ---------- - t: ndarray - Temperature (K) - p: ndarray - Pressure (Pa) - es: ndarray or None, optional - :func:`saturation_vapour_pressure` pre-computed for the given ``phase`` (Pa) - es_slope: ndarray or None, optional - :func:`saturation_vapour_pressure_slope` pre-computed for the given ``phase`` (Pa/K) - phase: str, optional - Defines the phase with respect to the computation will be performed. - It is either “water”, “ice” or “mixed”. See :func:`saturation_vapour_pressure` - for details. - eps: number - Where p - es < ``eps`` np.nan is returned. - - Returns - ------- - ndarray - Slope of saturation specific humidity (:math:`kg kg^{-1} K^{-1}`) - - - The computation is based on the following formula: - - .. math:: - - \frac{\partial q_{s}}{\partial t} = - \frac{\epsilon\; p}{(p+e_{s}(\epsilon - 1))^{2}} \frac{d e_{s}}{d t} - - where - - * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - * :math:`e_{s}` is the :func:`saturation_vapour_pressure` for the given ``phase`` - - """ - if eps <= 0: - raise ValueError(f"saturation_specific_humidity_slope(): eps={eps} must be > 0") - if es is None: - es = saturation_vapour_pressure(t, phase=phase) - if es_slope is None: - es_slope = saturation_vapour_pressure_slope(t, phase=phase) - - v = np.asarray(np.square(p + es * (constants.epsilon - 1.0))) - v[np.asarray(p - es) < eps] = np.nan - return constants.epsilon * es_slope * p / v - - -def temperature_from_saturation_vapour_pressure(es): - r"""Computes the temperature from saturation vapour pressure. - - Parameters - ---------- - es: ndarray - :func:`saturation_vapour_pressure` (Pa) - - Returns - ------- - ndarray - Temperature (K) - - - The computation is always based on the "water" phase of - the :func:`saturation_vapour_pressure` formulation irrespective of the - phase ``es`` was computed to. - - """ - return _EsComp().t_from_es(es) - - -def relative_humidity_from_dewpoint(t, td): - r"""Computes the relative humidity from dewpoint temperature. - - Parameters - ---------- - t: ndarray - Temperature (K) - td: ndarray - Dewpoint (K) - - Returns - ------- - ndarray - Relative humidity (%) - - - The computation is based on the following formula: - - .. math:: - - r = 100 \frac {e_{wsat}(td)}{e_{wsat}(t)} - - where :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water. - - """ - e = saturation_vapour_pressure(td, phase="water") - es = saturation_vapour_pressure(t, phase="water") - return 100.0 * e / es - - -def relative_humidity_from_specific_humidity(t, q, p): - r"""Computes the relative humidity from specific humidity. - - Parameters - ---------- - t: ndarray - Temperature (K) - q: ndarray - Specific humidity (kg/kg) - p: ndarray - Pressure (Pa) - - Returns - ------- - ndarray - Relative humidity (%) - - - The computation is based on the following formula: - - .. math:: - - r = 100 \frac {e(q, p)}{e_{msat}(t)} - - where: - - * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) - * :math:`e_{msat}` is the :func:`saturation_vapour_pressure` based on the "mixed" phase - - """ - svp = saturation_vapour_pressure(t) - e = vapour_pressure_from_specific_humidity(q, p) - return 100.0 * e / svp - - -def specific_humidity_from_dewpoint(td, p): - r"""Computes the specific humidity from dewpoint. - - Parameters - ---------- - td: ndarray - Dewpoint (K) - p: ndarray - Pressure (Pa) - - Returns - ------- - ndarray - Specific humidity (kg/kg) - - - The computation starts with determining the vapour pressure: - - .. math:: - - e(q, p) = e_{wsat}(td) - - where: - - * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) - * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water - * :math:`q` is the specific humidity - - Then `q` is computed from :math:`e` using :func:`specific_humidity_from_vapour_pressure`. - - """ - svp = saturation_vapour_pressure(td, phase="water") - return specific_humidity_from_vapour_pressure(svp, p) - - -def mixing_ratio_from_dewpoint(td, p): - r"""Computes the mixing ratio from dewpoint. - - Parameters - ---------- - td: ndarray - Dewpoint (K) - p: ndarray - Pressure (Pa) - - Returns - ------- - ndarray - Specific humidity (kg/kg) - - - The computation starts with determining the vapour pressure: - - .. math:: - - e(w, p) = e_{wsat}(td) - - where: - - * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_mixing_ratio`) - * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water - * :math:`w` is the mixing ratio - - Then `w` is computed from :math:`e` using :func:`mixing_ratio_from_vapour_pressure`. - - """ - svp = saturation_vapour_pressure(td, phase="water") - return mixing_ratio_from_vapour_pressure(svp, p) - - -def specific_humidity_from_relative_humidity(t, r, p): - r"""Computes the specific humidity from relative_humidity. - - Parameters - ---------- - t: ndarray - Temperature (K) - r: ndarray - Relative humidity(%) - p: ndarray - Pressure (Pa) - - Returns - ------- - ndarray - Specific humidity (kg/kg) units - - - The computation starts with determining the the vapour pressure: - - .. math:: - - e(q, p) = r\; \frac{e_{msat}(t)}{100} - - where: - - * :math:`e` is the vapour pressure (see :func:`vapour_pressure`) - * :math:`e_{msat}` is the :func:`saturation_vapour_pressure` based on the "mixed" phase - * :math:`q` is the specific humidity - - Then :math:`q` is computed from :math:`e` using :func:`specific_humidity_from_vapour_pressure`. - - """ - e = r * saturation_vapour_pressure(t) / 100.0 - return specific_humidity_from_vapour_pressure(e, p) - - -def dewpoint_from_relative_humidity(t, r): - r"""Computes the dewpoint temperature from relative humidity. - - Parameters - ---------- - t: ndarray - Temperature (K) - r: ndarray - Relative humidity (%) - - Returns - ------- - ndarray - Dewpoint (K) - - - The computation starts with determining the the saturation vapour pressure over - water at the dewpoint temperature: - - .. math:: - - e_{wsat}(td) = \frac{r\; e_{wsat}(t)}{100} - - where: - - * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water - * :math:`td` is the dewpoint. - - Then :math:`td` is computed from :math:`e_{wsat}(td)` by inverting the - equations used in :func:`saturation_vapour_pressure`. - - """ - es = saturation_vapour_pressure(t, phase="water") * r / 100.0 - return temperature_from_saturation_vapour_pressure(es) - - -def dewpoint_from_specific_humidity(q, p): - r"""Computes the dewpoint temperature from specific humidity. - - Parameters - ---------- - q: ndarray - Specific humidity (kg/kg) - p: ndarray - Pressure (Pa) - - Returns - ------- - ndarray - Dewpoint (K) - - - The computation starts with determining the the saturation vapour pressure over - water at the dewpoint temperature: - - .. math:: - - e_{wsat}(td) = e(q, p) - - where: - - * :math:`e` is the vapour pressure (see :func:`vapour_pressure_from_specific_humidity`) - * :math:`e_{wsat}` is the :func:`saturation_vapour_pressure` over water - * :math:`td` is the dewpoint - - Then :math:`td` is computed from :math:`e_{wsat}(td)` by inverting the equations - used in :func:`saturation_vapour_pressure`. - - """ - return temperature_from_saturation_vapour_pressure( - vapour_pressure_from_specific_humidity(q, p) - ) - - -def virtual_temperature(t, q): - r"""Computes the virtual temperature from temperature and specific humidity. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - q: number or ndarray - Specific humidity (kg/kg) - - Returns - ------- - number or ndarray - Virtual temperature (K) - - - The computation is based on the following formula [Wallace2006]_: - - .. math:: - - t_{v} = t (1 + \frac{1 - \epsilon}{\epsilon} q) - - with :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - c1 = (1.0 - constants.epsilon) / constants.epsilon - return t * (1.0 + c1 * q) - - -def virtual_potential_temperature(t, q, p): - r"""Computes the virtual potential temperature from temperature and specific humidity. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - q: number or ndarray - Specific humidity (kg/kg) - p: number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Virtual potential temperature (K) - - - The computation is based on the following formula: - - .. math:: - - \Theta_{v} = \theta (1 + \frac{1 - \epsilon}{\epsilon} q) - - where: - - * :math:`\Theta` is the :func:`potential_temperature` - * :math:`\epsilon = R_{d}/R_{v}` (see :data:`earthkit.meteo.constants.epsilon`). - - """ - c1 = (1.0 - constants.epsilon) / constants.epsilon - return potential_temperature(t, p) * (1.0 + c1 * q) - - -def potential_temperature(t, p): - r"""Computes the potential temperature. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - p: number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Potential temperature (K) - - - The computation is based on the following formula [Wallace2006]_: - - .. math:: - - \theta = t (\frac{10^{5}}{p})^{\kappa} - - with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). - - """ - return t * np.power(constants.p0 / p, constants.kappa) - - -def temperature_from_potential_temperature(th, p): - r"""Computes the temperature from potential temperature. - - Parameters - ---------- - th: number or ndarray - Potential temperature (K) - p: number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Temperature (K) - - - The computation is based on the following formula: - - .. math:: +from . import array - t = \theta (\frac{p}{10^{5}})^{\kappa} - with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). +def celsius_to_kelvin(*args, **kwargs): + return array.celsius_to_kelvin(*args, **kwargs) - """ - return th * np.power(p / constants.p0, constants.kappa) +def kelvin_to_celsius(*args, **kwargs): + return array.kelvin_to_celsius(*args, **kwargs) -def pressure_on_dry_adiabat(t, t_def, p_def): - r"""Computes the pressure on a dry adiabat. - Parameters - ---------- - t: number or ndarray - Temperature on the dry adiabat (K) - t_def: number or ndarray - Temperature defining the dry adiabat (K) - p_def: number or ndarray - Pressure defining the dry adiabat (Pa) +def specific_humidity_from_mixing_ratio(*args, **kwargs): + return array.specific_humidity_from_mixing_ratio(*args, **kwargs) - Returns - ------- - number or ndarray - Pressure on the dry adiabat (Pa) +def mixing_ratio_from_specific_humidity(*args, **kwargs): + return array.mixing_ratio_from_specific_humidity(*args, **kwargs) - The computation is based on the following formula: - .. math:: +def vapour_pressure_from_specific_humidity(*args, **kwargs): + return array.vapour_pressure_from_specific_humidity(*args, **kwargs) - p = p_{def} (\frac{t}{t_{def}})^{\frac{1}{\kappa}} - with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). +def vapour_pressure_from_mixing_ratio(*args, **kwargs): + return array.vapour_pressure_from_mixing_ratio(*args, **kwargs) - """ - return p_def * np.power(t / t_def, 1 / constants.kappa) +def specific_humidity_from_vapour_pressure(*args, **kwargs): + return array.specific_humidity_from_vapour_pressure(*args, **kwargs) -def temperature_on_dry_adiabat(p, t_def, p_def): - r"""Computes the temperature on a dry adiabat. - Parameters - ---------- - p: number or ndarray - Pressure on the dry adiabat (Pa) - t_def: number or ndarray - Temperature defining the dry adiabat (K) - p_def: number or ndarray - Pressure defining the dry adiabat (Pa) +def mixing_ratio_from_vapour_pressure(*args, **kwargs): + return array.mixing_ratio_from_vapour_pressure(*args, **kwargs) - Returns - ------- - number or ndarray - Temperature on the dry adiabat (K) +def saturation_vapour_pressure(*args, **kwargs): + return array.saturation_vapour_pressure(*args, **kwargs) - The computation is based on the following formula: - .. math:: +def saturation_mixing_ratio(*args, **kwargs): + return array.saturation_mixing_ratio(*args, **kwargs) - t = t_{def} (\frac{p}{p_{def}})^{\kappa} - with :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`). +def saturation_specific_humidity(*args, **kwargs): + return array.saturation_specific_humidity(*args, **kwargs) - """ - return t_def * np.power(p / p_def, constants.kappa) +def saturation_vapour_pressure_slope(*args, **kwargs): + return array.saturation_vapour_pressure_slope(*args, **kwargs) -def lcl_temperature(t, td, method="davies"): - r"""Computes the Lifting Condenstaion Level (LCL) temperature from dewpoint. - Parameters - ---------- - t: number or ndarray - Temperature at the start level (K) - td: number or ndarray - Dewpoint at the start level (K) - method: str, optional - The computation method: "davies" or "bolton". +def saturation_mixing_ratio_slope(*args, **kwargs): + return array.saturation_mixing_ratio_slope(*args, **kwargs) - Returns - ------- - number or ndarray - Temperature of the LCL (K) +def saturation_specific_humidity_slope(*args, **kwargs): + return array.saturation_specific_humidity_slope(*args, **kwargs) - The actual computation is based on the ``method``: - * "davies": the formula by [DaviesJones1983]_ is used (it is also used by the IFS model): +def temperature_from_saturation_vapour_pressure(*args, **kwargs): + return array.temperature_from_saturation_vapour_pressure(*args, **kwargs) - .. math:: - t_{LCL} = - td - (0.212 + 1.571\times 10^{-3} (td - t_{0}) - 4.36\times 10^{-4} (t - t_{0})) (t - td) +def relative_humidity_from_dewpoint(*args, **kwargs): + return array.relative_humidity_from_dewpoint(*args, **kwargs) - where :math:`t_{0}` is the triple point of water (see :data:`earthkit.meteo.constants.T0`). - * "bolton": the formula by [Bolton1980]_ is used: +def relative_humidity_from_specific_humidity(*args, **kwargs): + return array.relative_humidity_from_specific_humidity(*args, **kwargs) - .. math:: - t_{LCL} = 56.0 + \frac{1}{\frac{1}{td - 56} + \frac{log(\frac{t}{td})}{800}} +def specific_humidity_from_dewpoint(*args, **kwargs): + return array.specific_humidity_from_dewpoint(*args, **kwargs) - """ - # Davies-Jones formula - if method == "davies": - t_lcl = td - ( - 0.212 + 1.571e-3 * (td - constants.T0) - 4.36e-4 * (t - constants.T0) - ) * (t - td) - return t_lcl - # Bolton formula - elif method == "bolton": - return 56.0 + 1 / (1 / (td - 56) + np.log(t / td) / 800) - else: - raise ValueError(f"lcl_temperature: invalid method={method} specified!") +def mixing_ratio_from_dewpoint(*args, **kwargs): + return array.mixing_ratio_from_dewpoint(*args, **kwargs) -def lcl(t, td, p, method="davies"): - r"""Computes the temperature and pressure of the Lifting Condenstaion Level (LCL) from dewpoint. - Parameters - ---------- - t: number or ndarray - Temperature at the start level (K) - td: number or ndarray - Dewpoint at the start level (K) - p: number or ndarray - Pressure at the start level (Pa) - method: str - method: str, optional - The computation method: "davies" or "bolton". - - Returns - ------- - number or ndarray - Temperature of the LCL (K) - number or ndarray - Pressure of the LCL (Pa) - - - The LCL temperature is determined by :func:`lcl_temperature` with the given ``method`` - and the pressure is computed with :math:`t_{LCL}` using :func:`pressure_on_dry_adiabat`. - - """ - t_lcl = lcl_temperature(t, td, method=method) - p_lcl = pressure_on_dry_adiabat(t_lcl, t, p) - return t_lcl, p_lcl - - -class _ThermoState: - def __init__(self, t=None, td=None, q=None, w=None, p=None): - self.t = t - self.td = td - self.q = q - self.w = w - self.p = p - self.e = None - self.es = None - self.ws = None - self.qs = None - - -class _EptComp: - CM = {} - c_lambda = 1.0 / constants.kappa - - @staticmethod - def make(method): - return _EptComp.CM[method]() - - def is_mixing_ratio_based(self): - return True - - def compute_ept(self, t=None, td=None, q=None, p=None): - if td is None and q is None: - raise ValueError("ept: either td or q must have a valid value!") - - # compute dewpoint since all the methods require it - if td is None: - td = dewpoint_from_specific_humidity(q, p) - - ths = _ThermoState(t=t, td=td, q=q, p=p) - return self._ept(ths) - - def compute_ept_sat(self, t, p): - ths = _ThermoState(t=t, p=p) - return self._th_sat(ths) * np.exp(self._G_sat(ths)) - - def compute_wbpt(self, ept): - t0 = 273.16 - x = ept / t0 - a = [7.101574, -20.68208, 16.11182, 2.574631, -5.205688] - b = [1.0, -3.552497, 3.781782, -0.6899655, -0.5929340] - return ept - np.exp( - np.polynomial.polynomial.polyval(x, a) - / np.polynomial.polynomial.polyval(x, b) - ) - - def compute_t_on_ma_stipanuk(self, ept, p): - if isinstance(p, np.ndarray): - t = np.full(p.shape, constants.T0 - 20) - elif isinstance(ept, np.ndarray): - t = np.full(ept.shape, constants.T0 - 20) - else: - t = constants.T0 - 20 - max_iter = 12 - dt = 120.0 - for _ in range(max_iter): - ths = _ThermoState(t=t, p=p) - dt /= 2.0 - t += ( - np.sign(ept * np.exp(self._G_sat(ths, scale=-1.0)) - self._th_sat(ths)) - * dt - ) - # ths.t = t - # return ept - self._th_sat(ths) * np.exp(self._G_sat(ths)) - return t - - def compute_t_on_ma_davies(self, ept, p): - def _k1(pp): - """Function k1 in the article.""" - a = [-53.737, 137.81, -38.5] - return np.polynomial.polynomial.polyval(pp, a) - - def _k2(pp): - """Function k2 in the article.""" - a = [-0.384, 56.831, -4.392] - return np.polynomial.polynomial.polyval(pp, a) - - def _D(p): - """Function D in the article.""" - return 1.0 / (0.1859e-5 * p + 0.6512) - - max_iter = 1 - A = 2675 - t0 = 273.16 - - if not isinstance(p, np.ndarray): - p = np.full(ept.shape, p) - - tw = ept.copy() - pp = np.power(p / constants.p0, constants.kappa) - te = ept * pp - c_te = np.power(t0 / te, _EptComp.c_lambda) - - # initial guess - mask = c_te > _D(p) - if np.any(mask): - es = saturation_vapour_pressure(te[mask]) - ws = mixing_ratio_from_vapour_pressure(es, p[mask]) - d_es = saturation_vapour_pressure_slope(te[mask]) - tw[mask] = te[mask] - t0 - (A * ws) / (1 + A * ws * d_es / es) - - mask = (1 <= c_te) & (c_te <= _D(p)) - tw[mask] = _k1(pp[mask]) - _k2(pp[mask]) * c_te[mask] - - mask = (0.4 <= c_te) & (c_te < 1) - tw[mask] = (_k1(pp[mask]) - 1.21) - (_k2(pp[mask]) - 1.21) * c_te[mask] - - mask = c_te < 0.4 - tw[mask] = ( - (_k1(pp[mask]) - 2.66) - - (_k2(pp[mask]) - 1.21) * c_te[mask] - + 0.58 / c_te[mask] - ) - # tw has to be converted to Kelvin - tw = celsius_to_kelvin(tw) - - # np.place(tw, tw > 362.16, 362.1) - # print(f"ept={ept}") - # print(f"tw={tw}") - - for i in range(max_iter): - ths = _ThermoState(t=tw, p=p) - ths.c_tw = np.power(t0 / ths.t, _EptComp.c_lambda) - ths.es = saturation_vapour_pressure(ths.t) - if self.is_mixing_ratio_based(): - ths.ws = mixing_ratio_from_vapour_pressure(ths.es, ths.p) - else: - ths.qs = specific_humidity_from_vapour_pressure(ths.es, ths.p) - f_val = self._f(ths) - # print(f"{i}") - # print(f" p={p}") - # print(f" t={ths.t}") - # print(f" f_val={f_val}") - # print(f" f_val_cte={f_val - c_te}") - # print(f" d_lnf={self._d_lnf(ths)}") - # print(tw[mask].shape) - # print(f_val.shape) - tw -= (f_val - c_te) / (f_val * self._d_lnf(ths)) - # print(f"tw={tw}") - - # when ept is derived with the Bolton methods the iteration breaks down for extremeley - # hot and humid conditions (roughly about t > 50C and r > 95%) and results in negative - # tw values! - np.place(tw, tw <= 0, np.nan) - - # ths = _ThermoState(t=tw, p=p) - # return ept - self._th_sat(ths) * np.exp(self._G_sat(ths)) - return tw - - -class _EptCompIfs(_EptComp): - def __init__(self): - self.K0 = constants.Lv / constants.c_pd - - def is_mixing_ratio_based(self): - return False - - def _ept(self, ths): - th = potential_temperature(ths.t, ths.p) - t_lcl = lcl_temperature(ths.t, ths.td, method="davies") - if ths.q is None: - ths.q = specific_humidity_from_dewpoint(ths.td, ths.p) - return th * np.exp(self.K0 * ths.q / t_lcl) - - def _th_sat(self, ths): - return potential_temperature(ths.t, ths.p) - - def _G_sat(self, ths, scale=1.0): - qs = saturation_specific_humidity(ths.t, ths.p) - return (scale * self.K0) * qs / ths.t - - def _d_G_sat(self, ths): - if ths.qs is None: - ths.qs = saturation_specific_humidity(ths.t, ths.p) - return ( - -self.K0 * ths.qs / (ths.t**2) - + self.K0 * saturation_specific_humidity_slope(ths.t, ths.p) / ths.t - ) - - def _f(self, ths): - return ths.c_tw * np.exp(self._G_sat(ths, scale=-self.c_lambda)) - - def _d_lnf(self, ths): - return -self.c_lambda * (1 / ths.t + self._d_G_sat(ths)) - - -class _EptCompBolton35(_EptComp): - def __init__(self): - self.K0 = 2675.0 - self.K3 = 0.28 - - def _ept(self, ths): - t_lcl = lcl_temperature(ths.t, ths.td, method="bolton") - if ths.q is None: - w = mixing_ratio_from_dewpoint(ths.td, ths.p) - else: - w = mixing_ratio_from_specific_humidity(ths.q) - th = ths.t * np.power(constants.p0 / ths.p, constants.kappa * (1 - self.K3 * w)) - return th * np.exp(self.K0 * w / t_lcl) - - def _th_sat(self, ths): - if ths.ws is None: - ths.ws = saturation_mixing_ratio(ths.t, ths.p) - return ths.t * np.power( - constants.p0 / ths.p, constants.kappa * (1 - self.K3 * ths.ws) - ) - - def _G_sat(self, ths, scale=1.0): - if ths.ws is None: - ths.ws = saturation_mixing_ratio(ths.t, ths.p) - return (scale * self.K0) * ths.ws / ths.t - - def _d_G_sat(self, ths): - return ( - -self.K0 * ths.ws / np.square(ths.t) - + self.K0 * saturation_mixing_ratio_slope(ths.t, ths.p) / ths.t - ) - - def _f(self, ths): - # print(f" c_tw={ths.c_tw}") - # print(f" es_frac={ths.es / ths.p}") - # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}") - return ( - ths.c_tw - * np.power(ths.p / constants.p0, self.K3 * ths.ws) - * np.exp(self._G_sat(ths, scale=-self.c_lambda)) - ) - - def _d_lnf(self, ths): - return -self.c_lambda * ( - 1 / ths.t - + self.K3 - * np.log(ths.p / constants.p0) - * saturation_vapour_pressure_slope(ths.t) - + self._d_G_sat(ths) - ) - - -class _EptCompBolton39(_EptComp): - def __init__(self): - # Comment on the Bolton formulas. The constants used by Bolton to - # derive the formulas differ from the ones used by earthkit.meteo E.g. - # Bolton earthkit.meteo - # Rd 287.04 287.0597 - # cpd 1005.7 1004.79 - # kappa 0.2854 0.285691 - # epsilon 0.6220 0.621981 - - self.K0 = 3036.0 - self.K1 = 1.78 - self.K2 = 0.448 - self.K4 = 0.28 - - def _ept(self, ths): - t_lcl = lcl_temperature(ths.t, ths.td, method="bolton") - if ths.q is None: - w = mixing_ratio_from_dewpoint(ths.td, ths.p) - else: - w = mixing_ratio_from_specific_humidity(ths.q) - - e = vapour_pressure_from_mixing_ratio(w, ths.p) - th = potential_temperature(ths.t, ths.p - e) * np.power( - ths.t / t_lcl, self.K4 * w - ) - return th * np.exp((self.K0 / t_lcl - self.K1) * w * (1.0 + self.K2 * w)) - - def _th_sat(self, ths): - if ths.es is None: - ths.es = saturation_vapour_pressure(ths.t) - ths.es[ths.p - ths.es < 1e-4] = np.nan - return potential_temperature(ths.t, ths.p - ths.es) - - def _G_sat(self, ths, scale=1.0): - if ths.es is None: - ths.es = saturation_vapour_pressure(ths.t) - ths.es[ths.p - ths.es < 1e-4] = np.nan - # print(f" es={ths.es}") - ws = mixing_ratio_from_vapour_pressure(ths.es, ths.p) - # print(f" ws={ws}") - return ( - ((scale * self.K0) / ths.t - (scale * self.K1)) * ws * (1.0 + self.K2 * ws) - ) - - def _d_G_sat(self, ths): - # print(f" d_ws={saturation_mixing_ratio_slope(ths.t, ths.p)}") - return -self.K0 * (ths.ws + self.K2 * np.square(ths.ws)) / ( - np.square(ths.t) - ) + (self.K0 / ths.t - self.K1) * ( - 1 + (2 * self.K2) * ths.ws - ) * saturation_mixing_ratio_slope( - ths.t, ths.p - ) - - def _f(self, ths): - # print(f" c_tw={ths.c_tw}") - # print(f" es_frac={ths.es / ths.p}") - # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}") - return ( - ths.c_tw - * (1 - ths.es / ths.p) - * np.exp(self._G_sat(ths, scale=-self.c_lambda)) - ) - - def _d_lnf(self, ths): - return -self.c_lambda * ( - 1 / ths.t - + constants.kappa - * saturation_vapour_pressure_slope(ths.t) - / (ths.p - ths.es) - + self._d_G_sat(ths) - ) - - -_EptComp.CM = { - "ifs": _EptCompIfs, - "bolton35": _EptCompBolton35, - "bolton39": _EptCompBolton39, -} - - -def ept_from_dewpoint(t, td, p, method="ifs"): - r"""Computes the equivalent potential temperature from dewpoint. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - td: number or ndarray - Dewpoint (K) - p: number or ndarray - Pressure (Pa) - method: str, optional - Specifies the computation method. The possible values are: "ifs", "bolton35", "bolton39". - - Returns - ------- - number or ndarray - Equivalent potential temperature (K) - - - The actual computation is based on the value of ``method``: - - * "ifs": the formula from the IFS model [IFS-CY47R3-PhysicalProcesses]_ (Chapter 6.11) is used: - - .. math:: +def specific_humidity_from_relative_humidity(*args, **kwargs): + return array.specific_humidity_from_relative_humidity(*args, **kwargs) - \Theta_{e} = \Theta\; exp(\frac{L_{v}\; q}{c_{pd}\; t_{LCL}}) - * "bolton35": Eq (35) from [Bolton1980]_ is used: +def dewpoint_from_relative_humidity(*args, **kwargs): + return array.dewpoint_from_relative_humidity(*args, **kwargs) - .. math:: +def dewpoint_from_specific_humidity(*args, **kwargs): + return array.dewpoint_from_specific_humidity(*args, **kwargs) - \Theta_{e} = \Theta (\frac{10^{5}}{p})^{\kappa 0.28 w} exp(\frac{2675 w}{t_{LCL}}) - * "bolton39": Eq (39) from [Bolton1980]_ is used: +def virtual_temperature(*args, **kwargs): + return array.virtual_temperature(*args, **kwargs) - .. math:: - \Theta_{e} = - t (\frac{10^{5}}{p-e})^{\kappa} (\frac{t}{t_{LCL}})^{0.28 w} exp[(\frac{3036}{t_{LCL}} - - 1.78)w(1+0.448\; w)] +def virtual_potential_temperature(*args, **kwargs): + return array.virtual_potential_temperature(*args, **kwargs) - where: - * :math:`\Theta` is the :func:`potential_temperature` - * :math:`t_{LCL}` is the temperature at the Lifting Condestation Level computed - with :func:`lcl_temperature` using option: +def potential_temperature(*args, **kwargs): + return array.potential_temperature(*args, **kwargs) - * method="davis" when ``method`` is "ifs" - * method="bolton" when ``method`` is "bolton35" or "bolton39" - * :math:`q` is the specific humidity computed with :func:`specific_humidity_from_dewpoint` - * :math:`w`: is the mixing ratio computed with :func:`mixing_ratio_from_dewpoint` - * :math:`e` is the vapour pressure computed with :func:`vapour_pressure_from_mixing_ratio` - * :math:`L_{v}`: is the latent heat of vaporisation - (see :data:`earthkit.meteo.constants.Lv`) - * :math:`c_{pd}` is the specific heat of dry air on constant pressure - (see :data:`earthkit.meteo.constants.c_pd`) - * :math:`\kappa = R_{d}/c_{pd}` (see :data:`earthkit.meteo.constants.kappa`) - """ - return _EptComp.make(method).compute_ept(t=t, td=td, p=p) +def temperature_from_potential_temperature(*args, **kwargs): + return array.temperature_from_potential_temperature(*args, **kwargs) -def ept_from_specific_humidity(t, q, p, method="ifs"): - r"""Computes the equivalent potential temperature from specific humidity. +def pressure_on_dry_adiabat(*args, **kwargs): + return array.pressure_on_dry_adiabat(*args, **kwargs) - Parameters - ---------- - t: number or ndarray - Temperature (K) - q: number or ndarray - Specific humidity (kg/kg) - p: number or ndarray - Pressure (Pa) - method: str, optional - Specifies the computation method. The possible values are: "ifs", - "bolton35", "bolton39. See :func:`ept_from_dewpoint` for details. - Returns - ------- - number or ndarray - Equivalent potential temperature (K) +def temperature_on_dry_adiabat(*args, **kwargs): + return array.temperature_on_dry_adiabat(*args, **kwargs) - The computations are the same as in :func:`ept_from_dewpoint` - (the dewpoint is computed from q with :func:`dewpoint_from_specific_humidity`). +def lcl_temperature(*args, **kwargs): + return array.lcl_temperature(*args, **kwargs) - """ - return _EptComp.make(method).compute_ept(t=t, q=q, p=p) +def lcl(*args, **kwargs): + return array.lcl(*args, **kwargs) -def saturation_ept(t, p, method="ifs"): - r"""Computes the saturation equivalent potential temperature. - Parameters - ---------- - t: number or ndarray - Temperature (K) - p: number or ndarray - Pressure (Pa) - method: str, optional - Specifies the computation method. The possible values are: "ifs", "bolton35", "bolton39". +def ept_from_dewpoint(*args, **kwargs): + return array.ept_from_dewpoint(*args, **kwargs) - Returns - ------- - number or ndarray - Saturation equivalent potential temperature (K) +def ept_from_specific_humidity(*args, **kwargs): + return array.ept_from_specific_humidity(*args, **kwargs) - The actual computation is based on the ``method``: - * "ifs": The formula is based on the equivalent potential temperature definition used - in the IFS model [IFS-CY47R3-PhysicalProcesses]_ (see Chapter 6.11) : +def saturation_ept(*args, **kwargs): + return array.saturation_ept(*args, **kwargs) - .. math:: - \Theta_{esat} = \Theta\; exp(\frac{L_{v}\; q_{sat}}{c_{pd}\; t}) +def temperature_on_moist_adiabat(*args, **kwargs): + return array.temperature_on_moist_adiabat(*args, **kwargs) - * "bolton35": Eq (35) from [Bolton1980]_ is used: - .. math:: +def wet_bulb_temperature_from_dewpoint(*args, **kwargs): + return array.wet_bulb_temperature_from_dewpoint(*args, **kwargs) - \Theta_{e} = \Theta (\frac{10^{5}}{p})^{\kappa 0.28 w_{sat}}\; exp(\frac{2675\; w_{sat}}{t}) - * "bolton39": Eq (39) from [Bolton1980]_ is used: +def wet_bulb_temperature_from_specific_humidity(*args, **kwargs): + return array.wet_bulb_temperature_from_specific_humidity(*args, **kwargs) - .. math:: - \Theta_{e} = - t (\frac{10^{5}}{p-e_{sat}})^{\kappa} exp[(\frac{3036}{t} - 1.78)w_{sat}(1+0.448\; w_{sat})] +def wet_bulb_potential_temperature_from_dewpoint(*args, **kwargs): + return array.wet_bulb_potential_temperature_from_dewpoint(*args, **kwargs) - where: - * :math:`\Theta` is the :func:`potential_temperature` - * :math:`e_{sat}` is the :func:`saturation_vapor_pressure` - * :math:`q_{sat}` is the :func:`saturation_specific_humidity` - * :math:`w_{sat}` is the :func:`saturation_mixing_ratio` - * :math:`L_{v}` is the specific latent heat of vaporization (see :data:`earthkit.meteo.constants.Lv`) - * :math:`c_{pd}` is the specific heat of dry air on constant pressure - (see :data:`earthkit.meteo.constants.c_pd`) - - """ - return _EptComp.make(method).compute_ept_sat(t, p) - - -def temperature_on_moist_adiabat(ept, p, ept_method="ifs", t_method="bisect"): - r"""Computes the temperature on a moist adiabat (pseudoadiabat) - - Parameters - ---------- - ept: number or ndarray - Equivalent potential temperature defining the moist adiabat (K) - p: number or ndarray - Pressure on the moist adiabat (Pa) - ept_method: str, optional - Specifies the computation method that was used to compute ``ept``. The possible - values are: "ifs", "bolton35", "bolton39". - (See :func:`ept_from_dewpoint` for details.) - t_method: str, optional - Specifies the iteration method along the moist adiabat to find the temperature - for the given ``p`` pressure. The possible values are as follows: - - * "bisect": a bisection method is used as defined in [Stipanuk1973]_ - * "newton": Newtons's method is used as defined by Eq (2.6) in [DaviesJones2008]_. - For extremely hot and humid conditions (``ept`` > 800 K) depending on - ``ept_method`` the computation might not be carried out - and np.nan will be returned. - - - Returns - ------- - number or ndarray - Temperature on the moist adiabat (K). For values where the computation cannot - be carried out np.nan is returned. - - """ - cm = _EptComp.make(ept_method) - if t_method == "bisect": - return cm.compute_t_on_ma_stipanuk(ept, p) - elif t_method == "newton": - return cm.compute_t_on_ma_davies(ept, p) - else: - raise ValueError( - f"temperature_on_moist_adiabat: invalid t_method={t_method} specified!" - ) - - -def wet_bulb_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="bisect"): - r"""Computes the pseudo adiabatic wet bulb temperature from dewpoint. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - td: number or ndarray - Dewpoint (K) - p: number or ndarray - Pressure (Pa) - ept_method: str, optional - Specifies the computation method for the equivalent potential temperature. - The possible values are: "ifs", "bolton35", "bolton39". - (See :func:`ept_from_dewpoint` for details.) - t_method: str, optional - Specifies the method to find the temperature along the moist adiabat defined - by the equivalent potential temperature. The possible values are as follows: - - * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used - * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used - - Returns - ------- - number or ndarray - Wet bulb temperature (K) - - - The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): - - * first the equivalent potential temperature is computed with the given - ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. - * then the wet bulb potential temperature is determined as the temperature at - pressure ``p`` on the moist adiabat with the given ``t_method``. - - """ - ept = ept_from_dewpoint(t, td, p, method=ept_method) - return temperature_on_moist_adiabat( - ept, p, ept_method=ept_method, t_method=t_method - ) - - -def wet_bulb_temperature_from_specific_humidity( - t, q, p, ept_method="ifs", t_method="bisect" -): - r"""Computes the pseudo adiabatic wet bulb temperature from specific humidity. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - q: number or ndarray - Specific humidity (kg/kg) - p: number or ndarray - Pressure (Pa) - ept_method: str, optional - Specifies the computation method for the equivalent potential temperature. - The possible values are: "ifs", "bolton35", "bolton39". - (See :func:`ept_from_dewpoint` for details.) - t_method: str, optional - Specifies the method to find the temperature along the moist adiabat - defined by the equivalent potential temperature. The possible values are - as follows: - - * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used - * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used - - Returns - ------- - number or ndarray - Wet bulb temperature (K) - - - The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): - - * first the equivalent potential temperature is computed with the given - ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. - * then the wet bulb potential temperature is determined as the temperature at - pressure ``p`` on the moist adiabat with the given ``t_method``. - - """ - ept = ept_from_specific_humidity(t, q, p, method=ept_method) - return temperature_on_moist_adiabat( - ept, p, ept_method=ept_method, t_method=t_method - ) - - -def wet_bulb_potential_temperature_from_dewpoint( - t, td, p, ept_method="ifs", t_method="direct" -): - r"""Computes the pseudo adiabatic wet bulb potential temperature from dewpoint. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - td: number or ndarray - Dewpoint (K) - p: number or ndarray - Pressure (Pa) - ept_method: str, optional - Specifies the computation method for the equivalent potential temperature. - The possible values are: "ifs", "bolton35", "bolton39". - (See :func:`ept_from_dewpoint` for details.) - t_method: str, optional - Specifies the method to find the temperature along the moist adiabat defined - by the equivalent potential temperature. The possible values are as follows: - - * "direct": the rational formula defined by Eq (3.8) in [DaviesJones2008]_ is used - * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used - * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used - - Returns - ------- - number or ndarray - Wet bulb potential temperature (K) - - - The computation is based on Normand's rule [Wallace2006]_ (Chapter 3.5.6): - - * first the equivalent potential temperature is computed with the given - ``ept_method`` (using :func:`ept_from_dewpoint`). This defines the moist adiabat. - * then the wet bulb potential temperature is determined as the temperature at - pressure :math:`10^{5}` Pa on the moist adiabat with the given ``t_method``. - - """ - ept = ept_from_dewpoint(t, td, p, method=ept_method) - if t_method == "direct": - return _EptComp.make(ept_method).compute_wbpt(ept) - else: - return temperature_on_moist_adiabat( - ept, constants.p0, ept_method=ept_method, t_method=t_method - ) - - -def wet_bulb_potential_temperature_from_specific_humidity( - t, q, p, ept_method="ifs", t_method="direct" -): - r"""Computes the pseudo adiabatic wet bulb potential temperature from specific humidity. - - Parameters - ---------- - t: number or ndarray - Temperature (K) - q: number or ndarray - Specific humidity (kg/kg) - p: number or ndarray - Pressure (Pa) - ept_method: str, optional - Specifies the computation method for the equivalent potential temperature. - The possible values are: "ifs", "bolton35", "bolton39". - (See :func:`ept_from_dewpoint` for details.) - t_method: str, optional - Specifies the method to find the temperature along the moist adiabat - defined by the equivalent potential temperature. The possible values are as follows: - - * "direct": the rational formula defined by Eq (3.8) in [DaviesJones2008]_ is used - * "bisect": :func:`temperature_on_moist_adiabat` with ``t_method`` = "bisect" is used - * "newton": :func:`temperature_on_moist_adiabat` with ``t_method`` = "newton" is used - - Returns - ------- - number or ndarray - Wet bulb potential temperature (K) - - - The computations are the same as in - :func:`wet_bulb_potential_temperature_from_dewpoint` - (the dewpoint is computed from q with :func:`dewpoint_from_specific_humidity`). - - """ - ept = ept_from_specific_humidity(t, q, p, method=ept_method) - if t_method == "direct": - return _EptComp.make(ept_method).compute_wbpt(ept) - else: - return temperature_on_moist_adiabat( - ept, constants.p0, ept_method=ept_method, t_method=t_method - ) +def wet_bulb_potential_temperature_from_specific_humidity(*args, **kwargs): + return array.wet_bulb_potential_temperature_from_specific_humidity(*args, **kwargs) diff --git a/earthkit/meteo/wind/__init__.py b/earthkit/meteo/wind/__init__.py index b2cd29e..5aa9729 100644 --- a/earthkit/meteo/wind/__init__.py +++ b/earthkit/meteo/wind/__init__.py @@ -1 +1,18 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Wind related functions. + +The API is split into two levels. The low level functions are in the ``array`` submodule and they +can be used to operate on numpy arrays. The high level functions are still to be developed and +planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. +""" + from .wind import * # noqa diff --git a/earthkit/meteo/wind/array/__init__.py b/earthkit/meteo/wind/array/__init__.py new file mode 100644 index 0000000..1d61d5a --- /dev/null +++ b/earthkit/meteo/wind/array/__init__.py @@ -0,0 +1,14 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +""" +Wind related functions operating on numpy arrays. +""" + +from .wind import * # noqa diff --git a/earthkit/meteo/wind/array/wind.py b/earthkit/meteo/wind/array/wind.py new file mode 100644 index 0000000..f611ca9 --- /dev/null +++ b/earthkit/meteo/wind/array/wind.py @@ -0,0 +1,295 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + +from earthkit.meteo import constants + + +def speed(u, v): + r"""Computes the wind speed/vector magnitude. + + Parameters + ---------- + u: number or ndarray + u wind/x vector component + v: number or ndarray + v wind/y vector component (same units as ``u``) + + Returns + ------- + number or ndarray + Wind speed/magnitude (same units as ``u`` and ``v``) + + """ + return np.hypot(u, v) + + +def _direction_meteo(u, v): + minus_pi2 = -np.pi / 2.0 + d = np.arctan2(v, u) + d = np.asarray(d) + m = d <= minus_pi2 + d[m] = (minus_pi2 - d[m]) * constants.degree + m = ~m + d[m] = (1.5 * np.pi - d[m]) * constants.degree + return d + + +def _direction_polar(u, v, to_positive): + d = np.arctan2(v, u) * constants.degree + if to_positive: + d = np.asarray(d) + m = d < 0 + d[m] = 360.0 + d[m] + return d + + +def direction(u, v, convention="meteo", to_positive=True): + r"""Computes the direction/angle of a vector quantity. + + Parameters + ---------- + u: number or ndarray + u wind/x vector component + v: number or ndarray + v wind/y vector component (same units as ``u``) + convention: str, optional + Specifies how the direction/angle is interpreted. The possible values are as follows: + + * "meteo": the direction is the meteorological wind direction (see below for explanation) + * "polar": the direction is measured anti-clockwise from the x axis (East/right) to the vector + + positive: bool, optional + If it is True the resulting values are mapped into the [0, 360] range when + ``convention`` is "polar". Otherwise they lie in the [-180, 180] range. + + + Returns + ------- + number or ndarray + Direction/angle (degrees) + + + The meteorological wind direction is the direction from which the wind is + blowing. Wind direction increases clockwise such that a northerly wind is 0°, an easterly + wind is 90°, a southerly wind is 180°, and a westerly wind is 270°. The figure below illustrates + how it is related to the actual orientation of the wind vector: + + .. image:: /_static/wind_direction.png + :width: 400px + + """ + if convention == "meteo": + return _direction_meteo(u, v) + elif convention == "polar": + return _direction_polar(u, v, to_positive) + else: + raise ValueError(f"direction(): invalid convention={convention}!") + + +def xy_to_polar(x, y, convention="meteo"): + r"""Converts wind/vector data from xy representation to polar representation. + + Parameters + ---------- + x: number or ndarray + u wind/x vector component + y: number or ndarray + v wind/y vector component (same units as ``u``) + convention: str + Specifies how the direction/angle component of the target polar coordinate + system is interpreted. The possible values are as follows: + + * "meteo": the direction is the meteorological wind direction (see :func:`direction` for explanation) + * "polar": the direction is measured anti-clockwise from the x axis (East/right) to the vector + + + Returns + ------- + number or ndarray + Magnitude (same units as ``u``) + number or ndarray + Direction (degrees) + + + In the target xy representation the x axis points East while the y axis points North. + + """ + return speed(x, y), direction(x, y, convention=convention) + + +def _polar_to_xy_meteo(magnitude, direction): + a = (270.0 - direction) * constants.radian + return magnitude * np.cos(a), magnitude * np.sin(a) + + +def _polar_to_xy_polar(magnitude, direction): + a = direction * constants.radian + return magnitude * np.cos(a), magnitude * np.sin(a) + + +def polar_to_xy(magnitude, direction, convention="meteo"): + r"""Converts wind/vector data from polar representation to xy representation. + + Parameters + ---------- + magnitude: number or ndarray + Speed/magnitude of the vector + direction: number or ndarray + Direction of the vector (degrees) + convention: str + Specifies how ``direction`` is interpreted. The possible values are as follows: + + * "meteo": ``direction`` is the meteorological wind direction + (see :func:`direction` for explanation) + * "polar": ``direction`` is the angle measured anti-clockwise from the x axis + (East/right) to the vector + + Returns + ------- + number or ndarray + X vector component (same units as ``magnitude``) + number or ndarray + Y vector component (same units as ``magnitude``) + + + In the target xy representation the x axis points East while the y axis points North. + + """ + if convention == "meteo": + return _polar_to_xy_meteo(magnitude, direction) + elif convention == "polar": + return _polar_to_xy_polar(magnitude, direction) + else: + raise ValueError(f"polar_to_xy(): invalid convention={convention}!") + + +def w_from_omega(omega, t, p): + r"""Computes the hydrostatic vertical velocity from pressure velocity, temperature and pressure. + + Parameters + ---------- + omega : number or ndarray + Hydrostatic pressure velocity (Pa/s) + t : number or ndarray + Temperature (K) + p : number or ndarray + Pressure (Pa) + + Returns + ------- + number or ndarray + Hydrostatic vertical velocity (m/s) + + + The computation is based on the following hydrostatic formula: + + .. math:: + + w = - \frac{\omega\; t R_{d}}{p g} + + where + + * :math:`R_{d}` is the specific gas constant for dry air (see :data:`earthkit.meteo.constants.Rd`). + * :math:`g` is the gravitational acceleration (see :data:`earthkit.meteo.constants.g`) + + """ + return (-constants.Rd / constants.g) * (omega * t / p) + + +def coriolis(lat): + r"""Computes the Coriolis parameter. + + Parameters + ---------- + lat : number or ndarray + Latitude (degrees) + + Returns + ------- + number or ndarray + The Coriolis parameter (:math:`s^{-1}`) + + + The Coriolis parameter is defined by the following formula: + + .. math:: + + f = 2 \Omega sin(\phi) + + where :math:`\Omega` is the rotation rate of Earth + (see :data:`earthkit.meteo.constants.omega`) and :math:`\phi` is the latitude. + + """ + return 2 * constants.omega * np.sin(lat * constants.radian) + + +def windrose(speed, direction, sectors=16, speed_bins=[], percent=True): + """Generate windrose data. + + Parameters + ---------- + speed : number or ndarray + Speed + direction : number or ndarray + Meteorological wind direction (degrees). See :func:`direction` for details. + Values must be between 0 and 360. + sectors: number + Number of sectors the 360 degrees direction range is split into. See below for details. + speed_bin: list or ndarray + Speed bins + percent: bool + If False, returns the number of valid samples in each bin. If True, returns + the percentage of the number of samples in each bin with respect to the total + number of valid samples. + + + Returns + ------- + 2d-ndarray + The bi-dimensional histogram of ``speed`` and ``direction``. Values in + ``speed`` are histogrammed along the first dimension and values in ``direction`` + are histogrammed along the second dimension. + + ndarray + The direction bins (i.e. the sectors) (degrees) + + + The sectors do not start at 0 degrees (North) but are shifted by half a sector size. + E.g. if ``sectors`` is 4 the sectors are defined as: + + .. image:: /_static/wind_sector.png + :width: 350px + + """ + if len(speed_bins) < 2: + raise ValueError("windrose(): speed_bins must have at least 2 elements!") + + sectors = int(sectors) + if sectors < 1: + raise ValueError("windrose(): sectors must be greater than 1!") + + speed = np.atleast_1d(speed) + direction = np.atleast_1d(direction) + dir_step = 360.0 / sectors + dir_bins = np.linspace( + int(-dir_step / 2), int(360 + dir_step / 2), int(360 / dir_step) + 2 + ) + + res = np.histogram2d(speed, direction, bins=[speed_bins, dir_bins], density=False)[ + 0 + ] + + # unify the north bins + res[:, 0] = res[:, 0] + res[:, -1] + res = res[:, :-1] + dir_bins = dir_bins[:-1] + + return ((res * 100.0 / res.sum()) if percent else res), dir_bins diff --git a/earthkit/meteo/wind/wind.py b/earthkit/meteo/wind/wind.py index 1f55d87..3a28ae5 100644 --- a/earthkit/meteo/wind/wind.py +++ b/earthkit/meteo/wind/wind.py @@ -1,295 +1,38 @@ -# (C) Copyright 2021- ECMWF. +# +# # (C) Copyright 2021 ECMWF. +# # +# # This software is licensed under the terms of the Apache Licence Version 2.0 +# # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# # In applying this licence, ECMWF does not waive the privileges and immunities +# # granted to it by virtue of its status as an intergovernmental organisation +# # nor does it submit to any jurisdiction. +# ## +from . import array -# This software is licensed under the terms of the Apache Licence Version 2.0 -# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by -# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. +def speed(*args, **kwargs): + return array.speed(*args, **kwargs) -import numpy as np +def direction(*args, **kwargs): + return array.direction(*args, **kwargs) -from earthkit.meteo import constants +def xy_to_polar(*args, **kwargs): + return array.xy_to_polar(*args, **kwargs) -def speed(u, v): - r"""Computes the wind speed/vector magnitude. - Parameters - ---------- - u: number or ndarray - u wind/x vector component - v: number or ndarray - v wind/y vector component (same units as ``u``) +def polar_to_xy(*args, **kwargs): + return array.polar_to_xy(*args, **kwargs) - Returns - ------- - number or ndarray - Wind speed/magnitude (same units as ``u`` and ``v``) - """ - return np.hypot(u, v) +def w_from_omega(*args, **kwargs): + return array.w_from_omega(*args, **kwargs) -def _direction_meteo(u, v): - minus_pi2 = -np.pi / 2.0 - d = np.arctan2(v, u) - d = np.asarray(d) - m = d <= minus_pi2 - d[m] = (minus_pi2 - d[m]) * constants.degree - m = ~m - d[m] = (1.5 * np.pi - d[m]) * constants.degree - return d +def coriolis(*args, **kwargs): + return array.coriolis(*args, **kwargs) -def _direction_polar(u, v, to_positive): - d = np.arctan2(v, u) * constants.degree - if to_positive: - d = np.asarray(d) - m = d < 0 - d[m] = 360.0 + d[m] - return d - - -def direction(u, v, convention="meteo", to_positive=True): - r"""Computes the direction/angle of a vector quantity. - - Parameters - ---------- - u: number or ndarray - u wind/x vector component - v: number or ndarray - v wind/y vector component (same units as ``u``) - convention: str, optional - Specifies how the direction/angle is interpreted. The possible values are as follows: - - * "meteo": the direction is the meteorological wind direction (see below for explanation) - * "polar": the direction is measured anti-clockwise from the x axis (East/right) to the vector - - positive: bool, optional - If it is True the resulting values are mapped into the [0, 360] range when - ``convention`` is "polar". Otherwise they lie in the [-180, 180] range. - - - Returns - ------- - number or ndarray - Direction/angle (degrees) - - - The meteorological wind direction is the direction from which the wind is - blowing. Wind direction increases clockwise such that a northerly wind is 0°, an easterly - wind is 90°, a southerly wind is 180°, and a westerly wind is 270°. The figure below illustrates - how it is related to the actual orientation of the wind vector: - - .. image:: /_static/wind_direction.png - :width: 400px - - """ - if convention == "meteo": - return _direction_meteo(u, v) - elif convention == "polar": - return _direction_polar(u, v, to_positive) - else: - raise ValueError(f"direction(): invalid convention={convention}!") - - -def xy_to_polar(x, y, convention="meteo"): - r"""Converts wind/vector data from xy representation to polar representation. - - Parameters - ---------- - x: number or ndarray - u wind/x vector component - y: number or ndarray - v wind/y vector component (same units as ``u``) - convention: str - Specifies how the direction/angle component of the target polar coordinate - system is interpreted. The possible values are as follows: - - * "meteo": the direction is the meteorological wind direction (see :func:`direction` for explanation) - * "polar": the direction is measured anti-clockwise from the x axis (East/right) to the vector - - - Returns - ------- - number or ndarray - Magnitude (same units as ``u``) - number or ndarray - Direction (degrees) - - - In the target xy representation the x axis points East while the y axis points North. - - """ - return speed(x, y), direction(x, y, convention=convention) - - -def _polar_to_xy_meteo(magnitude, direction): - a = (270.0 - direction) * constants.radian - return magnitude * np.cos(a), magnitude * np.sin(a) - - -def _polar_to_xy_polar(magnitude, direction): - a = direction * constants.radian - return magnitude * np.cos(a), magnitude * np.sin(a) - - -def polar_to_xy(magnitude, direction, convention="meteo"): - r"""Converts wind/vector data from polar representation to xy representation. - - Parameters - ---------- - magnitude: number or ndarray - Speed/magnitude of the vector - direction: number or ndarray - Direction of the vector (degrees) - convention: str - Specifies how ``direction`` is interpreted. The possible values are as follows: - - * "meteo": ``direction`` is the meteorological wind direction - (see :func:`direction` for explanation) - * "polar": ``direction`` is the angle measured anti-clockwise from the x axis - (East/right) to the vector - - Returns - ------- - number or ndarray - X vector component (same units as ``magnitude``) - number or ndarray - Y vector component (same units as ``magnitude``) - - - In the target xy representation the x axis points East while the y axis points North. - - """ - if convention == "meteo": - return _polar_to_xy_meteo(magnitude, direction) - elif convention == "polar": - return _polar_to_xy_polar(magnitude, direction) - else: - raise ValueError(f"polar_to_xy(): invalid convention={convention}!") - - -def w_from_omega(omega, t, p): - r"""Computes the hydrostatic vertical velocity from pressure velocity, temperature and pressure. - - Parameters - ---------- - omega : number or ndarray - Hydrostatic pressure velocity (Pa/s) - t : number or ndarray - Temperature (K) - p : number or ndarray - Pressure (Pa) - - Returns - ------- - number or ndarray - Hydrostatic vertical velocity (m/s) - - - The computation is based on the following hydrostatic formula: - - .. math:: - - w = - \frac{\omega\; t R_{d}}{p g} - - where - - * :math:`R_{d}` is the specific gas constant for dry air (see :data:`earthkit.meteo.constants.Rd`). - * :math:`g` is the gravitational acceleration (see :data:`earthkit.meteo.constants.g`) - - """ - return (-constants.Rd / constants.g) * (omega * t / p) - - -def coriolis(lat): - r"""Computes the Coriolis parameter. - - Parameters - ---------- - lat : number or ndarray - Latitude (degrees) - - Returns - ------- - number or ndarray - The Coriolis parameter (:math:`s^{-1}`) - - - The Coriolis parameter is defined by the following formula: - - .. math:: - - f = 2 \Omega sin(\phi) - - where :math:`\Omega` is the rotation rate of Earth - (see :data:`earthkit.meteo.constants.omega`) and :math:`\phi` is the latitude. - - """ - return 2 * constants.omega * np.sin(lat * constants.radian) - - -def windrose(speed, direction, sectors=16, speed_bins=[], percent=True): - """Generate windrose data. - - Parameters - ---------- - speed : number or ndarray - Speed - direction : number or ndarray - Meteorological wind direction (degrees). See :func:`direction` for details. - Values must be between 0 and 360. - sectors: number - Number of sectors the 360 degrees direction range is split into. See below for details. - speed_bin: list or ndarray - Speed bins - percent: bool - If False, returns the number of valid samples in each bin. If True, returns - the percentage of the number of samples in each bin with respect to the total - number of valid samples. - - - Returns - ------- - 2d-ndarray - The bi-dimensional histogram of ``speed`` and ``direction``. Values in - ``speed`` are histogrammed along the first dimension and values in ``direction`` - are histogrammed along the second dimension. - - ndarray - The direction bins (i.e. the sectors) (degrees) - - - The sectors do not start at 0 degrees (North) but are shifted by half a sector size. - E.g. if ``sectors`` is 4 the sectors are defined as: - - .. image:: /_static/wind_sector.png - :width: 350px - - """ - if len(speed_bins) < 2: - raise ValueError("windrose(): speed_bins must have at least 2 elements!") - - sectors = int(sectors) - if sectors < 1: - raise ValueError("windrose(): sectors must be greater than 1!") - - speed = np.atleast_1d(speed) - direction = np.atleast_1d(direction) - dir_step = 360.0 / sectors - dir_bins = np.linspace( - int(-dir_step / 2), int(360 + dir_step / 2), int(360 / dir_step) + 2 - ) - - res = np.histogram2d(speed, direction, bins=[speed_bins, dir_bins], density=False)[ - 0 - ] - - # unify the north bins - res[:, 0] = res[:, 0] + res[:, -1] - res = res[:, :-1] - dir_bins = dir_bins[:-1] - - return ((res * 100.0 / res.sum()) if percent else res), dir_bins +def windrose(*args, **kwargs): + return array.windrose(*args, **kwargs) diff --git a/tests/test_extreme.py b/tests/extreme/test_extreme.py similarity index 94% rename from tests/test_extreme.py rename to tests/extreme/test_extreme.py index 5f81b2d..8412fae 100644 --- a/tests/test_extreme.py +++ b/tests/extreme/test_extreme.py @@ -1,10 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# import numpy as np @@ -504,7 +505,7 @@ def test_efi(): - efi = extreme.efi(clim, ens) + efi = extreme.array.efi(clim, ens) assert np.isclose(efi[0], -0.18384250406420133) @@ -513,34 +514,34 @@ def test_efi_sorted(): # ensures the algorithm is the same if we sort the data or not ens_perc = np.sort(ens) - efi = extreme.efi(clim, ens_perc) + efi = extreme.array.efi(clim, ens_perc) assert np.isclose(efi[0], -0.18384250406420133) def test_efi_eps(): - efi = extreme.efi(clim, ens, eps=1e-4) + efi = extreme.array.efi(clim, ens, eps=1e-4) assert np.isclose(efi[0], -0.18384250406420133) def test_efi_eps1(): - efi = extreme.efi(clim_eps, ens_eps, eps=1e-4) + efi = extreme.array.efi(clim_eps, ens_eps, eps=1e-4) # fortran code result is 0.4604220986366272 assert np.isclose(efi[0], 0.46039347745967046) def test_efi_eps2(): - efi = extreme.efi(clim_eps2, ens_eps2, eps=1e-4) + efi = extreme.array.efi(clim_eps2, ens_eps2, eps=1e-4) # fortran code result is assert np.isclose(efi[0], 0.6330071575726789) def test_sot(): - sot_upper = extreme.sot(clim, ens, 90) - sot_lower = extreme.sot(clim, ens, 10) + sot_upper = extreme.array.sot(clim, ens, 90) + sot_lower = extreme.array.sot(clim, ens, 10) print(sot_upper) print(sot_lower) @@ -555,7 +556,7 @@ def test_sot_missing(): qc = np.array([1.1, 1.0, 1.0, 1.00001]) qf = np.array([1.5, 1.2, 1.0, 0.9]) - sot = extreme.sot_func(qc_tail, qc, qf, eps=1e-4) + sot = extreme.array.sot_func(qc_tail, qc, qf, eps=1e-4) print(sot) print(np.isnan(sot)) @@ -571,7 +572,7 @@ def test_sot_bounds(): qc = np.array([1.1, 1.1]) qf = np.array([15, -15.0]) - sot = extreme.sot_func(qc_tail, qc, qf) + sot = extreme.array.sot_func(qc_tail, qc, qf) print(sot) assert np.allclose(sot, [-10, 10]) @@ -583,17 +584,17 @@ def test_sot_eps(): qc = np.array([0.1]) qf = np.array([0.2]) - sot = extreme.sot_func(qc_tail, qc, qf) + sot = extreme.array.sot_func(qc_tail, qc, qf) print(sot) assert np.isclose(sot[0], -3.0) - sot = extreme.sot_func(qc_tail, qc, qf, eps=0.15) + sot = extreme.array.sot_func(qc_tail, qc, qf, eps=0.15) print(sot) assert np.isnan(sot[0]) def test_sot_perc(): - sot = extreme.sot(clim_eps2, ens_eps2, 90, eps=1e4) + sot = extreme.array.sot(clim_eps2, ens_eps2, 90, eps=1e4) assert np.isnan(sot[0]) @@ -606,7 +607,7 @@ def test_efi_nan(): print(clim_nan) print(ens_nan) - efi = extreme.efi(clim_nan, ens_nan) + efi = extreme.array.efi(clim_nan, ens_nan) print(efi) assert np.isnan(efi[0]) @@ -617,21 +618,21 @@ def test_sot_nan(): qc = np.array([0.1]) qf = np.array([np.nan]) - sot = extreme.sot_func(qc_tail, qc, qf) + sot = extreme.array.sot_func(qc_tail, qc, qf) assert np.isnan(sot[0]) qc_tail = np.array([0.05]) qc = np.array([np.nan]) qf = np.array([0.1]) - sot = extreme.sot_func(qc_tail, qc, qf) + sot = extreme.array.sot_func(qc_tail, qc, qf) assert np.isnan(sot[0]) qc_tail = np.array([np.nan]) qc = np.array([0.1]) qf = np.array([0.2]) - sot = extreme.sot_func(qc_tail, qc, qf) + sot = extreme.array.sot_func(qc_tail, qc, qf) assert np.isnan(sot[0]) @@ -801,5 +802,24 @@ def test_sot_nan(): def test_cpf(): + cpf = extreme.array.cpf(cpf_clim, cpf_ens, sort_clim=True) + np.testing.assert_allclose(cpf, cpf_val) + + +def test_highlevel_efi(): + efi = extreme.efi(clim, ens) + + assert np.isclose(efi[0], -0.18384250406420133) + + +def test_highlevel_sot(): + sot_upper = extreme.sot(clim, ens, 90) + sot_lower = extreme.sot(clim, ens, 10) + + assert np.isclose(sot_upper[0], -2.14617638) + assert np.isclose(sot_lower[0], -1.3086723) + + +def test_highlevel_cpf(): cpf = extreme.cpf(cpf_clim, cpf_ens, sort_clim=True) np.testing.assert_allclose(cpf, cpf_val) diff --git a/tests/test_score.py b/tests/score/test_score.py similarity index 99% rename from tests/test_score.py rename to tests/score/test_score.py index 7deb852..aa80c04 100644 --- a/tests/test_score.py +++ b/tests/score/test_score.py @@ -1,10 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# import numpy as np diff --git a/tests/solar/test_solar.py b/tests/solar/test_solar.py new file mode 100644 index 0000000..0959128 --- /dev/null +++ b/tests/solar/test_solar.py @@ -0,0 +1,27 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import datetime + +import numpy as np +import pytest + +from earthkit.meteo import solar + + +@pytest.mark.parametrize( + "date,expected_value", + [ + (datetime.datetime(2024, 4, 22), 112.0), + (datetime.datetime(2024, 4, 22, 12, 0, 0), 112.5), + ], +) +def test_julian_day(date, expected_value): + v = solar.julian_day(date) + assert np.isclose(v, expected_value) diff --git a/tests/test_stats.py b/tests/stats/test_stats.py similarity index 98% rename from tests/test_stats.py rename to tests/stats/test_stats.py index 0978dd9..7709313 100644 --- a/tests/test_stats.py +++ b/tests/stats/test_stats.py @@ -1,10 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# import numpy as np import pytest diff --git a/tests/test_geo.py b/tests/test_geo.py deleted file mode 100644 index 682978c..0000000 --- a/tests/test_geo.py +++ /dev/null @@ -1,171 +0,0 @@ -# (C) Copyright 2023 ECMWF. -# -# This software is licensed under the terms of the Apache Licence Version 2.0 -# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -# In applying this licence, ECMWF does not waive the privileges and immunities -# granted to it by virtue of its status as an intergovernmental organisation -# nor does it submit to any jurisdiction. - - -import numpy as np - -from earthkit.meteo import geo - -# np.set_printoptions(formatter={"float_kind": "{:.10f}".format}) - - -def test_distance(): - # single values - lat_ref = 0 - lon_ref = 0 - - lats = np.array([0.0, 0, 0, 0, 90, -90, 48, 48, -48, -48, np.nan]) - lons = np.array([0, 90, -90, 180, 0, 0, 20, -20, -20, 20, 1.0]) - v_ref = np.array( - [ - 0.0000000000, - 10007903.1103691217, - 10007903.1103691217, - 20015806.2207382433, - 10007903.1103691217, - 10007903.1103691217, - 5675597.9227914885, - 5675597.9227914885, - 5675597.9227914885, - 5675597.9227914885, - np.nan, - ] - ) - ds = geo.distance(lat_ref, lon_ref, lats, lons) - np.testing.assert_allclose(ds, v_ref) - - lat_ref = 90.0 - lon_ref = 0 - v_ref = np.array( - [ - 10007903.1103691217, - 10007903.1103691217, - 10007903.1103691217, - 10007903.1103691217, - 0.0000000000, - 20015806.2207382433, - 4670354.7848389242, - 4670354.7848389242, - 15345451.4358993191, - 15345451.4358993191, - np.nan, - ] - ) - ds = geo.distance(lat_ref, lon_ref, lats, lons) - np.testing.assert_allclose(ds, v_ref) - - lat_ref = -90.0 - lon_ref = 0 - v_ref = np.array( - [ - 10007903.1103691217, - 10007903.1103691217, - 10007903.1103691217, - 10007903.1103691217, - 20015806.2207382433, - 0.0000000000, - 15345451.4358993191, - 15345451.4358993191, - 4670354.7848389242, - 4670354.7848389242, - np.nan, - ] - ) - ds = geo.distance(lat_ref, lon_ref, lats, lons) - np.testing.assert_allclose(ds, v_ref) - - lat_ref = 48.0 - lon_ref = 20 - v_ref = np.array( - [ - 5675597.9227914885, - 8536770.5279479641, - 11479035.6927902810, - 14340208.2979467567, - 4670354.7848389242, - 15345451.4358993191, - 0.0000000000, - 2942265.1648423159, - 11351195.8455829788, - 10675096.6510603968, - np.nan, - ] - ) - ds = geo.distance(lat_ref, lon_ref, lats, lons) - np.testing.assert_allclose(ds, v_ref) - - # arrays - lat_ref = np.array([48.0, 90.0]) - lon_ref = np.array([20, 0.0]) - v_ref = np.array([5675597.9227914885, 10007903.1103691217]) - - -def test_bearing(): - # single ref - multiple other points - lat_ref = 46 - lon_ref = 20 - lat = np.array([50.0, 46, 40, 46, 46, 46, -40, 80, 46, 50, 50, 40, 40, np.nan]) - lon = np.array([20.0, 24, 20, 16, -80, 100, 20, 20, 20, 22, 18, 18, 22, 22]) - v_ref = np.array( - [ - 0.0, - 90, - 180, - 270, - 270, - 90, - 180, - 0, - np.nan, - 19.0929472486, - 340.9070527514, - 193.0983348229, - 166.9016651771, - np.nan, - ] - ) - - b = geo.bearing(lat_ref, lon_ref, lat, lon) - np.testing.assert_allclose(b, v_ref) - - # multiple ref - multiple other points - lat_ref = [46, -14] - lon_ref = [20, -78] - lat = np.array([50.0, 46]) - lon = np.array([22.0, 24]) - v_ref = np.array([19.0929472486, 55.0620059697]) - b = geo.bearing(lat_ref, lon_ref, lat, lon) - np.testing.assert_allclose(b, v_ref) - - # multiple ref - single other points - lat_ref = [46, -14] - lon_ref = [20, -78] - lat = np.array([50.0]) - lon = np.array([22.0]) - v_ref = np.array([19.0929472486, 53.1446672968]) - b = geo.bearing(lat_ref, lon_ref, lat, lon) - np.testing.assert_allclose(b, v_ref) - - lat_ref = [46, -14] - lon_ref = [20, -78] - lat = np.array([46.0]) - lon = np.array([35.0]) - v_ref = np.array([90.0000000000, 54.7036374509]) - b = geo.bearing(lat_ref, lon_ref, lat, lon) - np.testing.assert_allclose(b, v_ref) - - # non-matching shapes - lat_ref = [46, -14] - lon_ref = [20, -78] - lat = np.array([46.0, 12, 12]) - lon = np.array([35.0, 22, 22]) - try: - b = geo.bearing(lat_ref, lon_ref, lat, lon) - assert False - except ValueError: - pass diff --git a/tests/test_version.py b/tests/test_version.py index 631a4b4..cd35996 100644 --- a/tests/test_version.py +++ b/tests/test_version.py @@ -1,10 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# import earthkit.meteo diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py new file mode 100644 index 0000000..4eda2fe --- /dev/null +++ b/tests/thermo/test_thermo.py @@ -0,0 +1,19 @@ +# (C) Copyright 2021 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np + +from earthkit.meteo import thermo + + +def test_high_level_celsius_to_kelvin(): + t = np.array([-10, 23.6]) + v = thermo.celsius_to_kelvin(t) + v_ref = np.array([263.16, 296.76]) + np.testing.assert_allclose(v, v_ref) diff --git a/tests/test_thermo.py b/tests/thermo/test_thermo_array.py similarity index 74% rename from tests/test_thermo.py rename to tests/thermo/test_thermo_array.py index 4cc5140..aa6a792 100644 --- a/tests/test_thermo.py +++ b/tests/thermo/test_thermo_array.py @@ -1,10 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. +# import os @@ -16,7 +17,7 @@ def data_file(name): - return os.path.join(os.path.dirname(__file__), "data", name) + return os.path.join(os.path.dirname(os.path.dirname(__file__)), "data", name) def save_test_reference(file_name, data): @@ -49,28 +50,28 @@ def __init__(self, file_name="t_hum_p_data.csv"): def test_celsius_to_kelvin(): t = np.array([-10, 23.6]) - v = thermo.celsius_to_kelvin(t) + v = thermo.array.celsius_to_kelvin(t) v_ref = np.array([263.16, 296.76]) np.testing.assert_allclose(v, v_ref) def test_kelvin_to_celsius(): t = np.array([263.16, 296.76]) - v = thermo.kelvin_to_celsius(t) + v = thermo.array.kelvin_to_celsius(t) v_ref = np.array([-10, 23.6]) np.testing.assert_allclose(v, v_ref) def test_mixing_ratio_from_specific_humidity(): q = np.array([0.008, 0.018]) - mr = thermo.mixing_ratio_from_specific_humidity(q) + mr = thermo.array.mixing_ratio_from_specific_humidity(q) v_ref = np.array([0.0080645161, 0.0183299389]) np.testing.assert_allclose(mr, v_ref, rtol=1e-7) def test_specific_humidity_from_mixing_ratio(): mr = np.array([0.0080645161, 0.0183299389]) - q = thermo.specific_humidity_from_mixing_ratio(mr) + q = thermo.array.specific_humidity_from_mixing_ratio(mr) v_ref = np.array([0.008, 0.018]) np.testing.assert_allclose(q, v_ref, rtol=1e-07) @@ -79,7 +80,7 @@ def test_vapour_pressure_from_specific_humidity(): q = np.array([0.008, 0.018]) p = np.array([700, 1000]) * 100 v_ref = np.array([895.992614, 2862.662152]) - vp = thermo.vapour_pressure_from_specific_humidity(q, p) + vp = thermo.array.vapour_pressure_from_specific_humidity(q, p) np.testing.assert_allclose(vp, v_ref) @@ -87,7 +88,7 @@ def test_vapour_pressure_from_mixing_ratio(): mr = np.array([0.0080645161, 0.0183299389]) p = np.array([700, 1000]) * 100 v_ref = np.array([895.992614, 2862.662152]) - vp = thermo.vapour_pressure_from_mixing_ratio(mr, p) + vp = thermo.array.vapour_pressure_from_mixing_ratio(mr, p) np.testing.assert_allclose(vp, v_ref) @@ -95,27 +96,27 @@ def test_specific_humidity_from_vapour_pressure(): vp = np.array([895.992614, 2862.662152, 10000]) p = np.array([700, 1000, 50]) * 100 v_ref = np.array([0.008, 0.018, np.nan]) - q = thermo.specific_humidity_from_vapour_pressure(vp, p) + q = thermo.array.specific_humidity_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) # single pressure value vp = np.array([895.992614, 2862.662152, 100000]) p = 700 * 100.0 v_ref = np.array([0.008, 0.0258354146, np.nan]) - q = thermo.specific_humidity_from_vapour_pressure(vp, p) + q = thermo.array.specific_humidity_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) # numbers vp = 895.992614 p = 700 * 100.0 v_ref = 0.008 - q = thermo.specific_humidity_from_vapour_pressure(vp, p) + q = thermo.array.specific_humidity_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) vp = 100000 p = 700 * 100.0 v_ref = np.nan - q = thermo.specific_humidity_from_vapour_pressure(vp, p) + q = thermo.array.specific_humidity_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) @@ -123,20 +124,20 @@ def test_mixing_ratio_from_vapour_pressure(): vp = np.array([895.992614, 2862.662152, 100000]) p = np.array([700, 1000, 50]) * 100 v_ref = np.array([0.0080645161, 0.0183299389, np.nan]) - mr = thermo.mixing_ratio_from_vapour_pressure(vp, p) + mr = thermo.array.mixing_ratio_from_vapour_pressure(vp, p) np.testing.assert_allclose(mr, v_ref, rtol=1e-07) # numbers vp = 895.992614 p = 700 * 100.0 v_ref = 0.0080645161 - q = thermo.mixing_ratio_from_vapour_pressure(vp, p) + q = thermo.array.mixing_ratio_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) vp = 100000.0 p = 50 * 100.0 v_ref = np.nan - q = thermo.mixing_ratio_from_vapour_pressure(vp, p) + q = thermo.array.mixing_ratio_from_vapour_pressure(vp, p) np.testing.assert_allclose(q, v_ref) @@ -144,9 +145,9 @@ def test_saturation_vapour_pressure(): ref_file = "sat_vp.csv" phases = ["mixed", "water", "ice"] - # o = {"t": thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 49))} + # o = {"t": thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 49))} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_vapour_pressure(o["t"], phase=phase) + # o[phase] = thermo.array.saturation_vapour_pressure(o["t"], phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -156,7 +157,7 @@ def test_saturation_vapour_pressure(): ) for phase in phases: - svp = thermo.saturation_vapour_pressure(d["t"], phase=phase) + svp = thermo.array.saturation_vapour_pressure(d["t"], phase=phase) np.testing.assert_allclose(svp, d[phase]) # scalar value @@ -167,7 +168,7 @@ def test_saturation_vapour_pressure(): "ice": 3.685208149695831139e02, } for phase, v in v_ref.items(): - svp = thermo.saturation_vapour_pressure(t, phase=phase) + svp = thermo.array.saturation_vapour_pressure(t, phase=phase) np.testing.assert_allclose(svp, v) @@ -175,14 +176,14 @@ def test_saturation_mixing_ratio(): ref_file = "sat_mr.csv" phases = ["mixed", "water", "ice"] - # t = thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) + # t = thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) # p = [1000, 950, 850, 700] # t_num = len(t) # t = np.repeat(t, repeats=len(p)) # p = np.array(p * t_num) * 100.0 # o = {"t": t, "p": p} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_mixing_ratio(t, p, phase=phase) + # o[phase] = thermo.array.saturation_mixing_ratio(t, p, phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -192,7 +193,7 @@ def test_saturation_mixing_ratio(): ) for phase in phases: - mr = thermo.saturation_mixing_ratio(d["t"], d["p"], phase=phase) + mr = thermo.array.saturation_mixing_ratio(d["t"], d["p"], phase=phase) np.testing.assert_allclose(mr, d[phase]) @@ -200,14 +201,14 @@ def test_saturation_specific_humidity(): ref_file = "sat_q.csv" phases = ["mixed", "water", "ice"] - # t = thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) + # t = thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) # p = [1000, 950, 850, 700] # t_num = len(t) # t = np.repeat(t, repeats=len(p)) # p = np.array(p * t_num) * 100.0 # o = {"t": t, "p": p} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_specific_humidity(t, p, phase=phase) + # o[phase] = thermo.array.saturation_specific_humidity(t, p, phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -217,7 +218,7 @@ def test_saturation_specific_humidity(): ) for phase in phases: - mr = thermo.saturation_specific_humidity(d["t"], d["p"], phase=phase) + mr = thermo.array.saturation_specific_humidity(d["t"], d["p"], phase=phase) np.testing.assert_allclose(mr, d[phase]) @@ -225,9 +226,9 @@ def test_saturation_vapour_pressure_slope(): ref_file = "sat_vp_slope.csv" phases = ["mixed", "water", "ice"] - # o = {"t": thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 49))} + # o = {"t": thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 49))} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_vapour_pressure_slope(o["t"], phase=phase) + # o[phase] = thermo.array.saturation_vapour_pressure_slope(o["t"], phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -237,7 +238,7 @@ def test_saturation_vapour_pressure_slope(): ) for phase in phases: - svp = thermo.saturation_vapour_pressure_slope(d["t"], phase=phase) + svp = thermo.array.saturation_vapour_pressure_slope(d["t"], phase=phase) np.testing.assert_allclose(svp, d[phase]) @@ -245,14 +246,14 @@ def test_saturation_mixing_ratio_slope(): ref_file = "sat_mr_slope.csv" phases = ["mixed", "water", "ice"] - # t = thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) + # t = thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) # p = [1000, 950, 850, 700] # t_num = len(t) # t = np.repeat(t, repeats=len(p)) # p = np.array(p * t_num) * 100.0 # o = {"t": t, "p": p} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_mixing_ratio_slope(t, p, phase=phase) + # o[phase] = thermo.array.saturation_mixing_ratio_slope(t, p, phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -262,13 +263,13 @@ def test_saturation_mixing_ratio_slope(): ) for phase in phases: - svp = thermo.saturation_mixing_ratio_slope(d["t"], d["p"], phase=phase) + svp = thermo.array.saturation_mixing_ratio_slope(d["t"], d["p"], phase=phase) np.testing.assert_allclose(svp, d[phase]) v_ref = np.array([np.nan]) - t = thermo.celsius_to_kelvin(np.array([200])) + t = thermo.array.celsius_to_kelvin(np.array([200])) p = np.array([1000]) - svp = thermo.saturation_mixing_ratio_slope(t, p, phase="mixed") + svp = thermo.array.saturation_mixing_ratio_slope(t, p, phase="mixed") np.testing.assert_allclose(svp, v_ref) # numbers @@ -276,7 +277,7 @@ def test_saturation_mixing_ratio_slope(): p = [1e5, 1e5] v_ref = [0.0005189819, np.nan] for i in range(len(t)): - svp = thermo.saturation_mixing_ratio_slope(t[i], p[i], phase="mixed") + svp = thermo.array.saturation_mixing_ratio_slope(t[i], p[i], phase="mixed") np.testing.assert_allclose(svp, v_ref[i]) @@ -284,14 +285,14 @@ def test_saturation_specific_humidity_slope(): ref_file = "sat_q_slope.csv" phases = ["mixed", "water", "ice"] - # t = thermo.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) + # t = thermo.array.celsius_to_kelvin(np.linspace(-40.0, 56.0, 25)) # p = [1000, 950, 850, 700] # t_num = len(t) # t = np.repeat(t, repeats=len(p)) # p = np.array(p * t_num) * 100.0 # o = {"t": t, "p": p} # for phase in ["mixed", "water", "ice"]: - # o[phase] = thermo.saturation_specific_humidity_slope(t, p, phase=phase) + # o[phase] = thermo.array.saturation_specific_humidity_slope(t, p, phase=phase) # save_test_reference(ref_file, o) d = np.genfromtxt( @@ -301,13 +302,15 @@ def test_saturation_specific_humidity_slope(): ) for phase in phases: - svp = thermo.saturation_specific_humidity_slope(d["t"], d["p"], phase=phase) + svp = thermo.array.saturation_specific_humidity_slope( + d["t"], d["p"], phase=phase + ) np.testing.assert_allclose(svp, d[phase]) v_ref = np.array([np.nan]) - t = thermo.celsius_to_kelvin(np.array([200])) + t = thermo.array.celsius_to_kelvin(np.array([200])) p = np.array([1000]) - svp = thermo.saturation_specific_humidity_slope(t, p, phase="mixed") + svp = thermo.array.saturation_specific_humidity_slope(t, p, phase="mixed") np.testing.assert_allclose(svp, v_ref) # numbers @@ -315,7 +318,7 @@ def test_saturation_specific_humidity_slope(): p = [1e5, 1e5] v_ref = [0.0005111349, np.nan] for i in range(len(t)): - svp = thermo.saturation_specific_humidity_slope(t[i], p[i], phase="mixed") + svp = thermo.array.saturation_specific_humidity_slope(t[i], p[i], phase="mixed") np.testing.assert_allclose(svp, v_ref[i]) @@ -327,13 +330,13 @@ def test_temperature_from_saturation_vapour_pressure(): names=True, ) - t = thermo.temperature_from_saturation_vapour_pressure(d["water"]) + t = thermo.array.temperature_from_saturation_vapour_pressure(d["water"]) np.testing.assert_allclose(t, d["t"]) def test_relative_humidity_from_dewpoint(): - t = thermo.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25])) - td = thermo.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3])) + t = thermo.array.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25])) + td = thermo.array.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3])) # reference was tested with an online relhum calculator at: # https://bmcnoldy.rsmas.miami.edu/Humidity.html v_ref = np.array( @@ -347,12 +350,12 @@ def test_relative_humidity_from_dewpoint(): 15.4779832381, ] ) - r = thermo.relative_humidity_from_dewpoint(t, td) + r = thermo.array.relative_humidity_from_dewpoint(t, td) np.testing.assert_allclose(r, v_ref, rtol=1e-05) def test_relative_humidity_from_specific_humidity(): - t = thermo.celsius_to_kelvin( + t = thermo.array.celsius_to_kelvin( np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]) ) p = np.array([300, 400, 500, 700, 850]) * 100.0 @@ -374,12 +377,12 @@ def test_relative_humidity_from_specific_humidity(): 73.41420898756694, ] ) - r = thermo.relative_humidity_from_specific_humidity(t, q, p) + r = thermo.array.relative_humidity_from_specific_humidity(t, q, p) np.testing.assert_allclose(r, v_ref) def test_specific_humidity_from_dewpoint(): - td = thermo.celsius_to_kelvin( + td = thermo.array.celsius_to_kelvin( np.array([21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916]) ) p = np.array([967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207]) * 100 @@ -393,12 +396,12 @@ def test_specific_humidity_from_dewpoint(): 0.0025150791, ] ) - q = thermo.specific_humidity_from_dewpoint(td, p) + q = thermo.array.specific_humidity_from_dewpoint(td, p) np.testing.assert_allclose(q, v_ref, rtol=1e-05) def test_specific_humidity_from_relative_humidity(): - t = thermo.celsius_to_kelvin( + t = thermo.array.celsius_to_kelvin( np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]) ) p = np.array([300, 400, 500, 700, 850]) * 100.0 @@ -420,12 +423,12 @@ def test_specific_humidity_from_relative_humidity(): 0.0114808182580539, ] ) - q = thermo.specific_humidity_from_relative_humidity(t, r, p) + q = thermo.array.specific_humidity_from_relative_humidity(t, r, p) np.testing.assert_allclose(q, v_ref) def test_dewpoint_from_relative_humidity(): - t = thermo.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25])) + t = thermo.array.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25])) r = np.array( [ 100.0000000000, @@ -437,10 +440,10 @@ def test_dewpoint_from_relative_humidity(): 15.4779832381, ] ) - v_ref = thermo.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3])) + v_ref = thermo.array.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3])) # reference was tested with an online relhum calculator at: # https://bmcnoldy.rsmas.miami.edu/Humidity.html - td = thermo.dewpoint_from_relative_humidity(t, r) + td = thermo.array.dewpoint_from_relative_humidity(t, r) np.testing.assert_allclose(td, v_ref) @@ -456,10 +459,10 @@ def test_dewpoint_from_specific_humidity(): 0.0025150791, ] ) - v_ref = thermo.celsius_to_kelvin( + v_ref = thermo.array.celsius_to_kelvin( np.array([21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916]) ) - td = thermo.dewpoint_from_specific_humidity(q, p) + td = thermo.array.dewpoint_from_specific_humidity(q, p) np.testing.assert_allclose(td, v_ref) @@ -468,7 +471,7 @@ def test_virtual_temperature(): # w = np.array([0.02, 0.03]) q = np.array([0.0196078431, 0.0291262136]) v_ref = np.array([289.8130240470, 298.5937453245]) - tv = thermo.virtual_temperature(t, q) + tv = thermo.array.virtual_temperature(t, q) np.testing.assert_allclose(tv, v_ref) @@ -478,7 +481,7 @@ def test_virtual_potential_temperature_temperature(): q = np.array([0.0196078431, 0.0291262136]) p = np.array([100300.0, 95000.0]) v_ref = np.array([289.5651110613, 303.0015650834]) - tv = thermo.virtual_potential_temperature(t, q, p) + tv = thermo.array.virtual_potential_temperature(t, q, p) np.testing.assert_allclose(tv, v_ref) @@ -486,7 +489,7 @@ def test_potential_temperature(): t = np.array([252.16, 298.16]) p = np.array([72350, 100500]) v_ref = np.array([276.588026, 297.735455]) - th = thermo.potential_temperature(t, p) + th = thermo.array.potential_temperature(t, p) np.testing.assert_allclose(th, v_ref) @@ -494,7 +497,7 @@ def test_temperature_from_potential_temperature(): p = np.array([72350, 100500]) th = np.array([276.588026, 297.735455]) v_ref = np.array([252.16, 298.16]) - t = thermo.temperature_from_potential_temperature(th, p) + t = thermo.array.temperature_from_potential_temperature(th, p) np.testing.assert_allclose(t, v_ref) @@ -503,17 +506,17 @@ def test_temperature_on_dry_adibat(): p_def = np.array([72350, 100500]) p = np.array([700, 500]) * 100 v_ref = np.array([249.792414, 244.246863]) - t = thermo.temperature_on_dry_adiabat(p, t_def, p_def) + t = thermo.array.temperature_on_dry_adiabat(p, t_def, p_def) np.testing.assert_allclose(t, v_ref) # cross checking - th1 = thermo.potential_temperature(t_def, p_def) - th2 = thermo.potential_temperature(t, p) + th1 = thermo.array.potential_temperature(t_def, p_def) + th2 = thermo.array.potential_temperature(t, p) np.testing.assert_allclose(th1, th2) # multiple values along a single adiabat v_ref = np.array([249.792414, 226.898581]) - t = thermo.temperature_on_dry_adiabat(p, t_def[0], p_def[0]) + t = thermo.array.temperature_on_dry_adiabat(p, t_def[0], p_def[0]) np.testing.assert_allclose(t, v_ref) @@ -522,17 +525,17 @@ def test_pressure_on_dry_adibat(): p_def = np.array([72350, 100500]) t = np.array([249.792414, 244.246863]) v_ref = np.array([700, 500]) * 100 - p = thermo.pressure_on_dry_adiabat(t, t_def, p_def) + p = thermo.array.pressure_on_dry_adiabat(t, t_def, p_def) np.testing.assert_allclose(p, v_ref) # cross checking - th1 = thermo.potential_temperature(t_def, p_def) - th2 = thermo.potential_temperature(t, p) + th1 = thermo.array.potential_temperature(t_def, p_def) + th2 = thermo.array.potential_temperature(t, p) np.testing.assert_allclose(th1, th2) # multiple values along a single adiabat v_ref = np.array([70000, 64709.699161]) - p = thermo.pressure_on_dry_adiabat(t, t_def[0], p_def[0]) + p = thermo.array.pressure_on_dry_adiabat(t, t_def[0], p_def[0]) np.testing.assert_allclose(p, v_ref) @@ -544,18 +547,18 @@ def test_lcl(): # davies t_ref = np.array([268.706024, 298.200936, 270.517934, 303.138144]) p_ref = np.array([79081.766347, 94862.350635, 57999.83367, 112654.210439]) - t_lcl = thermo.lcl_temperature(t, td, method="davies") + t_lcl = thermo.array.lcl_temperature(t, td, method="davies") np.testing.assert_allclose(t_lcl, t_ref) - t_lcl, p_lcl = thermo.lcl(t, td, p, method="davies") + t_lcl, p_lcl = thermo.array.lcl(t, td, p, method="davies") np.testing.assert_allclose(t_lcl, t_ref) np.testing.assert_allclose(p_lcl, p_ref) # bolton t_ref = np.array([268.683018, 298.182282, 270.531264, 303.199544]) p_ref = np.array([79058.068785, 94841.581142, 58009.838027, 112734.100243]) - t_lcl = thermo.lcl_temperature(t, td, method="bolton") + t_lcl = thermo.array.lcl_temperature(t, td, method="bolton") np.testing.assert_allclose(t_lcl, t_ref) - t_lcl, p_lcl = thermo.lcl(t, td, p, method="bolton") + t_lcl, p_lcl = thermo.array.lcl(t, td, p, method="bolton") np.testing.assert_allclose(t_lcl, t_ref) np.testing.assert_allclose(p_lcl, p_ref) @@ -568,8 +571,8 @@ def test_ept(): # o = {} # for m in methods: - # o[f"{m}_td"] = thermo.ept_from_dewpoint(data.t, data.td, data.p, method=m) - # o[f"{m}_q"] = thermo.ept_from_specific_humidity( + # o[f"{m}_td"] = thermo.array.ept_from_dewpoint(data.t, data.td, data.p, method=m) + # o[f"{m}_q"] = thermo.array.ept_from_specific_humidity( # data.t, data.q, data.p, method=m # ) # save_test_reference(ref_file, o) @@ -581,9 +584,9 @@ def test_ept(): ) for m in methods: - pt = thermo.ept_from_dewpoint(data.t, data.td, data.p, method=m) + pt = thermo.array.ept_from_dewpoint(data.t, data.td, data.p, method=m) np.testing.assert_allclose(pt, ref[m + "_td"], err_msg=f"method={m}") - pt = thermo.ept_from_specific_humidity(data.t, data.q, data.p, method=m) + pt = thermo.array.ept_from_specific_humidity(data.t, data.q, data.p, method=m) np.testing.assert_allclose(pt, ref[m + "_q"], err_msg=f"method={m}") @@ -595,7 +598,7 @@ def test_saturation_ept(): # o = {} # for m in methods: - # o[m] = thermo.saturation_ept(data.t, data.p, method=m) + # o[m] = thermo.array.saturation_ept(data.t, data.p, method=m) # save_test_reference(ref_file, o) ref = np.genfromtxt( @@ -605,7 +608,7 @@ def test_saturation_ept(): ) for m in methods: - pt = thermo.saturation_ept(data.t, data.p, method=m) + pt = thermo.array.saturation_ept(data.t, data.p, method=m) np.testing.assert_allclose(pt, ref[m], err_msg=f"method={m}") @@ -619,7 +622,7 @@ def test_temperature_on_moist_adiabat(): # o = {"ept": np.repeat(ept, repeats=len(p)), "p": np.array(p.tolist() * len(ept))} # for m_ept in ept_methods: # for m_t in t_methods: - # o[f"{m_ept}_{m_t}"] = thermo.temperature_on_moist_adiabat( + # o[f"{m_ept}_{m_t}"] = thermo.array.temperature_on_moist_adiabat( # o["ept"], o["p"], ept_method=m_ept, t_method=m_t # ) # save_test_reference(ref_file, o) @@ -632,7 +635,7 @@ def test_temperature_on_moist_adiabat(): for m_ept in ept_methods: for m_t in t_methods: - pt = thermo.thermo.temperature_on_moist_adiabat( + pt = thermo.array.temperature_on_moist_adiabat( ref["ept"], ref["p"], ept_method=m_ept, t_method=m_t ) np.testing.assert_allclose( @@ -650,10 +653,10 @@ def test_wet_bulb_temperature(): # o = {} # for m_ept in ept_methods: # for m_t in t_methods: - # o[f"{m_ept}_{m_t}_td"] = thermo.wet_bulb_temperature_from_dewpoint( + # o[f"{m_ept}_{m_t}_td"] = thermo.array.wet_bulb_temperature_from_dewpoint( # data.t, data.td, data.p, ept_method=m_ept, t_method=m_t # ) - # o[f"{m_ept}_{m_t}_q"] = thermo.wet_bulb_temperature_from_specific_humidity( + # o[f"{m_ept}_{m_t}_q"] = thermo.array.wet_bulb_temperature_from_specific_humidity( # data.t, data.q, data.p, ept_method=m_ept, t_method=m_t # ) # save_test_reference(ref_file, o) @@ -666,7 +669,7 @@ def test_wet_bulb_temperature(): for m_ept in ept_methods: for m_t in t_methods: - pt = thermo.wet_bulb_temperature_from_dewpoint( + pt = thermo.array.wet_bulb_temperature_from_dewpoint( data.t, data.td, data.p, ept_method=m_ept, t_method=m_t ) np.testing.assert_allclose( @@ -676,7 +679,7 @@ def test_wet_bulb_temperature(): atol=0, err_msg=f"method={m_ept}_{m_t}_td", ) - pt = thermo.wet_bulb_temperature_from_specific_humidity( + pt = thermo.array.wet_bulb_temperature_from_specific_humidity( data.t, data.q, data.p, ept_method=m_ept, t_method=m_t ) np.testing.assert_allclose( @@ -700,12 +703,12 @@ def test_wet_bulb_potential_temperature(): # for m_t in t_methods: # o[ # f"{m_ept}_{m_t}_td" - # ] = thermo.wet_bulb_potential_temperature_from_dewpoint( + # ] = thermo.array.wet_bulb_potential_temperature_from_dewpoint( # data.t, data.td, data.p, ept_method=m_ept, t_method=m_t # ) # o[ # f"{m_ept}_{m_t}_q" - # ] = thermo.wet_bulb_potential_temperature_from_specific_humidity( + # ] = thermo.array.wet_bulb_potential_temperature_from_specific_humidity( # data.t, data.q, data.p, ept_method=m_ept, t_method=m_t # ) # save_test_reference(ref_file, o) @@ -718,7 +721,7 @@ def test_wet_bulb_potential_temperature(): for m_ept in ept_methods: for m_t in t_methods: - pt = thermo.wet_bulb_potential_temperature_from_dewpoint( + pt = thermo.array.wet_bulb_potential_temperature_from_dewpoint( data.t, data.td, data.p, ept_method=m_ept, t_method=m_t ) np.testing.assert_allclose( @@ -728,7 +731,7 @@ def test_wet_bulb_potential_temperature(): atol=0, err_msg=f"method={m_ept}_{m_t}_td", ) - pt = thermo.wet_bulb_potential_temperature_from_specific_humidity( + pt = thermo.array.wet_bulb_potential_temperature_from_specific_humidity( data.t, data.q, data.p, ept_method=m_ept, t_method=m_t ) np.testing.assert_allclose( diff --git a/tests/test_wind.py b/tests/wind/test_wind.py similarity index 99% rename from tests/test_wind.py rename to tests/wind/test_wind.py index 0828c8f..09f1b1d 100644 --- a/tests/test_wind.py +++ b/tests/wind/test_wind.py @@ -1,11 +1,11 @@ -# (C) Copyright 2023 ECMWF. +# (C) Copyright 2021 ECMWF. # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. # In applying this licence, ECMWF does not waive the privileges and immunities # granted to it by virtue of its status as an intergovernmental organisation # nor does it submit to any jurisdiction. - +# import numpy as np