From ceea222c37ed254ed4f457c02ec24d54a118359e Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Wed, 30 Oct 2024 09:32:32 +0000 Subject: [PATCH 1/4] Add vertical submodule (#21) --- src/earthkit/meteo/vertical/__init__.py | 19 +++++++++++++++++++ src/earthkit/meteo/vertical/array/__init__.py | 14 ++++++++++++++ src/earthkit/meteo/vertical/array/vertical.py | 8 ++++++++ src/earthkit/meteo/vertical/vertical.py | 10 ++++++++++ 4 files changed, 51 insertions(+) create mode 100644 src/earthkit/meteo/vertical/__init__.py create mode 100644 src/earthkit/meteo/vertical/array/__init__.py create mode 100644 src/earthkit/meteo/vertical/array/vertical.py create mode 100644 src/earthkit/meteo/vertical/vertical.py diff --git a/src/earthkit/meteo/vertical/__init__.py b/src/earthkit/meteo/vertical/__init__.py new file mode 100644 index 0000000..5552f22 --- /dev/null +++ b/src/earthkit/meteo/vertical/__init__.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. +# + +""" +Vertical computation 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 .vertical import * # noqa diff --git a/src/earthkit/meteo/vertical/array/__init__.py b/src/earthkit/meteo/vertical/array/__init__.py new file mode 100644 index 0000000..9884ca5 --- /dev/null +++ b/src/earthkit/meteo/vertical/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. +# + +""" +Vertical computation functions operating on numpy arrays. +""" + +from .vertical import * # noqa diff --git a/src/earthkit/meteo/vertical/array/vertical.py b/src/earthkit/meteo/vertical/array/vertical.py new file mode 100644 index 0000000..975426c --- /dev/null +++ b/src/earthkit/meteo/vertical/array/vertical.py @@ -0,0 +1,8 @@ +# (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. +# diff --git a/src/earthkit/meteo/vertical/vertical.py b/src/earthkit/meteo/vertical/vertical.py new file mode 100644 index 0000000..5a1d5be --- /dev/null +++ b/src/earthkit/meteo/vertical/vertical.py @@ -0,0 +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 . import array From d6175a859dca3b110995068b133a87ea0ef38642 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Wed, 30 Oct 2024 11:03:06 +0000 Subject: [PATCH 2/4] Update pre-commit and ci configuration (#22) --- .github/workflows/ci.yml | 18 ++- .github/workflows/legacy-ci.yml | 14 +-- .github/workflows/python-pill-request.yml | 12 ++ .pre-commit-config.yaml | 56 +++++---- pyproject.toml | 58 +++++----- src/earthkit/meteo/extreme/array/cpf.py | 14 +-- src/earthkit/meteo/extreme/array/efi.py | 14 +-- src/earthkit/meteo/extreme/array/sot.py | 32 ++---- src/earthkit/meteo/score/array/crps.py | 4 +- src/earthkit/meteo/stats/array/quantiles.py | 8 +- src/earthkit/meteo/thermo/array/thermo.py | 120 +++++--------------- src/earthkit/meteo/wind/array/wind.py | 8 +- tests/score/test_score.py | 8 +- tests/solar/test_solar.py | 4 +- tests/stats/test_stats.py | 12 +- tests/thermo/test_thermo_array.py | 16 +-- tests/wind/test_wind.py | 16 +-- 17 files changed, 155 insertions(+), 259 deletions(-) create mode 100644 .github/workflows/python-pill-request.yml diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7783738..eb3a265 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -4,13 +4,19 @@ on: # Trigger the workflow on push to master or develop, except tag creation push: branches: - - 'main' - - 'develop' + - "main" + - "develop" tags-ignore: - - '**' + - "**" + paths-ignore: + - "docs/**" + - "README.md" # Trigger the workflow on pull request pull_request: + paths-ignore: + - "docs/**" + - "README.md" # Trigger the workflow manually workflow_dispatch: @@ -18,6 +24,9 @@ on: # Trigger after public PR approved for CI pull_request_target: types: [labeled] + paths-ignore: + - "docs/**" + - "README.md" jobs: # Run CI including downstream packages on self-hosted runners @@ -28,10 +37,9 @@ jobs: with: earthkit-meteo: ecmwf/earthkit-meteo@${{ github.event.pull_request.head.sha || github.sha }} codecov_upload: true - python_qa: true + python_qa: false secrets: inherit - # Build downstream packages on HPC downstream-ci-hpc: name: downstream-ci-hpc diff --git a/.github/workflows/legacy-ci.yml b/.github/workflows/legacy-ci.yml index b9619d8..77b015a 100644 --- a/.github/workflows/legacy-ci.yml +++ b/.github/workflows/legacy-ci.yml @@ -24,18 +24,6 @@ defaults: shell: bash -l {0} jobs: - pre-commit: - if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v3 - with: - ref: ${{ github.event.pull_request.head.sha || github.ref }} - - uses: actions/setup-python@v4 - with: - python-version: 3.x - - uses: pre-commit/action@v3.0.0 - documentation: if: ${{ !github.event.pull_request.head.repo.fork && github.event.action != 'labeled' || github.event.label.name == 'approved-for-ci' }} runs-on: ubuntu-latest @@ -48,7 +36,7 @@ jobs: with: ref: ${{ github.event.pull_request.head.sha || github.ref }} - name: Install Conda environment with Micromamba - uses: mamba-org/provision-with-micromamba@v12 + uses: mamba-org/setup-micromamba@v1 with: environment-file: tests/environment-unit-tests.yml environment-name: DEVELOP diff --git a/.github/workflows/python-pill-request.yml b/.github/workflows/python-pill-request.yml new file mode 100644 index 0000000..2c38063 --- /dev/null +++ b/.github/workflows/python-pill-request.yml @@ -0,0 +1,12 @@ +name: Code Quality checks for PRs + +on: + push: + pull_request: + types: [opened, synchronize, reopened] + +jobs: + quality: + uses: ecmwf-actions/reusable-workflows/.github/workflows/qa-precommit-run.yml@v2 + with: + skip-hooks: "no-commit-to-branch" diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5ffa512..260d26e 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,32 +1,48 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.4.0 + rev: v5.0.0 hooks: - - id: trailing-whitespace - - id: end-of-file-fixer + - id: trailing-whitespace # Trailing whitespace checker + - id: end-of-file-fixer # Ensure files end in a newline - id: check-json - - id: check-yaml + - id: check-yaml # Check YAML files for syntax errors only + args: [--unsafe, --allow-multiple-documents] - id: check-toml - # - id: check-added-large-files - - id: debug-statements + # - id: check-added-large-files + - id: debug-statements # Check for debugger imports and py37+ breakpoint() - id: mixed-line-ending + - id: no-commit-to-branch # Prevent committing to main / master + - id: check-merge-conflict # Check for files that contain merge conflict + exclude: /README\.rst$|^docs/.*\.rst$ - repo: https://github.com/PyCQA/isort - rev: 5.12.0 + rev: 5.13.2 hooks: - id: isort + args: + - -l 110 + - --force-single-line-imports + - --profile black - repo: https://github.com/psf/black - rev: 23.9.1 + rev: 24.8.0 hooks: - id: black + args: [--line-length=110] - repo: https://github.com/keewis/blackdoc rev: v0.3.8 hooks: - id: blackdoc additional_dependencies: [black==23.3.0] -- repo: https://github.com/PyCQA/flake8 - rev: 6.1.0 - hooks: - - id: flake8 + exclude: xr_engine_profile_rst\.py +- repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.6.9 + hooks: + - id: ruff + exclude: '(dev/.*|.*_)\.py$' + args: + - --line-length=110 + - --fix + - --exit-non-zero-on-fix + - --preview - repo: https://github.com/executablebooks/mdformat rev: 0.7.14 hooks: @@ -37,11 +53,11 @@ repos: hooks: - id: pretty-format-yaml args: [--autofix, --preserve-quotes] - - id: pretty-format-toml - args: [--autofix] -- repo: https://github.com/PyCQA/pydocstyle.git - rev: 6.1.1 - hooks: - - id: pydocstyle - additional_dependencies: [toml] - exclude: tests|docs +- repo: https://github.com/sphinx-contrib/sphinx-lint + rev: v1.0.0 + hooks: + - id: sphinx-lint +- repo: https://github.com/tox-dev/pyproject-fmt + rev: "2.2.4" + hooks: + - id: pyproject-fmt diff --git a/pyproject.toml b/pyproject.toml index 7279832..ad4cc62 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,58 +1,56 @@ [build-system] -requires = ["setuptools>=61", "setuptools-scm>=8.0"] +requires = [ "setuptools>=61", "setuptools-scm>=8" ] [project] +name = "earthkit-meteo" +description = "Meteorological computations" +readme = "README.md" +license = { text = "Apache License Version 2.0" } authors = [ - {name = "European Centre for Medium-Range Weather Forecasts (ECMWF)", email = "software.support@ecmwf.int"} + { name = "European Centre for Medium-Range Weather Forecasts (ECMWF)", email = "software.support@ecmwf.int" }, ] +requires-python = ">=3.8" + classifiers = [ "Development Status :: 2 - Pre-Alpha", "Intended Audience :: Developers", "License :: OSI Approved :: Apache Software License", - "Programming Language :: Python :: 3", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3 :: Only", "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", "Programming Language :: Python :: Implementation :: CPython", "Programming Language :: Python :: Implementation :: PyPy", - "Operating System :: OS Independent" ] +dynamic = [ "version" ] dependencies = [ - "numpy" + "numpy", ] -description = "Meteorological computations" -dynamic = ["version"] -license = {text = "Apache License Version 2.0"} -name = "earthkit-meteo" -readme = "README.md" -requires-python = ">= 3.8" - -[project.optional-dependencies] -test = [ +optional-dependencies.test = [ "pytest", - "pytest-cov" + "pytest-cov", ] +urls.Documentation = "https://earthkit-meteo.readthedocs.io/" +urls.Homepage = "https://github.com/ecmwf/earthkit-meteo/" +urls.Issues = "https://github.com/ecmwf/earthkit-meteo.issues" +urls.Repository = "https://github.com/ecmwf/earthkit-meteo/" -[project.urls] -Documentation = "https://earthkit-meteo.readthedocs.io/" -Homepage = "https://github.com/ecmwf/earthkit-meteo/" -Issues = "https://github.com/ecmwf/earthkit-meteo.issues" -Repository = "https://github.com/ecmwf/earthkit-meteo/" +[tool.setuptools.packages.find] +include = [ "earthkit.meteo" ] +where = [ "src/" ] -[tool.coverage.run] -branch = "true" +[tool.setuptools_scm] +version_file = "src/earthkit/meteo/_version.py" [tool.isort] profile = "black" +[tool.coverage.run] +branch = "true" + [tool.pydocstyle] -add_ignore = ["D1", "D200", "D205", "D400", "D401"] +add_ignore = [ "D1", "D200", "D205", "D400", "D401" ] convention = "numpy" - -[tool.setuptools.packages.find] -include = ["earthkit.meteo"] -where = ["src/"] - -[tool.setuptools_scm] -version_file = "src/earthkit/meteo/_version.py" diff --git a/src/earthkit/meteo/extreme/array/cpf.py b/src/earthkit/meteo/extreme/array/cpf.py index 16596de..1a737a9 100644 --- a/src/earthkit/meteo/extreme/array/cpf.py +++ b/src/earthkit/meteo/extreme/array/cpf.py @@ -64,10 +64,9 @@ def cpf(clim, ens, sort_clim=True, sort_ens=True): 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]) + 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) @@ -85,10 +84,9 @@ def cpf(clim, ens, sort_clim=True, sort_ens=True): 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]) + 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) diff --git a/src/earthkit/meteo/extreme/array/efi.py b/src/earthkit/meteo/extreme/array/efi.py index 455222b..6c08911 100644 --- a/src/earthkit/meteo/extreme/array/efi.py +++ b/src/earthkit/meteo/extreme/array/efi.py @@ -31,9 +31,7 @@ def efi(clim, ens, eps=-0.1): EFI values """ # locate missing values - missing_mask = np.logical_or( - np.sum(np.isnan(clim), axis=0), np.sum(np.isnan(ens), axis=0) - ) + 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 @@ -65,9 +63,7 @@ def efi(clim, ens, eps=-0.1): mask = clim[icl + 1, :] > eps dEFI = np.where( mask, - (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] - + acoef[icl] * dFdp[icl, :] - - proddiff[icl], + (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) @@ -77,11 +73,7 @@ def efi(clim, ens, eps=-0.1): efi /= efimax else: for icl in range(nclim - 1): - dEFI = ( - (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] - + acoef[icl] * dFdp[icl, :] - - proddiff[icl] - ) + dEFI = (2.0 * frac[icl, :] - 1.0) * acosdiff[icl] + acoef[icl] * dFdp[icl, :] - proddiff[icl] efi += dEFI efi *= 2.0 / np.pi ################################## diff --git a/src/earthkit/meteo/extreme/array/sot.py b/src/earthkit/meteo/extreme/array/sot.py index f556afa..8a7d933 100644 --- a/src/earthkit/meteo/extreme/array/sot.py +++ b/src/earthkit/meteo/extreme/array/sot.py @@ -39,9 +39,7 @@ def sot_func(qc_tail, qc, qf, eps=-1e-4, lower_bound=-10, upper_bound=10): 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 - ) + 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) @@ -78,20 +76,12 @@ def sot(clim, ens, perc, eps=-1e4): 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 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 - ) + "Climatology array should contain 101 percentiles, it has {} values".format(clim.shape) ) qc = clim[perc] @@ -136,20 +126,12 @@ def sot_unsorted(clim, ens, perc, eps=-1e4): 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 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 - ) + "Climatology array should contain 101 percentiles, it has {} values".format(clim.shape) ) if eps > 0: diff --git a/src/earthkit/meteo/score/array/crps.py b/src/earthkit/meteo/score/array/crps.py index 9595094..dcce923 100644 --- a/src/earthkit/meteo/score/array/crps.py +++ b/src/earthkit/meteo/score/array/crps.py @@ -48,9 +48,7 @@ def crps(x, y): 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 + 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 diff --git a/src/earthkit/meteo/stats/array/quantiles.py b/src/earthkit/meteo/stats/array/quantiles.py index 0f9109a..b590a1a 100644 --- a/src/earthkit/meteo/stats/array/quantiles.py +++ b/src/earthkit/meteo/stats/array/quantiles.py @@ -7,7 +7,9 @@ # nor does it submit to any jurisdiction. # -from typing import Iterable, List, Union +from typing import Iterable +from typing import List +from typing import Union import numpy as np @@ -41,9 +43,7 @@ def iter_quantiles( 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'" - ) + raise ValueError(f"Invalid method {method!r}, expected 'sort', 'numpy_bulk', or 'numpy'") if isinstance(which, int): n = which diff --git a/src/earthkit/meteo/thermo/array/thermo.py b/src/earthkit/meteo/thermo/array/thermo.py index f0c2b7a..54e045b 100644 --- a/src/earthkit/meteo/thermo/array/thermo.py +++ b/src/earthkit/meteo/thermo/array/thermo.py @@ -182,9 +182,7 @@ def specific_humidity_from_vapour_pressure(e, p, eps=1e-4): """ if eps <= 0: - raise ValueError( - f"specific_humidity_from_vapour_pressure(): eps={eps} must be > 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 @@ -294,24 +292,14 @@ def _es_mixed(self, t): # 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] - ) + 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) - ) + 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) - ) + 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) @@ -467,9 +455,7 @@ def saturation_vapour_pressure_slope(t, phase="mixed"): return _EsComp().compute_slope(t, phase) -def saturation_mixing_ratio_slope( - t, p, es=None, es_slope=None, phase="mixed", eps=1e-4 -): +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 @@ -519,9 +505,7 @@ def saturation_mixing_ratio_slope( 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 -): +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 @@ -835,9 +819,7 @@ def dewpoint_from_specific_humidity(q, p): used in :func:`saturation_vapour_pressure`. """ - return temperature_from_saturation_vapour_pressure( - vapour_pressure_from_specific_humidity(q, p) - ) + return temperature_from_saturation_vapour_pressure(vapour_pressure_from_specific_humidity(q, p)) def virtual_temperature(t, q): @@ -1057,9 +1039,7 @@ def lcl_temperature(t, td, method="davies"): """ # 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) + 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": @@ -1144,10 +1124,7 @@ def compute_wbpt(self, ept): 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) - ) + 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): @@ -1161,10 +1138,7 @@ def compute_t_on_ma_stipanuk(self, ept, p): 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 - ) + 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 @@ -1211,11 +1185,7 @@ def _D(p): 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[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) @@ -1306,9 +1276,7 @@ def _ept(self, ths): 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) - ) + 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: @@ -1334,9 +1302,7 @@ def _f(self, ths): 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.K3 * np.log(ths.p / constants.p0) * saturation_vapour_pressure_slope(ths.t) + self._d_G_sat(ths) ) @@ -1364,9 +1330,7 @@ def _ept(self, ths): 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 - ) + 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): @@ -1382,36 +1346,24 @@ def _G_sat(self, ths, scale=1.0): # 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) - ) + 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 - ) + 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)) - ) + 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) + + constants.kappa * saturation_vapour_pressure_slope(ths.t) / (ths.p - ths.es) + self._d_G_sat(ths) ) @@ -1606,9 +1558,7 @@ def temperature_on_moist_adiabat(ept, p, ept_method="ifs", t_method="bisect"): 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!" - ) + 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"): @@ -1648,14 +1598,10 @@ def wet_bulb_temperature_from_dewpoint(t, td, p, ept_method="ifs", t_method="bis """ 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 - ) + 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" -): +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 @@ -1693,14 +1639,10 @@ def wet_bulb_temperature_from_specific_humidity( """ 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 - ) + 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" -): +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 @@ -1741,14 +1683,10 @@ def wet_bulb_potential_temperature_from_dewpoint( 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 - ) + 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" -): +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 @@ -1786,6 +1724,4 @@ def wet_bulb_potential_temperature_from_specific_humidity( 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 - ) + return temperature_on_moist_adiabat(ept, constants.p0, ept_method=ept_method, t_method=t_method) diff --git a/src/earthkit/meteo/wind/array/wind.py b/src/earthkit/meteo/wind/array/wind.py index f611ca9..a20a529 100644 --- a/src/earthkit/meteo/wind/array/wind.py +++ b/src/earthkit/meteo/wind/array/wind.py @@ -279,13 +279,9 @@ def windrose(speed, direction, sectors=16, speed_bins=[], percent=True): 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 - ) + 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 - ] + res = np.histogram2d(speed, direction, bins=[speed_bins, dir_bins], density=False)[0] # unify the north bins res[:, 0] = res[:, 0] + res[:, -1] diff --git a/tests/score/test_score.py b/tests/score/test_score.py index aa80c04..ffa8cc1 100644 --- a/tests/score/test_score.py +++ b/tests/score/test_score.py @@ -43,13 +43,9 @@ def crps_quaver2(x, y): lcond = esarr[0, :] > anarr aa[0, lcond] = 1.0 bb[0, :] = np.where(lcond, esarr[0, :] - anarr, 0.0) - aa[1:-1, :] = np.where( - esarr[1:, :] <= anarr, esarr[1:, :] - esarr[:-1, :], anarr - esarr[:-1, :] - ) + aa[1:-1, :] = np.where(esarr[1:, :] <= anarr, esarr[1:, :] - esarr[:-1, :], anarr - esarr[:-1, :]) aa[1:-1, :][esarr[: n_ens - 1, :] > anarr] = 0.0 # this would be hard in xarray - bb[1:-1, :] = np.where( - esarr[:-1, :] > anarr, esarr[1:, :] - esarr[:-1, :], esarr[1:, :] - anarr - ) + bb[1:-1, :] = np.where(esarr[:-1, :] > anarr, esarr[1:, :] - esarr[:-1, :], esarr[1:, :] - anarr) bb[1:-1, :][esarr[1:, :] <= anarr] = 0.0 lcond = anarr > esarr[-1, :] aa[-1, :] = np.where(lcond, anarr - esarr[-1, :], 0.0) diff --git a/tests/solar/test_solar.py b/tests/solar/test_solar.py index bb362b8..64d00f4 100644 --- a/tests/solar/test_solar.py +++ b/tests/solar/test_solar.py @@ -58,9 +58,7 @@ def test_cos_solar_zenith_angle_integrated(): latitudes = np.array([40.0]) longitudes = np.array([18.0]) - v = solar.cos_solar_zenith_angle_integrated( - begin_date, end_date, latitudes, longitudes - ) + v = solar.cos_solar_zenith_angle_integrated(begin_date, end_date, latitudes, longitudes) assert np.isclose(v[0], 0.3110738757) diff --git a/tests/stats/test_stats.py b/tests/stats/test_stats.py index ba7aea4..03aa327 100644 --- a/tests/stats/test_stats.py +++ b/tests/stats/test_stats.py @@ -26,9 +26,7 @@ def test_nanaverage(): # replace nan values target_result[:, 1] = np.nansum(data, axis=-1)[:, 1] print(target_result) - np.testing.assert_equal( - stats.nanaverage(data, axis=-1, weights=weights), target_result - ) + np.testing.assert_equal(stats.nanaverage(data, axis=-1, weights=weights), target_result) @pytest.mark.parametrize("method", ["sort", "numpy_bulk", "numpy"]) @@ -71,10 +69,6 @@ def test_quantiles_nans(): arr = np.random.rand(100, 100, 100) arr.ravel()[np.random.choice(arr.size, 100000, replace=False)] = np.nan qs = [0.0, 0.25, 0.5, 0.75, 1.0] - sort = [ - quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="sort") - ] - numpy = [ - quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="numpy") - ] + sort = [quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="sort")] + numpy = [quantile for quantile in stats.iter_quantiles(arr.copy(), qs, method="numpy")] assert np.all(np.isclose(sort, numpy, equal_nan=True)) diff --git a/tests/thermo/test_thermo_array.py b/tests/thermo/test_thermo_array.py index aa6a792..8a1d75e 100644 --- a/tests/thermo/test_thermo_array.py +++ b/tests/thermo/test_thermo_array.py @@ -302,9 +302,7 @@ def test_saturation_specific_humidity_slope(): ) for phase in phases: - svp = thermo.array.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]) @@ -355,9 +353,7 @@ def test_relative_humidity_from_dewpoint(): def test_relative_humidity_from_specific_humidity(): - t = thermo.array.celsius_to_kelvin( - np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]) - ) + 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 q = np.array( [ @@ -401,9 +397,7 @@ def test_specific_humidity_from_dewpoint(): def test_specific_humidity_from_relative_humidity(): - t = thermo.array.celsius_to_kelvin( - np.array([-29.2884, -14.4118, -5.9235, 9.72339, 18.4514]) - ) + 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 r = np.array( [ @@ -638,9 +632,7 @@ def test_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( - pt, ref[f"{m_ept}_{m_t}"], err_msg=f"method={m_ept}_{m_t}" - ) + np.testing.assert_allclose(pt, ref[f"{m_ept}_{m_t}"], err_msg=f"method={m_ept}_{m_t}") def test_wet_bulb_temperature(): diff --git a/tests/wind/test_wind.py b/tests/wind/test_wind.py index 09f1b1d..3648065 100644 --- a/tests/wind/test_wind.py +++ b/tests/wind/test_wind.py @@ -42,9 +42,7 @@ def test_direction(): # meteo u = np.array([0, 1, 1, 1, 0, -1, -1, -1, 0, np.nan, 1, np.nan]) v = np.array([1, 1, 0, -1, -1, -1, 0, 1, 0, 1, np.nan, np.nan]) - v_ref = np.array( - [180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan] - ) + v_ref = np.array([180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan]) d = wind.direction(u, v, convention="meteo") np.testing.assert_allclose(d, v_ref) @@ -122,9 +120,7 @@ def test_xy_to_polar(): np.nan, ] ) - d_ref = np.array( - [180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan] - ) + d_ref = np.array([180.0, 225, 270, 315, 0, 45, 90, 135, 270, np.nan, np.nan, np.nan]) sp, d = wind.xy_to_polar(u, v) np.testing.assert_allclose(sp, sp_ref) np.testing.assert_allclose(d, d_ref) @@ -172,12 +168,8 @@ def test_coriolis(): def test_windrose(): - sp = np.array( - [3.5, 1, 1.1, 2.1, 0.1, 0.0, 2.4, 1.9, 1.7, 3.9, 3.1, 2.1, np.nan, np.nan] - ) - d = np.array( - [1.0, 29, 31, 93.0, 121, 171, 189, 245, 240.11, 311, 359.1, np.nan, 11, np.nan] - ) + sp = np.array([3.5, 1, 1.1, 2.1, 0.1, 0.0, 2.4, 1.9, 1.7, 3.9, 3.1, 2.1, np.nan, np.nan]) + d = np.array([1.0, 29, 31, 93.0, 121, 171, 189, 245, 240.11, 311, 359.1, np.nan, 11, np.nan]) sp_bins = [0, 1, 2, 3, 4] # count From 3b404c08130ce8389016fea6680cadd5a3c25fd8 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Mon, 4 Nov 2024 11:21:05 +0000 Subject: [PATCH 3/4] Control zero humidity in thermo computations (#23) * Control zero humidity in thermo computations --- src/earthkit/meteo/thermo/array/thermo.py | 93 ++++++++++-- tests/thermo/test_thermo_array.py | 177 +++++++++++++++++----- 2 files changed, 220 insertions(+), 50 deletions(-) diff --git a/src/earthkit/meteo/thermo/array/thermo.py b/src/earthkit/meteo/thermo/array/thermo.py index 54e045b..690f5e7 100644 --- a/src/earthkit/meteo/thermo/array/thermo.py +++ b/src/earthkit/meteo/thermo/array/thermo.py @@ -13,6 +13,10 @@ from earthkit.meteo import constants +def _valid_number(x): + return x is not None and not np.isnan(x) + + def celsius_to_kelvin(t): """Converts temperature values from Celsius to Kelvin. @@ -259,9 +263,21 @@ def compute_slope(self, t, phase): 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 t_from_es(self, es, eps=1e-8, out=None): + def _comp(x): + x = x / self.c1 + x = np.log(x) + return (x * self.c4w - self.c3w * self.t0) / (x - self.c3w) + + if out is not None: + v = np.asarray(es) + z_mask = v < eps + v_mask = ~z_mask + v[v_mask] = _comp(v[v_mask]) + v[z_mask] = out + return v + else: + return _comp(es) def _es_water(self, t): return self.c1 * np.exp(self.c3w * (t - self.t0) / (t - self.c4w)) @@ -556,13 +572,21 @@ def saturation_specific_humidity_slope(t, p, es=None, es_slope=None, phase="mixe return constants.epsilon * es_slope * p / v -def temperature_from_saturation_vapour_pressure(es): - r"""Computes the temperature from saturation vapour pressure. +def temperature_from_saturation_vapour_pressure(es, eps=1e-8, out=None): + r"""Compute the temperature from saturation vapour pressure. Parameters ---------- es: ndarray :func:`saturation_vapour_pressure` (Pa) + eps: number + If ``out`` is not None, return ``out`` when ``es`` < ``eps``. If out + is None, ``eps`` is ignored and return np.nan for ``es`` values very + close to zero. + out: number or None + If not None, return ``out`` when ``es`` < ``eps``. If None, ``eps`` is + ignored and return np.nan for ``es`` values very close to zero. + Returns ------- @@ -575,7 +599,7 @@ def temperature_from_saturation_vapour_pressure(es): phase ``es`` was computed to. """ - return _EsComp().t_from_es(es) + return _EsComp().t_from_es(es, eps=eps, out=out) def relative_humidity_from_dewpoint(t, td): @@ -750,8 +774,8 @@ def specific_humidity_from_relative_humidity(t, r, p): return specific_humidity_from_vapour_pressure(e, p) -def dewpoint_from_relative_humidity(t, r): - r"""Computes the dewpoint temperature from relative humidity. +def dewpoint_from_relative_humidity(t, r, eps=1e-8, out=None): + r"""Compute the dewpoint temperature from relative humidity. Parameters ---------- @@ -759,11 +783,18 @@ def dewpoint_from_relative_humidity(t, r): Temperature (K) r: ndarray Relative humidity (%) + eps: number + If ``out`` is not None, return ``out`` when ``r`` < ``eps``. + If out is None, ``eps`` is ignored and return np.nan for ``r`` + values very close to zero. + out: number or None + If not None, return ``out`` when ``r`` < ``eps``. If None, ``eps`` is + ignored and return np.nan for ``r`` values very close to zero. Returns ------- ndarray - Dewpoint (K) + Dewpoint temperature (K) The computation starts with determining the the saturation vapour pressure over @@ -782,11 +813,24 @@ def dewpoint_from_relative_humidity(t, r): equations used in :func:`saturation_vapour_pressure`. """ + # by masking upfront we avoid RuntimeWarning when calling log() in + # the computation of td when r is very small + if out is not None: + r = np.asarray(r) + mask = r < eps + r[mask] = np.nan + es = saturation_vapour_pressure(t, phase="water") * r / 100.0 - return temperature_from_saturation_vapour_pressure(es) + td = temperature_from_saturation_vapour_pressure(es) + if out is not None and not np.isnan(out): + td = np.asarray(td) + td[mask] = out -def dewpoint_from_specific_humidity(q, p): + return td + + +def dewpoint_from_specific_humidity(q, p, eps=1e-8, out=None): r"""Computes the dewpoint temperature from specific humidity. Parameters @@ -795,11 +839,18 @@ def dewpoint_from_specific_humidity(q, p): Specific humidity (kg/kg) p: ndarray Pressure (Pa) + eps: number + If ``out`` is not None, return ``out`` when ``q`` < ``eps``. + If out is None, ``eps`` is ignored and return np.nan for ``q`` + values very close to zero. + out: number or None + If not None, return ``out`` when ``q`` < ``eps``. If None, ``eps`` is + ignored and return np.nan for ``q`` values very close to zero. Returns ------- ndarray - Dewpoint (K) + Dewpoint temperature (K) The computation starts with determining the the saturation vapour pressure over @@ -819,7 +870,21 @@ def dewpoint_from_specific_humidity(q, p): used in :func:`saturation_vapour_pressure`. """ - return temperature_from_saturation_vapour_pressure(vapour_pressure_from_specific_humidity(q, p)) + # by masking upfront we avoid RuntimeWarning when calling log() in + # the computation of td when q is very small + if out is not None: + q = np.asarray(q) + mask = q < eps + q[mask] = np.nan + + es = vapour_pressure_from_specific_humidity(q, p) + td = temperature_from_saturation_vapour_pressure(es) + + if out is not None and not np.isnan(out): + td = np.asarray(td) + td[mask] = out + + return td def virtual_temperature(t, q): @@ -828,7 +893,7 @@ def virtual_temperature(t, q): Parameters ---------- t: number or ndarray - Temperature (K) + Temperature (K)s q: number or ndarray Specific humidity (kg/kg) diff --git a/tests/thermo/test_thermo_array.py b/tests/thermo/test_thermo_array.py index 8a1d75e..a87fcf8 100644 --- a/tests/thermo/test_thermo_array.py +++ b/tests/thermo/test_thermo_array.py @@ -10,6 +10,7 @@ import os import numpy as np +import pytest from earthkit.meteo import thermo @@ -320,7 +321,7 @@ def test_saturation_specific_humidity_slope(): np.testing.assert_allclose(svp, v_ref[i]) -def test_temperature_from_saturation_vapour_pressure(): +def test_temperature_from_saturation_vapour_pressure_1(): ref_file = "sat_vp.csv" d = np.genfromtxt( data_file(ref_file), @@ -329,7 +330,32 @@ def test_temperature_from_saturation_vapour_pressure(): ) t = thermo.array.temperature_from_saturation_vapour_pressure(d["water"]) - np.testing.assert_allclose(t, d["t"]) + assert np.allclose(t, d["t"], equal_nan=True) + + +@pytest.mark.parametrize( + "es,kwargs,expected_values", + [ + (4.2, {}, 219.7796336743947), + ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": 100}, [219.7796336743947, 100, 100, np.nan]), + ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": np.nan}, [219.7796336743947, np.nan, np.nan, np.nan]), + (0, {}, np.nan), + (0.001, {"eps": 1e-2, "out": 100}, 100.0), + (0.001, {"eps": 1e-2, "out": np.nan}, np.nan), + ], +) +def test_temperature_from_saturation_vapour_pressure_2(es, kwargs, expected_values): + + multi = isinstance(es, list) + if multi: + es = np.array(es) + expected_values = np.array(expected_values) + + t = thermo.array.temperature_from_saturation_vapour_pressure(es, **kwargs) + if multi: + np.testing.assert_allclose(t, expected_values, equal_nan=True) + else: + assert np.isclose(t, expected_values, equal_nan=True) def test_relative_humidity_from_dewpoint(): @@ -421,43 +447,122 @@ def test_specific_humidity_from_relative_humidity(): np.testing.assert_allclose(q, v_ref) -def test_dewpoint_from_relative_humidity(): - t = thermo.array.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25])) - r = np.array( - [ - 100.0000000000, - 52.5224541378, - 46.8714823296, - 84.5391163313, - 21.9244774232, - 46.1081101229, - 15.4779832381, - ] - ) - v_ref = thermo.array.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3])) +@pytest.mark.parametrize( + "t,r,kwargs,expected_values", + [ + ( + [20.0, 20, 0, 35, 5, -15, 25, 25], + [ + 100.0000000000, + 52.5224541378, + 46.8714823296, + 84.5391163313, + 21.9244774232, + 46.1081101229, + 15.4779832381, + 0, + ], + {}, + [20.0, 10, -10, 32, -15, -24, -3, np.nan], + ), + ( + [20.0, 20.0, 20.0], + [ + 52.5224541378, + 0.0, + 0.000001, + ], + {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, + [10, 100, 100], + ), + ( + [20.0, 20.0, 20.0], + [ + 52.5224541378, + 0.0, + 0.000001, + ], + {"eps": 1e-3, "out": np.nan}, + [10, np.nan, np.nan], + ), + (20.0, 52.5224541378, {}, 10.0), + (20.0, 0.0, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100), + (20.0, 0.000001, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100), + (20.0, 0.0, {"eps": 1e-3, "out": np.nan}, np.nan), + (20, 0.000001, {"eps": 1e-3, "out": np.nan}, np.nan), + ], +) +def test_dewpoint_from_relative_humidity(t, r, kwargs, expected_values): # reference was tested with an online relhum calculator at: # https://bmcnoldy.rsmas.miami.edu/Humidity.html - td = thermo.array.dewpoint_from_relative_humidity(t, r) - np.testing.assert_allclose(td, v_ref) - -def test_dewpoint_from_specific_humidity(): - p = np.array([967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207]) * 100 - q = np.array( - [ - 0.0169461501, - 0.0155840075, - 0.0134912382, - 0.0083409720, - 0.0057268584, - 0.0025150791, - ] - ) - v_ref = thermo.array.celsius_to_kelvin( - np.array([21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916]) - ) - td = thermo.array.dewpoint_from_specific_humidity(q, p) - np.testing.assert_allclose(td, v_ref) + multi = isinstance(t, list) + if multi: + t = np.array(t) + r = np.array(r) + expected_values = np.array(expected_values) + + t = thermo.array.celsius_to_kelvin(t) + v_ref = thermo.array.celsius_to_kelvin(expected_values) + + td = thermo.array.dewpoint_from_relative_humidity(t, r, **kwargs) + if multi: + assert np.allclose(td, v_ref, equal_nan=True) + else: + assert np.isclose(td, v_ref, equal_nan=True) + + +@pytest.mark.parametrize( + "q,p,kwargs,expected_values", + [ + ( + [0.0169461501, 0.0155840075, 0.0134912382, 0.0083409720, 0.0057268584, 0.0025150791, 0], + [967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207, 422.4207], + {}, + [21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916, np.nan], + ), + ( + [ + 0.0169461501, + 0.0, + 0.000001, + ], + [967.5085, 967.5085, 967.5085], + {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, + [21.78907, 100, 100], + ), + ( + [ + 0.0169461501, + 0.0, + 0.000001, + ], + [967.5085, 967.5085, 967.5085], + {"eps": 1e-3, "out": np.nan}, + [21.78907, np.nan, np.nan], + ), + (0.0169461501, 967.508, {}, 21.78907), + (0.0, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100), + (0.000001, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100), + (0.0, 967.5085, {"eps": 1e-3, "out": np.nan}, np.nan), + (0.000001, 967.5085, {"eps": 1e-3, "out": np.nan}, np.nan), + ], +) +def test_dewpoint_from_specific_humidity(q, p, kwargs, expected_values): + multi = isinstance(q, list) + if multi: + q = np.array(q) + p = np.array(p) + expected_values = np.array(expected_values) + + p = p * 100.0 + v_ref = thermo.array.celsius_to_kelvin(expected_values) + + td = thermo.array.dewpoint_from_specific_humidity(q, p, **kwargs) + if multi: + assert np.allclose(td, v_ref, equal_nan=True) + else: + assert np.isclose(td, v_ref, equal_nan=True) def test_virtual_temperature(): From ef58809ff822ae2419b9afdbb14cf86a563314e3 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Mon, 4 Nov 2024 11:59:57 +0000 Subject: [PATCH 4/4] Add release notes for version 0.2.0 (#24) --- docs/conf.py | 4 ++++ docs/release_notes/index.rst | 1 + docs/release_notes/version_0.2_updates.rst | 15 +++++++++++++++ docs/requirements.txt | 1 + environment.yml | 1 + src/earthkit/meteo/thermo/array/thermo.py | 11 +++++++---- tests/environment-unit-tests.yml | 1 + 7 files changed, 30 insertions(+), 4 deletions(-) create mode 100644 docs/release_notes/version_0.2_updates.rst diff --git a/docs/conf.py b/docs/conf.py index 4854c59..fa0b391 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -36,6 +36,7 @@ "sphinx.ext.autodoc", "sphinx.ext.napoleon", "autoapi.extension", + "sphinx_issues", ] # autodoc configuration @@ -62,6 +63,9 @@ napoleon_numpy_docstring = True napoleon_preprocess_types = True +# Path to GitHub repo {group}/{project} (note that `group` is the GitHub user or organization) +issues_github_path = "ecmwf/earthkit-meteo" + # Add any paths that contain templates here, relative to this directory. templates_path = ["_templates"] diff --git a/docs/release_notes/index.rst b/docs/release_notes/index.rst index 37ca6f6..3e54796 100644 --- a/docs/release_notes/index.rst +++ b/docs/release_notes/index.rst @@ -4,4 +4,5 @@ Release notes .. toctree:: :maxdepth: 1 + version_0.2_updates version_0.1_updates diff --git a/docs/release_notes/version_0.2_updates.rst b/docs/release_notes/version_0.2_updates.rst new file mode 100644 index 0000000..70dc830 --- /dev/null +++ b/docs/release_notes/version_0.2_updates.rst @@ -0,0 +1,15 @@ +Version 0.2 Updates +///////////////////////// + + +Version 0.2.0 +=============== + +New features ++++++++++++++++ + +- Added the ``eps`` and ``out`` keywords arguments to control zero humidity in the following methods (:pr:`23`): + + - :py:meth:`dewpoint_from_specific_humidity ` + - :py:meth:`dewpoint_from_relative_humidity ` + - :py:meth:`temperature_from_saturation_vapour_pressure ` diff --git a/docs/requirements.txt b/docs/requirements.txt index bcfdc7d..ae01e54 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -7,3 +7,4 @@ sphinx-rtd-theme nbsphinx setuptools sphinx-autoapi +sphinx-issues diff --git a/environment.yml b/environment.yml index ec9248d..6511388 100644 --- a/environment.yml +++ b/environment.yml @@ -17,6 +17,7 @@ dependencies: - sphinx-autoapi - sphinx_rtd_theme - sphinxcontrib-apidoc +- sphinx-issues - nbformat - nbconvert - nbsphinx diff --git a/src/earthkit/meteo/thermo/array/thermo.py b/src/earthkit/meteo/thermo/array/thermo.py index 690f5e7..4155d1a 100644 --- a/src/earthkit/meteo/thermo/array/thermo.py +++ b/src/earthkit/meteo/thermo/array/thermo.py @@ -582,10 +582,11 @@ def temperature_from_saturation_vapour_pressure(es, eps=1e-8, out=None): eps: number If ``out`` is not None, return ``out`` when ``es`` < ``eps``. If out is None, ``eps`` is ignored and return np.nan for ``es`` values very - close to zero. + close to zero. *New in version 0.2.0* out: number or None If not None, return ``out`` when ``es`` < ``eps``. If None, ``eps`` is ignored and return np.nan for ``es`` values very close to zero. + *New in version 0.2.0* Returns @@ -786,10 +787,11 @@ def dewpoint_from_relative_humidity(t, r, eps=1e-8, out=None): eps: number If ``out`` is not None, return ``out`` when ``r`` < ``eps``. If out is None, ``eps`` is ignored and return np.nan for ``r`` - values very close to zero. + values very close to zero. *New in version 0.2.0* out: number or None If not None, return ``out`` when ``r`` < ``eps``. If None, ``eps`` is ignored and return np.nan for ``r`` values very close to zero. + *New in version 0.2.0* Returns ------- @@ -831,7 +833,7 @@ def dewpoint_from_relative_humidity(t, r, eps=1e-8, out=None): def dewpoint_from_specific_humidity(q, p, eps=1e-8, out=None): - r"""Computes the dewpoint temperature from specific humidity. + r"""Compute the dewpoint temperature from specific humidity. Parameters ---------- @@ -842,10 +844,11 @@ def dewpoint_from_specific_humidity(q, p, eps=1e-8, out=None): eps: number If ``out`` is not None, return ``out`` when ``q`` < ``eps``. If out is None, ``eps`` is ignored and return np.nan for ``q`` - values very close to zero. + values very close to zero. *New in version 0.2.0* out: number or None If not None, return ``out`` when ``q`` < ``eps``. If None, ``eps`` is ignored and return np.nan for ``q`` values very close to zero. + *New in version 0.2.0* Returns ------- diff --git a/tests/environment-unit-tests.yml b/tests/environment-unit-tests.yml index 17eefa1..78e3148 100644 --- a/tests/environment-unit-tests.yml +++ b/tests/environment-unit-tests.yml @@ -16,6 +16,7 @@ dependencies: - sphinx-autoapi - sphinx_rtd_theme - sphinxcontrib-apidoc +- sphinx-issues - nbformat - nbconvert - nbsphinx