diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1ac1cf07..53fa47b3 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -8,41 +8,21 @@ repos: - id: check-yaml # Check YAML files for syntax errors only args: [--unsafe, --allow-multiple-documents] - id: check-toml - # - id: check-added-large-files + # - 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.13.2 - hooks: - - id: isort - args: - - -l 110 - - --force-single-line-imports - - --profile black -- repo: https://github.com/psf/black - 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] - exclude: xr_engine_profile_rst\.py - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.6.9 + rev: v0.15.4 hooks: - - id: ruff + - id: ruff-check exclude: '(dev/.*|.*_)\.py$' args: - - --line-length=110 - --fix - --exit-non-zero-on-fix - - --preview + - id: ruff-format - repo: https://github.com/executablebooks/mdformat rev: 0.7.14 hooks: @@ -53,11 +33,9 @@ repos: hooks: - id: pretty-format-yaml args: [--autofix, --preserve-quotes] + - id: pretty-format-toml + args: [--autofix] - 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: "v2.5.0" - hooks: - - id: pyproject-fmt diff --git a/docs/examples/hybrid_levels.ipynb b/docs/examples/hybrid_levels.ipynb index f48d229c..3414010d 100644 --- a/docs/examples/hybrid_levels.ipynb +++ b/docs/examples/hybrid_levels.ipynb @@ -131,10 +131,10 @@ "from earthkit.meteo.utils.sample import get_sample\n", "\n", "DATA = get_sample(\"vertical_hybrid_data\")\n", - "sp = DATA.p_surf # surface pressure [Pa]\n", - "zs = DATA.z_surf # surface geopotential [m2/s2]\n", - "t = DATA.t # temperature [K]\n", - "q = DATA.q # specific humidity [kg/kg]\n", + "sp = DATA.p_surf # surface pressure [Pa]\n", + "zs = DATA.z_surf # surface geopotential [m2/s2]\n", + "t = DATA.t # temperature [K]\n", + "q = DATA.q # specific humidity [kg/kg]\n", "\n", "sp.shape, zs.shape, t.shape, q.shape" ] @@ -329,9 +329,7 @@ } ], "source": [ - "p_full, p_half, alpha,delta = vertical.pressure_on_hybrid_levels(\n", - " A, B, sp, \n", - " output=(\"full\", \"half\",\"alpha\", \"delta\"))\n", + "p_full, p_half, alpha, delta = vertical.pressure_on_hybrid_levels(A, B, sp, output=(\"full\", \"half\", \"alpha\", \"delta\"))\n", "\n", "p_full.shape, p_half.shape, alpha.shape, delta.shape" ] @@ -376,9 +374,7 @@ } ], "source": [ - "p_full = vertical.pressure_on_hybrid_levels(\n", - " A, B, sp, levels=[135, 136, 137],\n", - " output=(\"full\"))\n", + "p_full = vertical.pressure_on_hybrid_levels(A, B, sp, levels=[135, 136, 137], output=(\"full\"))\n", "p_full" ] }, @@ -422,9 +418,7 @@ } ], "source": [ - "p_full = vertical.pressure_on_hybrid_levels(\n", - " A, B, sp, levels=[137, 136, 135],\n", - " output=(\"full\"))\n", + "p_full = vertical.pressure_on_hybrid_levels(A, B, sp, levels=[137, 136, 135], output=(\"full\"))\n", "p_full" ] }, @@ -470,8 +464,7 @@ }, "outputs": [], "source": [ - "z_thickness= vertical.relative_geopotential_thickness_on_hybrid_levels(t, q, \n", - " A, B, sp)" + "z_thickness = vertical.relative_geopotential_thickness_on_hybrid_levels(t, q, A, B, sp)" ] }, { @@ -514,7 +507,7 @@ } ], "source": [ - "z_thickness= vertical.relative_geopotential_thickness_on_hybrid_levels(t[-3:], q[-3:], A, B, sp)\n", + "z_thickness = vertical.relative_geopotential_thickness_on_hybrid_levels(t[-3:], q[-3:], A, B, sp)\n", "z_thickness" ] }, diff --git a/docs/examples/interpolate_hybrid_to_hl.ipynb b/docs/examples/interpolate_hybrid_to_hl.ipynb index d745f78f..b6066120 100644 --- a/docs/examples/interpolate_hybrid_to_hl.ipynb +++ b/docs/examples/interpolate_hybrid_to_hl.ipynb @@ -1,1047 +1,1074 @@ { - "cells": [ - { - "cell_type": "markdown", - "id": "150548d1-f025-47bf-9604-31a339c6bc91", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "## Interpolating from hybrid to height levels" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "653672d0-0d06-40bd-8945-7ec046bfc93a", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [], - "source": [ - "import earthkit.meteo.vertical.array as vertical" - ] - }, - { - "cell_type": "markdown", - "id": "10e033ec-653b-41f5-aac0-02200a27b8f6", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "#### Hybrid levels" - ] - }, - { - "cell_type": "raw", - "id": "f4d807e1-f3f1-4f0f-bcda-1f7718735898", - "metadata": { - "editable": true, - "raw_mimetype": "text/restructuredtext", - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "Hybrid levels are terrain following at the bottom and isobaric at the top of the atmosphere with a transition in between. This is the vertical coordinate system in the IFS model. \n", - "\n", - "For details about the hybrid levels see `IFS Documentation CY47R3 - Part IV Physical processes, Chapter 2, Section 2.2.1. `_. Also see the :ref:`/examples/hybrid_levels.ipynb` notebook for the related computations in earthkit-meteo." - ] - }, - { - "cell_type": "markdown", - "id": "cf162a36-b766-4467-8796-452aa7f70927", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "#### Getting the data" - ] - }, - { - "cell_type": "markdown", - "id": "5244b1be-9462-4cee-8fba-bcc105cc3fbb", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "First, we get some sample data containing 137 IFS model levels (hybrid full-levels) for two points. This data is in ascending order with regards to the model level number. So the first level is model level 1 (i.e. the top), while the last one is model level 137 (i.e. the bottom)." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "966eb14e-fde6-423a-b13d-a9a5274b832c", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "data": { - "text/plain": [ - "((2,), (2,), (137, 2), (137, 2))" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "from earthkit.meteo.utils.sample import get_sample\n", - "\n", - "DATA = get_sample(\"vertical_hybrid_data\")\n", - "sp = DATA.p_surf # surface pressure [Pa]\n", - "zs = DATA.z_surf # surface geopotential [m2/s2]\n", - "t = DATA.t # temperature [K]\n", - "q = DATA.q # specific humidity [kg/kg]\n", - "\n", - "sp.shape, zs.shape, t.shape, q.shape" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "72e1953c-a1da-438e-858d-d14b6b6c09c1", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[197.06175232, 198.54808044],\n", - " [230.29292297, 219.00386047],\n", - " [234.56352234, 229.85160828],\n", - " [264.58778381, 292.04481506]])" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "t[::40]" - ] - }, - { - "cell_type": "raw", - "id": "15bf8a34-a571-421d-82ea-f3a46fc4d753", - "metadata": { - "editable": true, - "raw_mimetype": "text/restructuredtext", - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "Next, get the hybrid level definitions using :py:meth:`~earthkit.meteo.vertical.array.hybrid_level_parameters`." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "02c083b4-d38f-4e6b-99d3-39fc6bdca336", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "data": { - "text/plain": [ - "((138,), (138,))" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "A, B = vertical.hybrid_level_parameters(137, model=\"ifs\")\n", - "A.shape, B.shape" - ] - }, - { - "cell_type": "markdown", - "id": "0420e042-eb20-4181-8750-b8c68db4c922", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "#### Using interpolate_hybrid_to_height_levels" - ] - }, - { - "cell_type": "raw", - "id": "3efd7623-b84a-48ff-a53c-642a796bc3ae", - "metadata": { - "editable": true, - "raw_mimetype": "text/restructuredtext", - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "We can use :py:meth:`~earthkit.meteo.vertical.array.interpolate_hybrid_to_height_levels` for hybrid to height level interpolation. The example below interpolates temperature to geometric heights above the ground." - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "a728b121-4642-46d6-a360-3727ede0b6aa", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [262.36938948 291.44520344]\n", - "5000.0 [236.77461002 265.49528592]\n" - ] - } - ], - "source": [ - "target_h = [1000., 5000.] # # geometric height, above the ground\n", - "\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\",\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "b025c0ce-12a8-40d1-8311-9fff11873b33", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Heights above sea level" - ] - }, - { - "cell_type": "markdown", - "id": "e78b8f4e-8873-4659-b154-be28e23d84a3", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "The reference height can also be the sea level." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "a5bf08f5-342b-4955-9dc8-9f7d4f7cc766", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [265.83447529 291.04194846]\n", - "5000.0 [239.80992741 264.96290891]\n" - ] - } - ], - "source": [ - "target_h = [1000., 5000.] # # geometric height, above the sea\n", - "\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"sea\",\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "f999af82-ca46-49bf-9ef3-b9fde15628b9", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Geopotential heights" - ] - }, - { - "cell_type": "markdown", - "id": "3594a0ed-361a-4fc7-93bd-57eeea43a653", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "The height type can be geopotential height." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "46e94b8d-e002-44d4-a5c3-aa544cc4a867", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [262.36579436 291.44471712]\n", - "5000.0 [236.7517288 265.47139844]\n" - ] - } - ], - "source": [ - "target_h = [1000., 5000.] # # geopotential height, above the ground\n", - "\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " t, q, zs, A, B, sp, \n", - " h_type=\"geopotential\", \n", - " h_reference=\"ground\",\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "4e89d7b9-823a-406d-a739-f8835757c9ad", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Using a subset of levels" - ] - }, - { - "cell_type": "markdown", - "id": "e651ed62-ea64-4797-b0f9-94de2a13794e", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "It is possible to use only a **subset of the model levels** in the input data. This can significantly speed up the computations and reduce memory usage. \n", - "\n", - "In this case the model level range in the input must be contiguous and include the bottom-most level. The following example only uses the lowest 50 model levels above the surface. " - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "8e215841-9355-43fd-8108-9fa0462579c4", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [262.36938948 291.44520344]\n", - "5000.0 [236.77461002 265.49528592]\n" - ] - } - ], - "source": [ - "# using the 50 lowest model levels as input\n", - "# Please note we still need to use the full A and B arrays\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t[-50:], # data to interpolate\n", - " target_h, \n", - " t[-50:], q[-50:], zs, A, B, sp,\n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\",\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "dae03c01-1a3a-4d31-b2d9-ed820f06946a", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Using aux levels" - ] - }, - { - "cell_type": "markdown", - "id": "a831c226-fbd2-471d-8e0b-c5edff42bf02", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "Since, the lowest (full) model level is always above the surface the interpolation fails if the target height is between the surface and the lowest model level. We can compute the height of the lowest model level in our data:" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "91cc251d-88ee-472f-9dd5-2315c1e23204", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[ 9.34500108, 10.21640974]])" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "vertical.height_on_hybrid_levels(t[-1:], q[-1:], zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\")" - ] - }, - { - "cell_type": "markdown", - "id": "756e06c4-294c-40e4-ab67-8e8bbaa347f7", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "As a consequence, interpolation to e.g. 5 m above the ground does not work:" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "id": "69b6fb60-f40e-4471-8ca7-edc680b47194", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "5.0 [nan nan]\n" - ] - } - ], - "source": [ - "target_h = [5.] # # geometric height, above the ground\n", - "\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\",\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "6330b804-b753-47c6-b4c0-9fe061e9da15", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "To overcome this problem we can define \"aux\" levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the ``aux_bottom_*`` kwargs." - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "ccf15534-d40b-42f8-b2b9-84071c134a67", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "5.0 [274.04815857 296.50007348]\n" - ] - } - ], - "source": [ - "target_h = [5.] # # geometric height, above the ground\n", - "\n", - "t_res = vertical.interpolate_hybrid_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\",\n", - " aux_bottom_h=0.,\n", - " aux_bottom_data=[280., 300.],\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "926c8528-6429-4302-98f3-e643012b8938", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "#### Using interpolate_monotonic" - ] - }, - { - "cell_type": "raw", - "id": "88a37616-4861-498a-b94e-746b428be2fa", - "metadata": { - "editable": true, - "raw_mimetype": "text/restructuredtext", - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "We can also use :py:meth:`~earthkit.meteo.vertical.array.interpolate_monotonic` to carry out the interpolation. This is a generic method and we need an extra step to compute the height on (full) model levels using :py:meth:`~meteo.vertical.array.height_on_hybrid_levels`. This height data then can be used to interpolate multiple input data arrays." - ] - }, - { - "cell_type": "markdown", - "id": "d407f188-c3dd-4f0f-8965-0f7f08ae29cc", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "First, compute the geometric height above the ground on all the model levels in the input data. " - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "cdd71bf5-4755-437a-b0d4-333e84afb337", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "data": { - "text/plain": [ - "((137, 2),\n", - " array([[81782.73207519, 80656.03851153],\n", - " [23758.58686088, 24154.77147989],\n", - " [ 9064.8524558 , 10238.94032239],\n", - " [ 654.70913487, 736.54960026]]))" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "h = vertical.height_on_hybrid_levels(t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\")\n", - "h.shape, h[::40]" - ] - }, - { - "cell_type": "markdown", - "id": "e4d2d855-b0d7-463f-9437-0ca0882f24d4", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "Next, interpolate the temperature to the target height levels." - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "id": "7d3107c7-36b5-4684-b081-812e90f26213", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [262.36938948 291.44520344]\n", - "5000.0 [236.77461002 265.49528592]\n" - ] - } - ], - "source": [ - "target_h=[1000., 5000.] # geometric height (m), above the ground\n", - "\n", - "t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "472a9462-c037-4a3f-8fda-f8997ec43cbb", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "Repeat the same computation for specific humidity." - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "id": "da1dfe47-31e2-4ff8-97db-b381a3cc86e9", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [0.00112569 0.00748045]\n", - "5000.0 [7.55903873e-05 1.37677882e-03]\n" - ] - } - ], - "source": [ - "q_res = vertical.interpolate_monotonic(q, h, target_h, interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, q_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "27ec0863-136a-44f8-8df0-bac928786117", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Heights above sea level" - ] - }, - { - "cell_type": "raw", - "id": "fe2b3bac-730c-49d1-b8a9-ce059f2a1668", - "metadata": { - "editable": true, - "raw_mimetype": "text/restructuredtext", - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "The reference height can also be the sea level. This time we need to pass the proper surface geopotential to :py:meth:`~meteo.vertical.array.height_on_hybrid_levels`." - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "id": "3ec0ce27-1b4d-4ae5-98e8-3d0dde568769", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [265.83447529 291.04194846]\n", - "5000.0 [239.80992741 264.96290891]\n" - ] - } - ], - "source": [ - "h_sea = vertical.height_on_hybrid_levels(t, q, zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"sea\")\n", - "\n", - "target_h_sea=[1000., 5000.] # geometric height, above the sea\n", - "\n", - "t_res = vertical.interpolate_monotonic(t, h_sea, target_h_sea, interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h_sea, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "b991af5a-251b-4cb1-83d1-15621e315c6b", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Using a subset of levels" - ] - }, - { - "cell_type": "markdown", - "id": "b06f9596-af2d-4419-8e74-66295c1d7563", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "We can also use a **subset of levels**. This can significantly speed up the computations and reduce memory usage.\n", - "\n", - "In this case the model level range in the input data must be contiguous and include the bottom-most level. The example below only uses the lowest 50 model levels above the surface." - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "id": "a874cfee-fd79-4842-8893-0dc038a3a4db", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1000.0 [262.36938948 291.44520344]\n", - "5000.0 [236.77461002 265.49528592]\n" - ] - } - ], - "source": [ - "# compute height on the 50 lowest model levels\n", - "# Please note we still need to use the full A and B arrays.\n", - "h_sub = vertical.height_on_hybrid_levels(t[-50:], q[-50:], zs, A, B, sp, \n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\")\n", - "\n", - "# interpolate to height levels\n", - "t_res = vertical.interpolate_monotonic(t[-50:], h_sub, target_h, interpolation=\"linear\")\n", - "\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "b97579dd-35e8-4d9a-a4d9-e1649ae648da", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "##### Using aux levels" - ] - }, - { - "cell_type": "markdown", - "id": "806558b7-a7a4-4db0-a71d-1f0176379284", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "As described above, the lowest (full) model level is always above the surface so the interpolation fails if the target height is between the surface and the lowest model level." - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "id": "6b536af5-56e2-41c5-858f-28db8205b775", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "5.0 [nan nan]\n" - ] - } - ], - "source": [ - "target_h=[5.] # geometric height, above the ground\n", - "\n", - "t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "markdown", - "id": "b154ee25-5802-4215-95cc-6e5376373f69", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "source": [ - "To overcome this problem we can define \"aux\" levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the ``aux_min_level_*`` kwargs." - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "id": "1f5e41a6-9835-4afd-b971-95e37af67e1a", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "5.0 [274.04815857 296.50007348]\n" - ] - } - ], - "source": [ - "target_h=[5.] # geometric height, above the ground\n", - "\n", - "t_res = vertical.interpolate_monotonic(t, h, target_h, \n", - " aux_min_level_coord=0,\n", - " aux_min_level_data=[280., 300.],\n", - " interpolation=\"linear\")\n", - "\n", - "for hv, rv in zip(target_h, t_res):\n", - " print(hv, rv)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bd2fc482-88a5-4edb-8b9c-77e8496202d8", - "metadata": { - "editable": true, - "slideshow": { - "slide_type": "" - }, - "tags": [] - }, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "dev", - "language": "python", - "name": "dev" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.13.1" - } + "cells": [ + { + "cell_type": "markdown", + "id": "150548d1-f025-47bf-9604-31a339c6bc91", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" }, - "nbformat": 4, - "nbformat_minor": 5 + "tags": [] + }, + "source": [ + "## Interpolating from hybrid to height levels" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "653672d0-0d06-40bd-8945-7ec046bfc93a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "import earthkit.meteo.vertical.array as vertical" + ] + }, + { + "cell_type": "markdown", + "id": "10e033ec-653b-41f5-aac0-02200a27b8f6", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "#### Hybrid levels" + ] + }, + { + "cell_type": "raw", + "id": "f4d807e1-f3f1-4f0f-bcda-1f7718735898", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "Hybrid levels are terrain following at the bottom and isobaric at the top of the atmosphere with a transition in between. This is the vertical coordinate system in the IFS model. \n", + "\n", + "For details about the hybrid levels see `IFS Documentation CY47R3 - Part IV Physical processes, Chapter 2, Section 2.2.1. `_. Also see the :ref:`/examples/hybrid_levels.ipynb` notebook for the related computations in earthkit-meteo." + ] + }, + { + "cell_type": "markdown", + "id": "cf162a36-b766-4467-8796-452aa7f70927", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "#### Getting the data" + ] + }, + { + "cell_type": "markdown", + "id": "5244b1be-9462-4cee-8fba-bcc105cc3fbb", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "First, we get some sample data containing 137 IFS model levels (hybrid full-levels) for two points. This data is in ascending order with regards to the model level number. So the first level is model level 1 (i.e. the top), while the last one is model level 137 (i.e. the bottom)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "966eb14e-fde6-423a-b13d-a9a5274b832c", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((2,), (2,), (137, 2), (137, 2))" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from earthkit.meteo.utils.sample import get_sample\n", + "\n", + "DATA = get_sample(\"vertical_hybrid_data\")\n", + "sp = DATA.p_surf # surface pressure [Pa]\n", + "zs = DATA.z_surf # surface geopotential [m2/s2]\n", + "t = DATA.t # temperature [K]\n", + "q = DATA.q # specific humidity [kg/kg]\n", + "\n", + "sp.shape, zs.shape, t.shape, q.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "72e1953c-a1da-438e-858d-d14b6b6c09c1", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[197.06175232, 198.54808044],\n", + " [230.29292297, 219.00386047],\n", + " [234.56352234, 229.85160828],\n", + " [264.58778381, 292.04481506]])" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "t[::40]" + ] + }, + { + "cell_type": "raw", + "id": "15bf8a34-a571-421d-82ea-f3a46fc4d753", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "Next, get the hybrid level definitions using :py:meth:`~earthkit.meteo.vertical.array.hybrid_level_parameters`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "02c083b4-d38f-4e6b-99d3-39fc6bdca336", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((138,), (138,))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A, B = vertical.hybrid_level_parameters(137, model=\"ifs\")\n", + "A.shape, B.shape" + ] + }, + { + "cell_type": "markdown", + "id": "0420e042-eb20-4181-8750-b8c68db4c922", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "#### Using interpolate_hybrid_to_height_levels" + ] + }, + { + "cell_type": "raw", + "id": "3efd7623-b84a-48ff-a53c-642a796bc3ae", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can use :py:meth:`~earthkit.meteo.vertical.array.interpolate_hybrid_to_height_levels` for hybrid to height level interpolation. The example below interpolates temperature to geometric heights above the ground." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a728b121-4642-46d6-a360-3727ede0b6aa", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [262.36938948 291.44520344]\n", + "5000.0 [236.77461002 265.49528592]\n" + ] + } + ], + "source": [ + "target_h = [1000.0, 5000.0] # # geometric height, above the ground\n", + "\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t, # data to interpolate\n", + " target_h,\n", + " t,\n", + " q,\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "b025c0ce-12a8-40d1-8311-9fff11873b33", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Heights above sea level" + ] + }, + { + "cell_type": "markdown", + "id": "e78b8f4e-8873-4659-b154-be28e23d84a3", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "The reference height can also be the sea level." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "a5bf08f5-342b-4955-9dc8-9f7d4f7cc766", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [265.83447529 291.04194846]\n", + "5000.0 [239.80992741 264.96290891]\n" + ] + } + ], + "source": [ + "target_h = [1000.0, 5000.0] # # geometric height, above the sea\n", + "\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t, # data to interpolate\n", + " target_h,\n", + " t,\n", + " q,\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geometric\",\n", + " h_reference=\"sea\",\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "f999af82-ca46-49bf-9ef3-b9fde15628b9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Geopotential heights" + ] + }, + { + "cell_type": "markdown", + "id": "3594a0ed-361a-4fc7-93bd-57eeea43a653", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "The height type can be geopotential height." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "46e94b8d-e002-44d4-a5c3-aa544cc4a867", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [262.36579436 291.44471712]\n", + "5000.0 [236.7517288 265.47139844]\n" + ] + } + ], + "source": [ + "target_h = [1000.0, 5000.0] # # geopotential height, above the ground\n", + "\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t, # data to interpolate\n", + " target_h,\n", + " t,\n", + " q,\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geopotential\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "4e89d7b9-823a-406d-a739-f8835757c9ad", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Using a subset of levels" + ] + }, + { + "cell_type": "markdown", + "id": "e651ed62-ea64-4797-b0f9-94de2a13794e", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "It is possible to use only a **subset of the model levels** in the input data. This can significantly speed up the computations and reduce memory usage. \n", + "\n", + "In this case the model level range in the input must be contiguous and include the bottom-most level. The following example only uses the lowest 50 model levels above the surface. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8e215841-9355-43fd-8108-9fa0462579c4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [262.36938948 291.44520344]\n", + "5000.0 [236.77461002 265.49528592]\n" + ] + } + ], + "source": [ + "# using the 50 lowest model levels as input\n", + "# Please note we still need to use the full A and B arrays\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t[-50:], # data to interpolate\n", + " target_h,\n", + " t[-50:],\n", + " q[-50:],\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "dae03c01-1a3a-4d31-b2d9-ed820f06946a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Using aux levels" + ] + }, + { + "cell_type": "markdown", + "id": "a831c226-fbd2-471d-8e0b-c5edff42bf02", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "Since, the lowest (full) model level is always above the surface the interpolation fails if the target height is between the surface and the lowest model level. We can compute the height of the lowest model level in our data:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "91cc251d-88ee-472f-9dd5-2315c1e23204", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 9.34500108, 10.21640974]])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vertical.height_on_hybrid_levels(t[-1:], q[-1:], zs, A, B, sp, h_type=\"geometric\", h_reference=\"ground\")" + ] + }, + { + "cell_type": "markdown", + "id": "756e06c4-294c-40e4-ab67-8e8bbaa347f7", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "As a consequence, interpolation to e.g. 5 m above the ground does not work:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "69b6fb60-f40e-4471-8ca7-edc680b47194", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.0 [nan nan]\n" + ] + } + ], + "source": [ + "target_h = [5.0] # # geometric height, above the ground\n", + "\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t, # data to interpolate\n", + " target_h,\n", + " t,\n", + " q,\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "6330b804-b753-47c6-b4c0-9fe061e9da15", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "To overcome this problem we can define \"aux\" levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the ``aux_bottom_*`` kwargs." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "ccf15534-d40b-42f8-b2b9-84071c134a67", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.0 [274.04815857 296.50007348]\n" + ] + } + ], + "source": [ + "target_h = [5.0] # # geometric height, above the ground\n", + "\n", + "t_res = vertical.interpolate_hybrid_to_height_levels(\n", + " t, # data to interpolate\n", + " target_h,\n", + " t,\n", + " q,\n", + " zs,\n", + " A,\n", + " B,\n", + " sp,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " aux_bottom_h=0.0,\n", + " aux_bottom_data=[280.0, 300.0],\n", + " interpolation=\"linear\",\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "926c8528-6429-4302-98f3-e643012b8938", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "#### Using interpolate_monotonic" + ] + }, + { + "cell_type": "raw", + "id": "88a37616-4861-498a-b94e-746b428be2fa", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can also use :py:meth:`~earthkit.meteo.vertical.array.interpolate_monotonic` to carry out the interpolation. This is a generic method and we need an extra step to compute the height on (full) model levels using :py:meth:`~meteo.vertical.array.height_on_hybrid_levels`. This height data then can be used to interpolate multiple input data arrays." + ] + }, + { + "cell_type": "markdown", + "id": "d407f188-c3dd-4f0f-8965-0f7f08ae29cc", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "First, compute the geometric height above the ground on all the model levels in the input data. " + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "cdd71bf5-4755-437a-b0d4-333e84afb337", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((137, 2),\n", + " array([[81782.73207519, 80656.03851153],\n", + " [23758.58686088, 24154.77147989],\n", + " [ 9064.8524558 , 10238.94032239],\n", + " [ 654.70913487, 736.54960026]]))" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "h = vertical.height_on_hybrid_levels(t, q, zs, A, B, sp, h_type=\"geometric\", h_reference=\"ground\")\n", + "h.shape, h[::40]" + ] + }, + { + "cell_type": "markdown", + "id": "e4d2d855-b0d7-463f-9437-0ca0882f24d4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "Next, interpolate the temperature to the target height levels." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "7d3107c7-36b5-4684-b081-812e90f26213", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [262.36938948 291.44520344]\n", + "5000.0 [236.77461002 265.49528592]\n" + ] + } + ], + "source": [ + "target_h = [1000.0, 5000.0] # geometric height (m), above the ground\n", + "\n", + "t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation=\"linear\")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "472a9462-c037-4a3f-8fda-f8997ec43cbb", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "Repeat the same computation for specific humidity." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "da1dfe47-31e2-4ff8-97db-b381a3cc86e9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [0.00112569 0.00748045]\n", + "5000.0 [7.55903873e-05 1.37677882e-03]\n" + ] + } + ], + "source": [ + "q_res = vertical.interpolate_monotonic(q, h, target_h, interpolation=\"linear\")\n", + "\n", + "for hv, rv in zip(target_h, q_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "27ec0863-136a-44f8-8df0-bac928786117", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Heights above sea level" + ] + }, + { + "cell_type": "raw", + "id": "fe2b3bac-730c-49d1-b8a9-ce059f2a1668", + "metadata": { + "editable": true, + "raw_mimetype": "text/restructuredtext", + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "The reference height can also be the sea level. This time we need to pass the proper surface geopotential to :py:meth:`~meteo.vertical.array.height_on_hybrid_levels`." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "3ec0ce27-1b4d-4ae5-98e8-3d0dde568769", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [265.83447529 291.04194846]\n", + "5000.0 [239.80992741 264.96290891]\n" + ] + } + ], + "source": [ + "h_sea = vertical.height_on_hybrid_levels(t, q, zs, A, B, sp, h_type=\"geometric\", h_reference=\"sea\")\n", + "\n", + "target_h_sea = [1000.0, 5000.0] # geometric height, above the sea\n", + "\n", + "t_res = vertical.interpolate_monotonic(t, h_sea, target_h_sea, interpolation=\"linear\")\n", + "\n", + "for hv, rv in zip(target_h_sea, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "b991af5a-251b-4cb1-83d1-15621e315c6b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Using a subset of levels" + ] + }, + { + "cell_type": "markdown", + "id": "b06f9596-af2d-4419-8e74-66295c1d7563", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can also use a **subset of levels**. This can significantly speed up the computations and reduce memory usage.\n", + "\n", + "In this case the model level range in the input data must be contiguous and include the bottom-most level. The example below only uses the lowest 50 model levels above the surface." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "a874cfee-fd79-4842-8893-0dc038a3a4db", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000.0 [262.36938948 291.44520344]\n", + "5000.0 [236.77461002 265.49528592]\n" + ] + } + ], + "source": [ + "# compute height on the 50 lowest model levels\n", + "# Please note we still need to use the full A and B arrays.\n", + "h_sub = vertical.height_on_hybrid_levels(t[-50:], q[-50:], zs, A, B, sp, h_type=\"geometric\", h_reference=\"ground\")\n", + "\n", + "# interpolate to height levels\n", + "t_res = vertical.interpolate_monotonic(t[-50:], h_sub, target_h, interpolation=\"linear\")\n", + "\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "b97579dd-35e8-4d9a-a4d9-e1649ae648da", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "##### Using aux levels" + ] + }, + { + "cell_type": "markdown", + "id": "806558b7-a7a4-4db0-a71d-1f0176379284", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "As described above, the lowest (full) model level is always above the surface so the interpolation fails if the target height is between the surface and the lowest model level." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "6b536af5-56e2-41c5-858f-28db8205b775", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.0 [nan nan]\n" + ] + } + ], + "source": [ + "target_h = [5.0] # geometric height, above the ground\n", + "\n", + "t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation=\"linear\")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "markdown", + "id": "b154ee25-5802-4215-95cc-6e5376373f69", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "To overcome this problem we can define \"aux\" levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the ``aux_min_level_*`` kwargs." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "1f5e41a6-9835-4afd-b971-95e37af67e1a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.0 [274.04815857 296.50007348]\n" + ] + } + ], + "source": [ + "target_h = [5.0] # geometric height, above the ground\n", + "\n", + "t_res = vertical.interpolate_monotonic(\n", + " t, h, target_h, aux_min_level_coord=0, aux_min_level_data=[280.0, 300.0], interpolation=\"linear\"\n", + ")\n", + "\n", + "for hv, rv in zip(target_h, t_res):\n", + " print(hv, rv)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd2fc482-88a5-4edb-8b9c-77e8496202d8", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "dev", + "language": "python", + "name": "dev" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 } diff --git a/docs/examples/interpolate_hybrid_to_pl.ipynb b/docs/examples/interpolate_hybrid_to_pl.ipynb index 14939f34..43c39660 100644 --- a/docs/examples/interpolate_hybrid_to_pl.ipynb +++ b/docs/examples/interpolate_hybrid_to_pl.ipynb @@ -116,9 +116,9 @@ "from earthkit.meteo.utils.sample import get_sample\n", "\n", "DATA = get_sample(\"vertical_hybrid_data\")\n", - "sp = DATA.p_surf # surface pressure [Pa]\n", - "t = DATA.t # temperature [K]\n", - "q = DATA.q # specific humidity [kg/kg]\n", + "sp = DATA.p_surf # surface pressure [Pa]\n", + "t = DATA.t # temperature [K]\n", + "q = DATA.q # specific humidity [kg/kg]\n", "\n", "sp.shape, t.shape, q.shape" ] @@ -247,13 +247,16 @@ } ], "source": [ - "target_p = [85000., 50000.] # Pa\n", + "target_p = [85000.0, 50000.0] # Pa\n", "\n", "t_res = vertical.interpolate_hybrid_to_pressure_levels(\n", - " t, # data to interpolate\n", + " t, # data to interpolate\n", " target_p,\n", - " A, B, sp, \n", - " interpolation=\"linear\")\n", + " A,\n", + " B,\n", + " sp,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for pv, rv in zip(target_p, t_res):\n", " print(pv, rv)" @@ -314,10 +317,13 @@ "# using the 50 lowest model levels as input\n", "# Please note we still need to use the full A and B arrays\n", "t_res = vertical.interpolate_hybrid_to_pressure_levels(\n", - " t[-50:], # data to interpolate\n", - " target_p, \n", - " A, B, sp,\n", - " interpolation=\"linear\")\n", + " t[-50:], # data to interpolate\n", + " target_p,\n", + " A,\n", + " B,\n", + " sp,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for pv, rv in zip(target_p, t_res):\n", " print(pv, rv)" @@ -414,13 +420,16 @@ } ], "source": [ - "target_p = [95100.] # Pa\n", + "target_p = [95100.0] # Pa\n", "\n", "t_res = vertical.interpolate_hybrid_to_pressure_levels(\n", - " t, # data to interpolate\n", + " t, # data to interpolate\n", " target_p,\n", - " A, B, sp, \n", - " interpolation=\"linear\")\n", + " A,\n", + " B,\n", + " sp,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for pv, rv in zip(target_p, t_res):\n", " print(pv, rv)" @@ -461,15 +470,18 @@ } ], "source": [ - "target_p = [95100.] # Pa\n", + "target_p = [95100.0] # Pa\n", "\n", "t_res = vertical.interpolate_hybrid_to_pressure_levels(\n", - " t, # data to interpolate\n", + " t, # data to interpolate\n", " target_p,\n", - " A, B, sp, \n", - " aux_bottom_p=[ 95178.337944 , 102659.81019512],\n", + " A,\n", + " B,\n", + " sp,\n", + " aux_bottom_p=[95178.337944, 102659.81019512],\n", " aux_bottom_data=[270.0, 293.0],\n", - " interpolation=\"linear\")\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for pv, rv in zip(target_p, t_res):\n", " print(pv, rv)" @@ -582,7 +594,7 @@ } ], "source": [ - "target_p = [85000., 50000.] # Pa\n", + "target_p = [85000.0, 50000.0] # Pa\n", "\n", "t_res = vertical.interpolate_monotonic(t, p, target_p, interpolation=\"linear\")\n", "\n", @@ -686,9 +698,9 @@ "source": [ "# compute pressure on the 50 lowest model levels\n", "# Please note we still need to use the full A and B arrays.\n", - "levels = list(range(138-50, 138)) # 97-137\n", + "levels = list(range(138 - 50, 138)) # 97-137\n", "p_sub = vertical.pressure_on_hybrid_levels(A, B, sp, levels=levels, output=\"full\")\n", - " \n", + "\n", "\n", "# interpolate to pressure levels\n", "t_res = vertical.interpolate_monotonic(t[-50:], p_sub, target_p, interpolation=\"linear\")\n", @@ -746,8 +758,8 @@ } ], "source": [ - "# this pressure is larger than the pressure on the lowest model level in point 1 \n", - "target_p=[95100.] # Pa\n", + "# this pressure is larger than the pressure on the lowest model level in point 1\n", + "target_p = [95100.0] # Pa\n", "\n", "t_res = vertical.interpolate_monotonic(t, p, target_p, interpolation=\"linear\")\n", "\n", @@ -790,12 +802,16 @@ } ], "source": [ - "target_p=[95100.] # Pa\n", + "target_p = [95100.0] # Pa\n", "\n", - "t_res = vertical.interpolate_monotonic(t, p, target_p, \n", - " aux_max_level_coord=[ 95178.337944 , 102659.81019512],\n", - " aux_max_level_data=[270.0, 293.0],\n", - " interpolation=\"linear\")\n", + "t_res = vertical.interpolate_monotonic(\n", + " t,\n", + " p,\n", + " target_p,\n", + " aux_max_level_coord=[95178.337944, 102659.81019512],\n", + " aux_max_level_data=[270.0, 293.0],\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for pv, rv in zip(target_p, t_res):\n", " print(pv, rv)" diff --git a/docs/examples/interpolate_pl_to_hl.ipynb b/docs/examples/interpolate_pl_to_hl.ipynb index 6e2b1a21..b7de77b6 100644 --- a/docs/examples/interpolate_pl_to_hl.ipynb +++ b/docs/examples/interpolate_pl_to_hl.ipynb @@ -101,11 +101,11 @@ "from earthkit.meteo.utils.sample import get_sample\n", "\n", "DATA = get_sample(\"vertical_pl_data\")\n", - "p = DATA.p # pressure [Pa]\n", - "t = DATA.t # temperature [K]\n", - "z = DATA.z # geopotential [m2/s2]\n", - "sp = DATA.p_surf # surface pressure [Pa]\n", - "zs = DATA.z_surf # surface geopotential [m2/s2]\n", + "p = DATA.p # pressure [Pa]\n", + "t = DATA.t # temperature [K]\n", + "z = DATA.z # geopotential [m2/s2]\n", + "sp = DATA.p_surf # surface pressure [Pa]\n", + "zs = DATA.z_surf # surface geopotential [m2/s2]\n", "\n", "p, t.shape" ] @@ -161,15 +161,17 @@ } ], "source": [ - "target_h=[1000., 5000.] # geometric height (m), above the ground\n", + "target_h = [1000.0, 5000.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_pressure_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " z, zs,\n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\", \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " target_h,\n", + " z,\n", + " zs,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -225,17 +227,19 @@ } ], "source": [ - "target_h=[1000., 5000.] # geometric height (m), above the sea\n", + "target_h = [1000.0, 5000.0] # geometric height (m), above the sea\n", "\n", "# zs is not used in interpolate_pressure_to_height_levels()\n", "# when h_reference=\"sea\", so we can pass None for it\n", "t_res = vertical.interpolate_pressure_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " z, None, \n", - " h_type=\"geometric\", \n", - " h_reference=\"sea\", \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " target_h,\n", + " z,\n", + " None,\n", + " h_type=\"geometric\",\n", + " h_reference=\"sea\",\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -291,15 +295,17 @@ } ], "source": [ - "target_h=[1000., 5000.] # geopotential height (m), above the ground\n", + "target_h = [1000.0, 5000.0] # geopotential height (m), above the ground\n", "\n", "t_res = vertical.interpolate_pressure_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " z, zs,\n", - " h_type=\"geopotential\", \n", - " h_reference=\"ground\", \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " target_h,\n", + " z,\n", + " zs,\n", + " h_type=\"geopotential\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -395,15 +401,17 @@ } ], "source": [ - "target_h=[2.] # geometric height (m), above the ground\n", + "target_h = [2.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_pressure_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " z, zs,\n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\", \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " target_h,\n", + " z,\n", + " zs,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -444,17 +452,19 @@ } ], "source": [ - "target_h=[2.] # geometric height (m), above the ground\n", + "target_h = [2.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_pressure_to_height_levels(\n", - " t, # data to interpolate\n", - " target_h, \n", - " z, zs,\n", - " h_type=\"geometric\", \n", - " h_reference=\"ground\", \n", - " aux_bottom_h=0.,\n", - " aux_bottom_data=[304., 306.],\n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " target_h,\n", + " z,\n", + " zs,\n", + " h_type=\"geometric\",\n", + " h_reference=\"ground\",\n", + " aux_bottom_h=0.0,\n", + " aux_bottom_data=[304.0, 306.0],\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -516,7 +526,7 @@ }, "outputs": [], "source": [ - "# geometric height is a non-linear function of geopotential. So we need to \n", + "# geometric height is a non-linear function of geopotential. So we need to\n", "# do the subtraction as follows:\n", "h = vertical.geometric_height_from_geopotential(z) - vertical.geometric_height_from_geopotential(zs)" ] @@ -557,13 +567,14 @@ } ], "source": [ - "target_h=[1000., 5000.] # geometric height (m), above the ground\n", + "target_h = [1000.0, 5000.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_monotonic(\n", - " t, # data to interpolate\n", - " h, \n", - " target_h, \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " h,\n", + " target_h,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -622,13 +633,14 @@ "# compute geometric height above sea level\n", "h_sea = vertical.geometric_height_from_geopotential(z)\n", "\n", - "target_h=[1000., 5000.] # geometric height (m), above the sea\n", + "target_h = [1000.0, 5000.0] # geometric height (m), above the sea\n", "\n", "t_res = vertical.interpolate_monotonic(\n", - " t, # data to interpolate\n", - " h_sea, \n", - " target_h, \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " h_sea,\n", + " target_h,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -683,13 +695,14 @@ } ], "source": [ - "target_h=[2.] # geometric height (m), above the ground\n", + "target_h = [2.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_monotonic(\n", - " t, # data to interpolate\n", + " t, # data to interpolate\n", " h,\n", - " target_h, \n", - " interpolation=\"linear\")\n", + " target_h,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" @@ -730,15 +743,16 @@ } ], "source": [ - "target_h=[2.] # geometric height (m), above the ground\n", + "target_h = [2.0] # geometric height (m), above the ground\n", "\n", "t_res = vertical.interpolate_monotonic(\n", - " t, # data to interpolate\n", + " t, # data to interpolate\n", " h,\n", - " target_h, \n", - " aux_min_level_coord=0.,\n", - " aux_min_level_data=[304., 306.],\n", - " interpolation=\"linear\")\n", + " target_h,\n", + " aux_min_level_coord=0.0,\n", + " aux_min_level_data=[304.0, 306.0],\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hv, rv in zip(target_h, t_res):\n", " print(hv, rv)" diff --git a/docs/examples/interpolate_pl_to_pl.ipynb b/docs/examples/interpolate_pl_to_pl.ipynb index cd537c40..0273c5fa 100644 --- a/docs/examples/interpolate_pl_to_pl.ipynb +++ b/docs/examples/interpolate_pl_to_pl.ipynb @@ -101,8 +101,8 @@ "from earthkit.meteo.utils.sample import get_sample\n", "\n", "DATA = get_sample(\"vertical_pl_data\")\n", - "p = DATA.p # pressure [Pa]\n", - "t = DATA.t # temperature [K]\n", + "p = DATA.p # pressure [Pa]\n", + "t = DATA.t # temperature [K]\n", "\n", "p, t.shape" ] @@ -158,13 +158,14 @@ } ], "source": [ - "target_p=[87000., 60000.] # Pa\n", + "target_p = [87000.0, 60000.0] # Pa\n", "\n", "t_res = vertical.interpolate_monotonic(\n", - " t, # data to interpolate\n", - " p, \n", - " target_p, \n", - " interpolation=\"linear\")\n", + " t, # data to interpolate\n", + " p,\n", + " target_p,\n", + " interpolation=\"linear\",\n", + ")\n", "\n", "for hp, rv in zip(target_p, t_res):\n", " print(hp, rv)" diff --git a/docs/examples/return_period.ipynb b/docs/examples/return_period.ipynb index 441985fc..9cdb3968 100644 --- a/docs/examples/return_period.ipynb +++ b/docs/examples/return_period.ipynb @@ -14,9 +14,10 @@ "outputs": [], "source": [ "import earthkit.data as ekd\n", - "import earthkit.meteo.stats as ekm_stats\n", "import earthkit.plots as ekp\n", "\n", + "import earthkit.meteo.stats as ekm_stats\n", + "\n", "ekd.settings.set(\"cache-policy\", \"user\")" ] }, @@ -41,16 +42,48 @@ " product_type=[\"era5\"],\n", " temporal_aggregation=[\"yearly\"],\n", " period=[\n", - " \"1979\", \"1980\", \"1981\", \"1982\", \"1983\",\n", - " \"1984\", \"1985\", \"1986\", \"1987\", \"1988\",\n", - " \"1989\", \"1990\", \"1991\", \"1992\", \"1993\",\n", - " \"1994\", \"1995\", \"1996\", \"1997\", \"1998\",\n", - " \"1999\", \"2000\", \"2001\", \"2002\", \"2003\",\n", - " \"2004\", \"2005\", \"2006\", \"2007\", \"2008\",\n", - " \"2009\", \"2010\", \"2011\", \"2012\", \"2013\",\n", - " \"2014\", \"2015\", \"2016\", \"2017\", \"2018\",\n", - " \"2019\"\n", - " ]\n", + " \"1979\",\n", + " \"1980\",\n", + " \"1981\",\n", + " \"1982\",\n", + " \"1983\",\n", + " \"1984\",\n", + " \"1985\",\n", + " \"1986\",\n", + " \"1987\",\n", + " \"1988\",\n", + " \"1989\",\n", + " \"1990\",\n", + " \"1991\",\n", + " \"1992\",\n", + " \"1993\",\n", + " \"1994\",\n", + " \"1995\",\n", + " \"1996\",\n", + " \"1997\",\n", + " \"1998\",\n", + " \"1999\",\n", + " \"2000\",\n", + " \"2001\",\n", + " \"2002\",\n", + " \"2003\",\n", + " \"2004\",\n", + " \"2005\",\n", + " \"2006\",\n", + " \"2007\",\n", + " \"2008\",\n", + " \"2009\",\n", + " \"2010\",\n", + " \"2011\",\n", + " \"2012\",\n", + " \"2013\",\n", + " \"2014\",\n", + " \"2015\",\n", + " \"2016\",\n", + " \"2017\",\n", + " \"2018\",\n", + " \"2019\",\n", + " ],\n", ")" ] }, @@ -83,7 +116,7 @@ "metadata": {}, "outputs": [], "source": [ - "threshold = 30.\n", + "threshold = 30.0\n", "\n", "return_period = ekm_stats.value_to_return_period(dist, threshold)" ] diff --git a/docs/examples/seven_weather_regimes.ipynb b/docs/examples/seven_weather_regimes.ipynb index c264c218..0c13801a 100644 --- a/docs/examples/seven_weather_regimes.ipynb +++ b/docs/examples/seven_weather_regimes.ipynb @@ -17,9 +17,9 @@ "metadata": {}, "outputs": [], "source": [ - "import pooch\n", "import numpy as np\n", "import pandas as pd\n", + "import pooch\n", "import xarray as xr\n", "\n", "from earthkit.meteo import regimes" @@ -45,9 +45,10 @@ "files = pooch.retrieve(\n", " url=\"doi:10.5281/zenodo.18154492/wr_data_package_V1.1.zip\",\n", " known_hash=\"dc942ff2a1b3da6dedd3b0b2fadda017fff9e8fc10228ace31c2209a9be7dc62\",\n", - " processor=pooch.Unzip()\n", + " processor=pooch.Unzip(),\n", ")\n", "\n", + "\n", "def get_file(name):\n", " for file in files:\n", " if file.endswith(name):\n", @@ -90,9 +91,10 @@ "metadata": {}, "outputs": [], "source": [ - "mod = 1. / mod_ds[\"normwgt\"].to_series()\n", + "mod = 1.0 / mod_ds[\"normwgt\"].to_series()\n", "mod.index = pd.Index(zip(mod.index.month, mod.index.day, mod.index.hour))\n", "\n", + "\n", "def pattern_normalisation_weight(date):\n", " date = np.asarray(date)\n", " shp = date.shape\n", @@ -143,10 +145,10 @@ " labels=pattern_ds.attrs[\"ClassNames\"].split(),\n", " grid={\n", " \"grid\": [0.5, 0.5],\n", - " \"area\": [90, -80, 30, 40] # 30-90°N, 80°W-40°E\n", + " \"area\": [90, -80, 30, 40], # 30-90°N, 80°W-40°E\n", " },\n", " base_patterns=pattern_ds[\"Z0500_mean\"].values,\n", - " modulator=pattern_normalisation_weight\n", + " modulator=pattern_normalisation_weight,\n", ")" ] }, @@ -740,7 +742,7 @@ "# Reconstruct spatial coordinates\n", "ds_field = ds_field.assign_coords({\n", " \"lat\": np.linspace(ds_field.attrs[\"domymin\"], ds_field.attrs[\"domymax\"], ds_field[\"Z0\"].shape[0]),\n", - " \"lon\": np.linspace(ds_field.attrs[\"domxmin\"], ds_field.attrs[\"domxmax\"], ds_field[\"Z0\"].shape[1])\n", + " \"lon\": np.linspace(ds_field.attrs[\"domxmin\"], ds_field.attrs[\"domxmax\"], ds_field[\"Z0\"].shape[1]),\n", "})\n", "\n", "# Make time coordinate a dimension so it can be used by the pattern generator\n", @@ -824,9 +826,10 @@ " std = f.readline().strip().split()[1:]\n", " return (\n", " xr.DataArray([float(v) for v in mean], coords={\"pattern\": name}),\n", - " xr.DataArray([float(v) for v in std], coords={\"pattern\": name})\n", + " xr.DataArray([float(v) for v in std], coords={\"pattern\": name}),\n", " )\n", - " \n", + "\n", + "\n", "with open(get_file(\"wr_data/WRI_std_params.txt\"), \"r\") as f:\n", " norm_mean, norm_std = read_wri_std_parameters(f)" ] diff --git a/docs/release_notes/include/migrated_hybrid_pressure_at_height_levels.py b/docs/release_notes/include/migrated_hybrid_pressure_at_height_levels.py index 856bcdce..76d82ca4 100644 --- a/docs/release_notes/include/migrated_hybrid_pressure_at_height_levels.py +++ b/docs/release_notes/include/migrated_hybrid_pressure_at_height_levels.py @@ -40,9 +40,7 @@ # Option 2 # compute the pressure on full hybrid levels -p_full, alpha, delta = vertical.pressure_on_hybrid_levels( - A, B, sp, alpha_top="ifs", output=("full", "alpha", "delta") -) +p_full, alpha, delta = vertical.pressure_on_hybrid_levels(A, B, sp, alpha_top="ifs", output=("full", "alpha", "delta")) z = vertical.relative_geopotential_thickness_on_hybrid_levels_from_alpha_delta(t, q, alpha, delta) h = vertical.geopotential_height_from_geopotential(z) diff --git a/pyproject.toml b/pyproject.toml index f4337512..7d0ba4bc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,17 +1,11 @@ [build-system] build-backend = "setuptools.build_meta" -requires = [ "setuptools>=61", "setuptools-scm>=8", "wheel" ] +requires = ["setuptools>=61", "setuptools-scm>=8", "wheel"] [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.10" - classifiers = [ "Development Status :: 2 - Pre-Alpha", "Intended Audience :: Developers", @@ -23,20 +17,25 @@ classifiers = [ "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", "Programming Language :: Python :: Implementation :: CPython", - "Programming Language :: Python :: Implementation :: PyPy", + "Programming Language :: Python :: Implementation :: PyPy" ] -dynamic = [ "version" ] dependencies = [ "deprecation", "earthkit-utils>=0.2", "numpy", - "xarray", + "xarray" ] +description = "Meteorological computations" +dynamic = ["version"] +license = {text = "Apache License Version 2.0"} +name = "earthkit-meteo" +readme = "README.md" +requires-python = ">=3.10" optional-dependencies.all = [ - "earthkit-meteo[samples,scores]", + "earthkit-meteo[samples,scores]" ] optional-dependencies.dev = [ - "earthkit-meteo[all,docs,test]", + "earthkit-meteo[all,docs,test]" ] optional-dependencies.docs = [ "nbsphinx", @@ -46,43 +45,71 @@ optional-dependencies.docs = [ "sphinx-copybutton", "sphinx-issues", "sphinx-rtd-theme", - "sphinx-tabs", + "sphinx-tabs" ] optional-dependencies.gpu = [ "cupy", - "torch", + "torch" ] optional-dependencies.samples = [ - "requests", + "requests" ] optional-dependencies.scores = [ - "scores", + "scores" ] optional-dependencies.test = [ "earthkit-meteo[samples,scores]", "netcdf4", "pytest", "pytest-cov", - "xarray", + "xarray" ] 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.Issues = "https://github.com/ecmwf/earthkit-meteo/issues" urls.Repository = "https://github.com/ecmwf/earthkit-meteo/" -[tool.setuptools.packages.find] -include = [ "earthkit.meteo" ] -where = [ "src/" ] - -[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.ruff] +line-length = 120 +preview = true + +[tool.ruff.lint] +ignore = [ + "D1", # pydocstyle: Missing Docstrings + "D107", # pydocstyle: numpy convention + "D203", + "D205", + "D212", + "D213", + "D401", + "D402", + "D413", + "D415", + "D416", + "D417", + "D420" +] +select = [ + "F", # pyflakes + "E", # pycodestyle + "W", # pycodestyle warnings + "I", # isort + "D" # pydocstyle +] + +[tool.ruff.lint.per-file-ignores] +"src/earthkit/meteo/score/xarray/ensemble.py" = ["E501"] + +[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/__init__.py b/src/earthkit/meteo/extreme/array/__init__.py index dcf22428..4629fcac 100644 --- a/src/earthkit/meteo/extreme/array/__init__.py +++ b/src/earthkit/meteo/extreme/array/__init__.py @@ -7,9 +7,7 @@ # nor does it submit to any jurisdiction. # -""" -Extreme index functions operating on numpy arrays. -""" +"""Extreme index functions operating on numpy arrays.""" from .cpf import * # noqa from .efi import * # noqa diff --git a/src/earthkit/meteo/extreme/array/cpf.py b/src/earthkit/meteo/extreme/array/cpf.py index e7d17c9e..58e6c5a8 100644 --- a/src/earthkit/meteo/extreme/array/cpf.py +++ b/src/earthkit/meteo/extreme/array/cpf.py @@ -100,7 +100,7 @@ def cpf( symmetric=False, from_zero=False, ): - """Compute Crossing Point Forecast (CPF) + """Compute Crossing Point Forecast (CPF). WARNING: this code is experimental, use at your own risk! diff --git a/src/earthkit/meteo/extreme/array/efi.py b/src/earthkit/meteo/extreme/array/efi.py index b533e09e..8d65d74d 100644 --- a/src/earthkit/meteo/extreme/array/efi.py +++ b/src/earthkit/meteo/extreme/array/efi.py @@ -14,7 +14,7 @@ def efi(clim, ens, eps=-0.1): - """Compute Extreme Forecast Index (EFI) + """Compute Extreme Forecast Index (EFI). Parameters ---------- @@ -30,7 +30,6 @@ def efi(clim, ens, eps=-0.1): array-like (npoints) EFI values """ - xp = array_namespace(clim, ens) clim = xp.asarray(clim) ens = xp.asarray(ens) diff --git a/src/earthkit/meteo/extreme/array/sot.py b/src/earthkit/meteo/extreme/array/sot.py index da5e5159..72008807 100644 --- a/src/earthkit/meteo/extreme/array/sot.py +++ b/src/earthkit/meteo/extreme/array/sot.py @@ -11,7 +11,7 @@ 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 + """Compute basic Shift of Tails (SOT) using already computed percentiles. Parameters ---------- @@ -51,7 +51,7 @@ def sot_func(qc_tail, qc, qf, eps=-1e-4, lower_bound=-10, upper_bound=10): def sot(clim, ens, perc, eps=-1e4): """Compute Shift of Tails (SOT) from climatology percentiles (sorted) - and ensemble forecast (not sorted) + and ensemble forecast (not sorted). Parameters ---------- @@ -78,9 +78,7 @@ def sot(clim, ens, perc, eps=-1e4): 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) - ) + 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 @@ -94,9 +92,7 @@ def sot(clim, ens, perc, eps=-1e4): 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" - ) + 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) @@ -106,7 +102,7 @@ def sot(clim, ens, perc, eps=-1e4): def sot_unsorted(clim, ens, perc, eps=-1e4): """Compute Shift of Tails (SOT) from climatology percentiles (sorted) - and ensemble forecast (not sorted) + and ensemble forecast (not sorted). Parameters ---------- @@ -133,9 +129,7 @@ def sot_unsorted(clim, ens, perc, eps=-1e4): 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) - ) + raise Exception("Climatology array should contain 101 percentiles, it has {} values".format(clim.shape)) if eps > 0: ens = xp.where(ens < eps, 0.0, ens) @@ -148,9 +142,7 @@ def sot_unsorted(clim, ens, perc, eps=-1e4): 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" - ) + raise Exception("Percentile value to be computed cannot be 50 for sot, has to be in the upper or lower half") qc_tail = xp.percentile(clim, q=perc_tail, axis=0) sot = sot_func(qc_tail, qc, qf, eps=eps) diff --git a/src/earthkit/meteo/regimes/__init__.py b/src/earthkit/meteo/regimes/__init__.py index 39b7f81a..edf4e8e4 100644 --- a/src/earthkit/meteo/regimes/__init__.py +++ b/src/earthkit/meteo/regimes/__init__.py @@ -27,10 +27,7 @@ """ from . import array -from .patterns import ConstantPatterns -from .patterns import ModulatedPatterns -from .patterns import Patterns +from .patterns import ConstantPatterns, ModulatedPatterns, Patterns # Offer xarray implementations at high-level (TODO: support fieldlist) -from .xarray import project -from .xarray import regime_index +from .xarray import project, regime_index diff --git a/src/earthkit/meteo/regimes/array/index.py b/src/earthkit/meteo/regimes/array/index.py index d7bb0a7c..21cd0de0 100644 --- a/src/earthkit/meteo/regimes/array/index.py +++ b/src/earthkit/meteo/regimes/array/index.py @@ -36,9 +36,7 @@ def project(field, patterns, weights, **patterns_extra_coords): """ ndim_field = len(patterns.shape) if field.shape[-ndim_field:] != patterns.shape: - raise ValueError( - f"shape of input fields {field.shape} incompatible with shape of patterns {patterns.shape}" - ) + raise ValueError(f"shape of input fields {field.shape} incompatible with shape of patterns {patterns.shape}") if weights is None: # TODO generate area-based weights from grid of patterns with earthkit-geo diff --git a/src/earthkit/meteo/regimes/patterns.py b/src/earthkit/meteo/regimes/patterns.py index fe4c25fd..aa86312b 100644 --- a/src/earthkit/meteo/regimes/patterns.py +++ b/src/earthkit/meteo/regimes/patterns.py @@ -96,13 +96,9 @@ def _patterns_iterxr(self, reference_da, patterns_extra_coords): dims = [*extra_dims, *reference_da.dims[-self.ndim :]] coords = {dim: reference_da.coords[dim] for dim in dims} # Cartesian product of coordinates for patterns generator - extra_coords_arrs = dict( - zip(extra_dims, self.xp.meshgrid(*(coords[dim] for dim in extra_dims), indexing="ij")) - ) + extra_coords_arrs = dict(zip(extra_dims, self.xp.meshgrid(*(coords[dim] for dim in extra_dims), indexing="ij"))) # Rearrange to match provided kwarg-coord mapping - extra_coords = { - kwarg: extra_coords_arrs[patterns_extra_coords[kwarg]] for kwarg in patterns_extra_coords - } + extra_coords = {kwarg: extra_coords_arrs[patterns_extra_coords[kwarg]] for kwarg in patterns_extra_coords} # Delegate the pattern generation and package the patterns as DataArrays for name, patterns in self.patterns(**extra_coords).items(): yield name, xr.DataArray(patterns, coords=coords, dims=dims) diff --git a/src/earthkit/meteo/score/array/deterministic.py b/src/earthkit/meteo/score/array/deterministic.py index a94d9f89..42baf5da 100644 --- a/src/earthkit/meteo/score/array/deterministic.py +++ b/src/earthkit/meteo/score/array/deterministic.py @@ -11,7 +11,7 @@ def pearson_correlation(x, y, axis=0): - """Compute pearson correlation coefficient + """Compute pearson correlation coefficient. Parameters ---------- diff --git a/src/earthkit/meteo/score/xarray/deterministic.py b/src/earthkit/meteo/score/xarray/deterministic.py index a4b56eb7..f92dd881 100644 --- a/src/earthkit/meteo/score/xarray/deterministic.py +++ b/src/earthkit/meteo/score/xarray/deterministic.py @@ -1,5 +1,4 @@ -from typing import Literal -from typing import TypeVar +from typing import Literal, TypeVar import numpy as np import xarray as xr @@ -44,7 +43,9 @@ def error( .. seealso:: - This function leverages the `scores.continuous.additive_bias `_ function. + This function leverages the + `scores.continuous.additive_bias + `_ function. Parameters ---------- @@ -106,7 +107,9 @@ def mean_error( .. seealso:: - This function leverages the `scores.continuous.additive_bias `_ function. + This function leverages the + `scores.continuous.additive_bias + `_ function. Parameters ---------- @@ -160,7 +163,8 @@ def abs_error( .. seealso:: - This function leverages the `scores.continuous.mae `_ function. + This function leverages the + `scores.continuous.mae `_ function. Parameters ---------- @@ -221,7 +225,8 @@ def mean_abs_error( .. seealso:: - This function leverages the `scores.continuous.mae `_ function. + This function leverages the + `scores.continuous.mae `_ function. Parameters ---------- @@ -271,7 +276,8 @@ def squared_error( .. seealso:: - This function leverages the `scores.continuous.mse `_ function. + This function leverages the + `scores.continuous.mse `_ function. Parameters ---------- @@ -332,7 +338,8 @@ def mean_squared_error( .. seealso:: - This function leverages the `scores.continuous.mse `_ function. + This function leverages the + `scores.continuous.mse `_ function. Parameters ---------- @@ -389,7 +396,8 @@ def root_mean_squared_error( .. seealso:: - This function leverages the `scores.continuous.mse `_ function. + This function leverages the + `scores.continuous.mse `_ function. Parameters ---------- @@ -670,11 +678,13 @@ def kge( - :math:`\rho` = Pearson's correlation coefficient between observed and forecast values. - :math:`f` and :math:`o` are forecast and observed values, respectively - :math:`\mu_f` and :math:`\mu_o` are the means of forecast and observed values, respectively - - :math:`\sigma_f` and :math:`\sigma_o` are the standard deviations of forecast and observed values, respectively + - :math:`\sigma_f` and :math:`\sigma_o` are the standard deviations of forecast and observed values, + respectively .. seealso:: - This function leverages the `scores.continuous.kge `_ function. + This function leverages the + `scores.continuous.kge `_ function. Parameters ---------- @@ -685,9 +695,11 @@ def kge( over : str or list of str The dimension(s) over which to compute the kge. method : str, optional - The method to compute the variability term :math:`\alpha`. Can be either "original" or "modified". Default is "modified". + The method to compute the variability term :math:`\alpha`. + Can be either "original" or "modified". Default is "modified". return_components : bool, optional - Whether to return the individual components (:math:`\rho`, :math:`\alpha` (or :math:`\gamma`), :math:`\beta`) along with the KGE value. Default is False. + Whether to return the individual components (:math:`\rho`, :math:`\alpha` (or :math:`\gamma`), + :math:`\beta`) along with the KGE value. Default is False. Returns ------- @@ -706,9 +718,7 @@ def kge( if method == "modified": kge["alpha"] = kge["alpha"] / kge["beta"] - kge["kge"] = 1 - ( - ((kge["rho"] - 1) ** 2 + ((kge["alpha"] - 1) ** 2) + ((kge["beta"] - 1) ** 2)) ** 0.5 - ) + kge["kge"] = 1 - (((kge["rho"] - 1) ** 2 + ((kge["alpha"] - 1) ** 2) + ((kge["beta"] - 1) ** 2)) ** 0.5) kge = kge.rename({"alpha": "gamma"}) diff --git a/src/earthkit/meteo/score/xarray/ensemble.py b/src/earthkit/meteo/score/xarray/ensemble.py index 75b940d7..7ba57063 100644 --- a/src/earthkit/meteo/score/xarray/ensemble.py +++ b/src/earthkit/meteo/score/xarray/ensemble.py @@ -1,5 +1,4 @@ -from typing import Literal -from typing import TypeVar +from typing import Literal, TypeVar import numpy as np import xarray as xr @@ -50,7 +49,6 @@ def spread(fcst: T, over: str | list[str], reference: T | None = None) -> T: xarray object The spread of the forecast compared to the reference. """ - # TODO: this could call the rmse function if reference is None: reference = fcst.mean(dim=over) @@ -129,10 +127,15 @@ def crps_from_gaussian(fcst: xr.Dataset, obs: xr.DataArray) -> xr.DataArray: - :math:`\mathcal{N}(\mu, \sigma^2)` is the probabilistic (Gaussian) forecast, - :math:`o` are the observations, - - :math:`\phi\left( (o - \mu)/\sigma \right)` denotes the probability density function of the normal distribution with mean 0 and variance 1 evaluated at the normalised prediction error, :math:`(o - \mu)/\sigma`, - - :math:`\Phi\left( (o - \mu)/\sigma \right)` denotes the cumulative distribution function of the normal distribution with mean 0 and variance 1 evaluated at the normalised prediction error, :math:`(o - \mu)/\sigma`. + - :math:`\phi\left( (o - \mu)/\sigma \right)` denotes the probability density function of the normal + distribution with mean 0 and variance 1 evaluated at the normalised prediction error, + :math:`(o - \mu)/\sigma`, + - :math:`\Phi\left( (o - \mu)/\sigma \right)` denotes the cumulative distribution function of the normal + distribution with mean 0 and variance 1 evaluated at the normalised prediction error, + :math:`(o - \mu)/\sigma`. - Reference: Gneiting, Tilmann, et al. "Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation." Monthly weather review 133.5 (2005): 1098-1118. + Reference: Gneiting, Tilmann, et al. "Calibrated probabilistic forecasting using ensemble model output + statistics and minimum CRPS estimation." Monthly weather review 133.5 (2005): 1098-1118. Parameters ---------- @@ -149,9 +152,7 @@ def crps_from_gaussian(fcst: xr.Dataset, obs: xr.DataArray) -> xr.DataArray: if not isinstance(fcst, xr.Dataset): raise TypeError(f"Expected fcst to be an xarray.Dataset object, got {type(fcst)}") if not {"mean", "stdev"}.issubset(fcst.data_vars): - raise ValueError( - f"Expected fcst to have 'mean' and 'stdev' data variables, got {list(fcst.data_vars)}" - ) + raise ValueError(f"Expected fcst to have 'mean' and 'stdev' data variables, got {list(fcst.data_vars)}") if not isinstance(obs, xr.DataArray): raise TypeError(f"Expected obs to be an xarray.DataArray object, got {type(obs)}") @@ -161,9 +162,7 @@ def crps_from_gaussian(fcst: xr.Dataset, obs: xr.DataArray) -> xr.DataArray: c2 = np.sqrt(2.0 / np.pi) za = (obs - fcst["mean"]) / fcst["stdev"] return fcst["stdev"] * ( - (2.0 * scipy.stats.norm().cdf(za.values) - 1.0) * za - + c2 * np.exp(-(za**2) / 2.0) - - 1.0 / np.sqrt(np.pi) + (2.0 * scipy.stats.norm().cdf(za.values) - 1.0) * za + c2 * np.exp(-(za**2) / 2.0) - 1.0 / np.sqrt(np.pi) ) @@ -195,11 +194,13 @@ def crps_from_ensemble( - :math:`o` are the observations, - :math:`K=M^2` for the 'ecdf' method and :math:`M(M-1)` for the 'fair' method, - With `return_components=True`, this function returns an ``xr.Dataset`` with variables for the decompositions defined below. + With `return_components=True`, this function returns an ``xr.Dataset`` + with variables for the decompositions defined below. If the `decomposition_method="underover"`, the ``xr.Dataset`` variables values are - ``underforecast_penalty``, ``overforecast_penalty``, ``spread`` and either ``fcrps`` if `method="fair"` or ``crps`` if `method="ecdf"` (ordering is not - guaranteed and might differ). The overall CRPS is given by + ``underforecast_penalty``, ``overforecast_penalty``, ``spread`` and either ``fcrps`` if + `method="fair"` or ``crps`` if `method="ecdf"` (ordering is not guaranteed and might differ). + The overall CRPS is given by ``underforecast_penalty + overforecast_penalty - spread``. .. math:: @@ -217,10 +218,13 @@ def crps_from_ensemble( S(f, f) &= \frac{1}{2K} \sum_{i=1}^{M} \sum_{j=1}^{M} |f_i - f_j| \quad& \text{(forecast spread term)} \end{align*} - If the decomposition method is `decomposition_method="hersbach"`, the ``xr.Dataset`` variables values are - ``alpha``, ``beta``, ``crps`` and additionally also ``fcrps`` if `method="fair"` (ordering is not guaranteed and might differ). + If the decomposition method is `decomposition_method="hersbach"`, the ``xr.Dataset`` + variables values are ``alpha``, ``beta``, ``crps`` and additionally also ``fcrps`` + if `method="fair"` (ordering is not guaranteed and might differ). - We denote by :math:`x_1 \le x_2 \le \dots \le x_M` the members of the ensemble forecast :math:`f` after sorting. The unfair CRPS decomposition for `decomposition_method="hersbach"` is then given by + We denote by :math:`x_1 \le x_2 \le \dots \le x_M` the members of the ensemble forecast + :math:`f` after sorting. + The unfair CRPS decomposition for `decomposition_method="hersbach"` is then given by .. math:: :nowrap: @@ -319,38 +323,33 @@ def crps_from_ensemble( include_components=return_components, ) if return_components: - return scores_xr.to_dataset(dim="component").rename( - {"total": "crps" if method == "ecdf" else "fcrps"} - ) + return scores_xr.to_dataset(dim="component").rename({"total": "crps" if method == "ecdf" else "fcrps"}) else: return scores_xr else: valid_mask, alpha, beta, crps, fcrps = _crps_from_ensemble_hersbach(fcst, obs, over) if return_components: if method == "fair": - return xr.Dataset( - { - "alpha": alpha.where(valid_mask), - "beta": beta.where(valid_mask), - "crps": crps.where(valid_mask), - "fcrps": fcrps.where(valid_mask), - } - ) + return xr.Dataset({ + "alpha": alpha.where(valid_mask), + "beta": beta.where(valid_mask), + "crps": crps.where(valid_mask), + "fcrps": fcrps.where(valid_mask), + }) else: - return xr.Dataset( - { - "alpha": alpha.where(valid_mask), - "beta": beta.where(valid_mask), - "crps": crps.where(valid_mask), - } - ) + return xr.Dataset({ + "alpha": alpha.where(valid_mask), + "beta": beta.where(valid_mask), + "crps": crps.where(valid_mask), + }) else: return fcrps.where(valid_mask) if method == "fair" else crps.where(valid_mask) # TODO: does this work when over is a list of dimensions? -# TODO: decide on the nan distribution strategy and make sure it's consistent with other functions (e.g. crps_from_ensemble) +# TODO: decide on the nan distribution strategy and make sure it's consistent with +# other functions (e.g. crps_from_ensemble) def _crps_from_ensemble_hersbach( fcst: T, obs: T, @@ -361,9 +360,7 @@ def _crps_from_ensemble_hersbach( if components_coords is None: components_coords = np.arange(1, ens_size + 2) else: - assert ( - len(components_coords) == ens_size + 1 - ), "component_coords must have the length of ensemble size + 1" + assert len(components_coords) == ens_size + 1, "component_coords must have the length of ensemble size + 1" # sort forecast values along the ensemble dimension fcst_sorted = _sorted_ensemble(fcst, over) alpha = xr.concat([xr.zeros_like(fcst_sorted[{over: 0}])] * (ens_size + 1), dim=over) @@ -375,7 +372,8 @@ def _crps_from_ensemble_hersbach( beta[{over: 0}] = (fcst_sorted[{over: 0}] - obs).where(obs_below_ens, 0.0) rhs = ( - fcst_sorted.diff(dim=over) + fcst_sorted + .diff(dim=over) .where( fcst_sorted[{over: slice(1, None)}] <= obs, -fcst_sorted[{over: slice(None, -1)}] + obs, @@ -386,7 +384,8 @@ def _crps_from_ensemble_hersbach( alpha[{over: slice(1, -1)}] = rhs rhs = ( - fcst_sorted.diff(dim=over) + fcst_sorted + .diff(dim=over) .where( fcst_sorted[{over: slice(None, -1)}] > obs, fcst_sorted[{over: slice(1, None)}] - obs, @@ -482,7 +481,6 @@ def crps_from_cdf( xarray.DataArray or xarray.Dataset The CRPS of the CDF compared to the observations. """ - scores = _import_scores_or_prompt_install() reduce_dim = [over] if return_components: diff --git a/src/earthkit/meteo/solar/array/__init__.py b/src/earthkit/meteo/solar/array/__init__.py index da2e210b..9f5a1e11 100644 --- a/src/earthkit/meteo/solar/array/__init__.py +++ b/src/earthkit/meteo/solar/array/__init__.py @@ -7,8 +7,6 @@ # nor does it submit to any jurisdiction. # -""" -Solar computation functions operating on numpy arrays. -""" +"""Solar computation functions operating on numpy arrays.""" from .solar import * # noqa diff --git a/src/earthkit/meteo/solar/array/solar.py b/src/earthkit/meteo/solar/array/solar.py index c130217c..298de6c2 100644 --- a/src/earthkit/meteo/solar/array/solar.py +++ b/src/earthkit/meteo/solar/array/solar.py @@ -125,22 +125,18 @@ def _integrate( elif integration_order == 4: # slower, more accurate (4 points) _C1 = xp.sqrt(xp.asarray(6.0 / 5.0)) _C2 = xp.sqrt(xp.asarray(30)) - E = xp.asarray( - [ - -xp.sqrt(3.0 / 7.0 + 2.0 / 7.0 * _C1), - -xp.sqrt(3.0 / 7.0 - 2.0 / 7.0 * _C1), - xp.sqrt(3.0 / 7.0 - 2.0 / 7.0 * _C1), - xp.sqrt(3.0 / 7.0 + 2.0 / 7.0 * _C1), - ] - ) - W = xp.asarray( - [ - (18 - _C2) / 36, - (18 + _C2) / 36, - (18 + _C2) / 36, - (18 - _C2) / 36, - ] - ) + E = xp.asarray([ + -xp.sqrt(3.0 / 7.0 + 2.0 / 7.0 * _C1), + -xp.sqrt(3.0 / 7.0 - 2.0 / 7.0 * _C1), + xp.sqrt(3.0 / 7.0 - 2.0 / 7.0 * _C1), + xp.sqrt(3.0 / 7.0 + 2.0 / 7.0 * _C1), + ]) + W = xp.asarray([ + (18 - _C2) / 36, + (18 + _C2) / 36, + (18 + _C2) / 36, + (18 - _C2) / 36, + ]) else: raise ValueError("Invalid integration order %d", integration_order) diff --git a/src/earthkit/meteo/stats/array/__init__.py b/src/earthkit/meteo/stats/array/__init__.py index 9ba6d05d..45a3ae05 100644 --- a/src/earthkit/meteo/stats/array/__init__.py +++ b/src/earthkit/meteo/stats/array/__init__.py @@ -7,9 +7,7 @@ # nor does it submit to any jurisdiction. # -""" -Statistical functions operating on numpy arrays. -""" +"""Statistical functions operating on numpy arrays.""" from .extreme_values import * # noqa from .numpy_extended import * # noqa diff --git a/src/earthkit/meteo/stats/array/quantiles.py b/src/earthkit/meteo/stats/array/quantiles.py index 3e56f14e..d5351be5 100644 --- a/src/earthkit/meteo/stats/array/quantiles.py +++ b/src/earthkit/meteo/stats/array/quantiles.py @@ -7,9 +7,7 @@ # nor does it submit to any jurisdiction. # -from typing import Iterable -from typing import List -from typing import Union +from typing import Iterable, List, Union import numpy as np from earthkit.utils.array import array_namespace @@ -21,7 +19,7 @@ def iter_quantiles( axis: int = 0, method: str = "sort", ) -> Iterable[np.ndarray]: - """Iterate over the quantiles of a large array + """Iterate over the quantiles of a large array. Parameters ---------- @@ -43,7 +41,6 @@ def iter_quantiles( 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'") diff --git a/src/earthkit/meteo/thermo/__init__.py b/src/earthkit/meteo/thermo/__init__.py index 270a8409..bc85986a 100644 --- a/src/earthkit/meteo/thermo/__init__.py +++ b/src/earthkit/meteo/thermo/__init__.py @@ -15,5 +15,4 @@ planned to work with objects like *earthkit.data FieldLists* or *xarray DataSets*. """ - from .thermo import * # noqa diff --git a/src/earthkit/meteo/thermo/array/__init__.py b/src/earthkit/meteo/thermo/array/__init__.py index 5507b1b5..07d568c4 100644 --- a/src/earthkit/meteo/thermo/array/__init__.py +++ b/src/earthkit/meteo/thermo/array/__init__.py @@ -7,8 +7,6 @@ # nor does it submit to any jurisdiction. # -""" -Thermodynamic functions operating on numpy arrays. -""" +"""Thermodynamic functions operating on numpy arrays.""" from .thermo import * # noqa diff --git a/src/earthkit/meteo/thermo/array/thermo.py b/src/earthkit/meteo/thermo/array/thermo.py index 1adfab44..d72a4424 100644 --- a/src/earthkit/meteo/thermo/array/thermo.py +++ b/src/earthkit/meteo/thermo/array/thermo.py @@ -1184,10 +1184,7 @@ def _G_sat(self, ths, scale=1.0): 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 - ) + return -self.K0 * ths.qs / (ths.t**2) + self.K0 * saturation_specific_humidity_slope(ths.t, ths.p) / ths.t def _f(self, ths): xp = ths.ns @@ -1225,10 +1222,7 @@ def _G_sat(self, ths, scale=1.0): def _d_G_sat(self, ths): xp = ths.ns - return ( - -self.K0 * ths.ws / xp.square(ths.t) - + self.K0 * saturation_mixing_ratio_slope(ths.t, ths.p) / ths.t - ) + return -self.K0 * ths.ws / xp.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}") @@ -1236,9 +1230,7 @@ def _f(self, ths): # print(f" exp={self._G_sat(ths, scale=-self.c_lambda)}") xp = ths.ns return ( - ths.c_tw - * xp.pow(ths.p / constants.p0, self.K3 * ths.ws) - * xp.exp(self._G_sat(ths, scale=-self.c_lambda)) + ths.c_tw * xp.pow(ths.p / constants.p0, self.K3 * ths.ws) * xp.exp(self._G_sat(ths, scale=-self.c_lambda)) ) def _d_lnf(self, ths): @@ -1297,9 +1289,9 @@ def _G_sat(self, ths, scale=1.0): def _d_G_sat(self, ths): xp = ths.ns # print(f" d_ws={saturation_mixing_ratio_slope(ths.t, ths.p)}") - return -self.K0 * (ths.ws + self.K2 * xp.square(ths.ws)) / (xp.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 * xp.square(ths.ws)) / (xp.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): xp = ths.ns @@ -1470,7 +1462,7 @@ def saturation_ept(t, p, method="ifs"): def temperature_on_moist_adiabat(ept, p, ept_method="ifs", t_method="bisect"): - r"""Compute the temperature on a moist adiabat (pseudoadiabat) + r"""Compute the temperature on a moist adiabat (pseudoadiabat). Parameters ---------- diff --git a/src/earthkit/meteo/utils/testing.py b/src/earthkit/meteo/utils/testing.py index d4a79495..aba685bf 100644 --- a/src/earthkit/meteo/utils/testing.py +++ b/src/earthkit/meteo/utils/testing.py @@ -78,7 +78,7 @@ def read_test_data_file(path): def save_test_data_reference(file_name, data): - """Helper function to save test reference data into csv""" + """Helper function to save test reference data into csv.""" import numpy as np np.savetxt( diff --git a/src/earthkit/meteo/vertical/__init__.py b/src/earthkit/meteo/vertical/__init__.py index 5552f228..2b96f461 100644 --- a/src/earthkit/meteo/vertical/__init__.py +++ b/src/earthkit/meteo/vertical/__init__.py @@ -15,5 +15,4 @@ 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 index f0939320..ae86a93f 100644 --- a/src/earthkit/meteo/vertical/array/__init__.py +++ b/src/earthkit/meteo/vertical/array/__init__.py @@ -7,9 +7,7 @@ # nor does it submit to any jurisdiction. # -""" -Vertical computation functions operating on numpy arrays. -""" +"""Vertical computation functions operating on numpy arrays.""" from .hybrid import * # noqa from .vertical import * # noqa diff --git a/src/earthkit/meteo/vertical/array/hybrid.py b/src/earthkit/meteo/vertical/array/hybrid.py index 18387b1e..29c49344 100644 --- a/src/earthkit/meteo/vertical/array/hybrid.py +++ b/src/earthkit/meteo/vertical/array/hybrid.py @@ -7,8 +7,7 @@ # nor does it submit to any jurisdiction. # -from typing import Any -from typing import Tuple +from typing import Any, Tuple import numpy as np from numpy.typing import NDArray @@ -41,7 +40,7 @@ def hybrid_level_parameters(n_levels: int, model: str = "ifs") -> Tuple[NDArray[ r"""Get the A and B parameters of hybrid levels for a given configuration. Parameters - ----------- + ---------- n_levels: int Number of (full) hybrid levels. Currently, only ``n_levels`` 91 and 137 are supported. model : str @@ -95,8 +94,6 @@ def hybrid_level_parameters(n_levels: int, model: str = "ifs") -> Tuple[NDArray[ c = _CONF[model][n_levels] return c["A"], c["B"] else: - raise ValueError( - f"Hybrid level parameters not available for {n_levels} levels in model '{model}'." - ) + raise ValueError(f"Hybrid level parameters not available for {n_levels} levels in model '{model}'.") else: raise ValueError(f"Model '{model}' not recognized for hybrid level parameters.") diff --git a/src/earthkit/meteo/vertical/array/monotonic.py b/src/earthkit/meteo/vertical/array/monotonic.py index 463549cb..90be74f7 100644 --- a/src/earthkit/meteo/vertical/array/monotonic.py +++ b/src/earthkit/meteo/vertical/array/monotonic.py @@ -96,7 +96,8 @@ def __call__( if data.shape[0] != coord.shape[0]: raise ValueError( - f"The first dimension of data and that of coord must match! {data.shape=} {coord.shape=} {data.shape[0]} != {coord.shape[0]}" + f"The first dimension of data and that of coord must match! " + f"{data.shape=} {coord.shape=} {data.shape[0]} != {coord.shape[0]}" ) self.data_is_scalar = data[0].ndim == 0 @@ -107,18 +108,15 @@ def __call__( if same_shape: if self.data_is_scalar and not self.target_is_scalar: raise ValueError("If values and p have the same shape, they cannot both be scalars.") - if ( - not self.data_is_scalar - and not self.target_is_scalar - and data.shape[1:] != target_coord.shape[1:] - ): + if not self.data_is_scalar and not self.target_is_scalar and data.shape[1:] != target_coord.shape[1:]: raise ValueError( "When values and target_p have different shapes, target_p must be a scalar or a 1D array." ) if not same_shape and coord.ndim != 1: raise ValueError( - f"When values and p have different shapes, p must be a scalar or a 1D array. {data.shape=} {coord.shape=} {coord.ndim}" + f"When values and p have different shapes, p must be a scalar or a 1D array. " + f"{data.shape=} {coord.shape=} {coord.ndim}" ) # initialize the output array @@ -135,7 +133,8 @@ def __call__( assert self.coord_is_scalar # print(f"scalar_info.target: {scalar_info.target}") # if scalar_info.target: - # return _to_level_1(data, coord, nlev, target_coord, interpolation, scalar_info, xp, res, aux_bottom, aux_top) + # return _to_level_1(data, coord, nlev, target_coord, interpolation, + # scalar_info, xp, res, aux_bottom, aux_top) # else: # coord = xp.broadcast_to(coord, (nlev,) + data.shape[1:]).T @@ -187,7 +186,6 @@ def compute(self, res): # vertical position in the atmosphere. Of course, if the coordinate is pressure these # two definitions coincide. for target_idx, tc in enumerate(self.target_coord): - # find the level below the target idx_bottom = (self.coord > tc).sum(0) idx_bottom = xp.atleast_1d(idx_bottom) diff --git a/src/earthkit/meteo/vertical/array/vertical.py b/src/earthkit/meteo/vertical/array/vertical.py index b270df51..a130ca46 100644 --- a/src/earthkit/meteo/vertical/array/vertical.py +++ b/src/earthkit/meteo/vertical/array/vertical.py @@ -7,15 +7,12 @@ # nor does it submit to any jurisdiction. -from typing import Any -from typing import Tuple -from typing import Union +from typing import Any, Tuple, Union import deprecation import numpy as np from earthkit.utils.array import array_namespace -from numpy.typing import ArrayLike -from numpy.typing import NDArray +from numpy.typing import ArrayLike, NDArray from earthkit.meteo import constants @@ -40,7 +37,8 @@ def pressure_at_model_levels( sp : number or ndarray Surface pressure (Pa) alpha_top : str, optional - Option to initialise alpha on the top of the model atmosphere (first half-level in vertical coordinate system). The possible values are: + Option to initialise alpha on the top of the model atmosphere + (first half-level in vertical coordinate system). The possible values are: - "ifs": alpha is set to log(2). See [IFS-CY47R3-Dynamics]_ (page 7) for details. - "arpege": alpha is set to 1.0 @@ -57,6 +55,11 @@ def pressure_at_model_levels( Alpha at full-levels + See Also + -------- + pressure_at_height_levels + relative_geopotential_thickness + Notes ----- ``A`` and ``B`` must contain the same model half-levels in ascending order with @@ -83,11 +86,6 @@ def pressure_at_model_levels( - :math:`A_{k+1/2}` and :math:`B_{k+1/2}` are the A- and B-coefficients defining the model levels. - See also - -------- - pressure_at_height_levels - relative_geopotential_thickness - """ # constants PRESSURE_TOA = 0.1 # safety when highest pressure level = 0.0 @@ -124,35 +122,25 @@ def pressure_at_model_levels( # calculate alpha alpha = np.zeros(new_shape_full) - alpha[1:, ...] = ( - 1.0 - p_half_level[1:-1, ...] / (p_half_level[2:, ...] - p_half_level[1:-1, ...]) * delta[1:, ...] - ) + alpha[1:, ...] = 1.0 - p_half_level[1:-1, ...] / (p_half_level[2:, ...] - p_half_level[1:-1, ...]) * delta[1:, ...] # pressure at highest half-level <= 0.1 if np.any(p_half_level[0, ...] <= PRESSURE_TOA): alpha[0, ...] = alpha_top # pressure at highest half-level > 0.1 else: - alpha[0, ...] = ( - 1.0 - p_half_level[0, ...] / (p_half_level[1, ...] - p_half_level[0, ...]) * delta[0, ...] - ) + alpha[0, ...] = 1.0 - p_half_level[0, ...] / (p_half_level[1, ...] - p_half_level[0, ...]) * delta[0, ...] # calculate pressure on model full-levels # TODO: is there a faster way to calculate the averages? # TODO: introduce option to calculate full-levels in more complicated way - p_full_level = np.apply_along_axis( - lambda m: np.convolve(m, np.ones(2) / 2, mode="valid"), axis=0, arr=p_half_level - ) + p_full_level = np.apply_along_axis(lambda m: np.convolve(m, np.ones(2) / 2, mode="valid"), axis=0, arr=p_half_level) return p_full_level, p_half_level, delta, alpha -@deprecation.deprecated( - deprecated_in="0.7", details="Use relative_geopotential_thickness_on_hybrid_levels instead." -) -def relative_geopotential_thickness( - alpha: ArrayLike, delta: ArrayLike, t: ArrayLike, q: ArrayLike -) -> ArrayLike: +@deprecation.deprecated(deprecated_in="0.7", details="Use relative_geopotential_thickness_on_hybrid_levels instead.") +def relative_geopotential_thickness(alpha: ArrayLike, delta: ArrayLike, t: ArrayLike, q: ArrayLike) -> ArrayLike: """Calculate the geopotential thickness with respect to the surface on hybrid (IFS model) full-levels. *Deprecated in version 0.7.0* @@ -175,6 +163,10 @@ def relative_geopotential_thickness( array-like Geopotential thickness (m2/s2) of hybrid (IFS model) full-levels with respect to the surface + See Also + -------- + pressure_at_model_levels + Notes ----- ``t`` and ``q`` must contain the same levels in ascending order with respect to @@ -187,10 +179,6 @@ def relative_geopotential_thickness( The computations are described in [IFS-CY47R3-Dynamics]_ Chapter 2, Section 2.2.1. - See also - -------- - pressure_at_model_levels - """ from earthkit.meteo.thermo import specific_gas_constant @@ -251,6 +239,11 @@ def pressure_at_height_levels( number or ndarray pressure at the given height level (Pa) + See Also + -------- + pressure_at_model_levels + relative_geopotential_thickness + Notes ----- ``t`` and ``q`` must contain the same model levels in ascending order with respect to @@ -265,12 +258,6 @@ def pressure_at_height_levels( The pressure at height level is calculated by finding the model level above and below the specified height and interpolating the pressure with linear interpolation. - See also - -------- - pressure_at_model_levels - relative_geopotential_thickness - - """ A = np.asarray(A) B = np.asarray(B) @@ -317,7 +304,8 @@ def pressure_at_height_levels( dphi_below = dphi[below] # print( - # f"tdphi: {tdphi} above: {above} below: {below} dphi_above: {dphi_above} dphi_below {dphi_below} p_full[above]: {p_full[above]} p_full[below]: {p_full[below]}" + # f"tdphi: {tdphi} above: {above} below: {below} dphi_above: {dphi_above} " + # f"dphi_below {dphi_below} p_full[above]: {p_full[above]} p_full[below]: {p_full[below]}" # ) # calculate the interpolation factor @@ -601,15 +589,15 @@ def pressure_on_hybrid_levels( For more details see [IFS-CY47R3-Dynamics]_ Chapter 2, Section 2.2.1. + See Also + -------- + relative_geopotential_thickness_on_hybrid_levels + Examples -------- - :ref:`/examples/hybrid_levels.ipynb` - See also - -------- - relative_geopotential_thickness_on_hybrid_levels - """ if isinstance(output, str): output = (output,) @@ -619,9 +607,7 @@ def pressure_on_hybrid_levels( for out in output: if out not in ["full", "half", "alpha", "delta"]: - raise ValueError( - f"Unknown output type '{out}'. Allowed values are 'full', 'half', 'alpha' or 'delta'." - ) + raise ValueError(f"Unknown output type '{out}'. Allowed values are 'full', 'half', 'alpha' or 'delta'.") if alpha_top not in ["ifs", "arpege"]: raise ValueError(f"Unknown method '{alpha_top}' for pressure calculation. Use 'ifs' or 'arpege'.") @@ -693,9 +679,7 @@ def pressure_on_hybrid_levels( alpha[0, ...] = alpha_top # pressure at highest half-level > 0.1 else: - alpha[0, ...] = ( - 1.0 - p_half_level[0, ...] / (p_half_level[1, ...] - p_half_level[0, ...]) * delta[0, ...] - ) + alpha[0, ...] = 1.0 - p_half_level[0, ...] / (p_half_level[1, ...] - p_half_level[0, ...]) * delta[0, ...] if "full" in output: # calculate pressure on model full-levels @@ -790,7 +774,7 @@ def _compute_relative_geopotential_thickness_on_hybrid_levels( - :ref:`/examples/hybrid_levels.ipynb` - See also + See Also -------- pressure_on_hybrid_levels @@ -863,12 +847,11 @@ def relative_geopotential_thickness_on_hybrid_levels_from_alpha_delta( - :ref:`/examples/hybrid_levels.ipynb` - See also + See Also -------- pressure_on_hybrid_levels """ - xp = array_namespace(alpha, delta, q, t) alpha = xp.asarray(alpha) delta = xp.asarray(delta) @@ -954,7 +937,7 @@ def relative_geopotential_thickness_on_hybrid_levels( -------- - :ref:`/examples/hybrid_levels.ipynb` - See also + See Also -------- pressure_on_hybrid_levels relative_geopotential_thickness_on_hybrid_levels_from_alpha_delta @@ -969,9 +952,7 @@ def relative_geopotential_thickness_on_hybrid_levels( levels = _hybrid_subset(t, A, B, vertical_axis) - alpha, delta = pressure_on_hybrid_levels( - A, B, sp, alpha_top=alpha_top, levels=levels, output=("alpha", "delta") - ) + alpha, delta = pressure_on_hybrid_levels(A, B, sp, alpha_top=alpha_top, levels=levels, output=("alpha", "delta")) # return relative_geopotential_thickness_on_hybrid_levels_from_alpha_delta( # t, q, alpha, delta, vertical_axis=vertical_axis @@ -1055,7 +1036,7 @@ def geopotential_on_hybrid_levels( - :ref:`/examples/hybrid_levels.ipynb` - See also + See Also -------- pressure_on_hybrid_levels relative_geopotential_thickness_on_hybrid_levels @@ -1152,14 +1133,13 @@ def height_on_hybrid_levels( -------- - :ref:`/examples/hybrid_levels.ipynb` - See also + See Also -------- hybrid_level_parameters pressure_on_hybrid_levels geopotential_on_hybrid_levels relative_geopotential_thickness_on_hybrid_levels """ - if h_reference not in ["sea", "ground"]: raise ValueError(f"Unknown '{h_reference=}'. Use 'sea' or 'ground'.") @@ -1299,7 +1279,7 @@ def interpolate_hybrid_to_pressure_levels( -------- - :ref:`/examples/interpolate_hybrid_to_pl.ipynb` - See also + See Also -------- interpolate_monotonic @@ -1455,7 +1435,7 @@ def interpolate_hybrid_to_height_levels( - :ref:`/examples/interpolate_hybrid_to_hl.ipynb` - See also + See Also -------- interpolate_monotonic @@ -1585,7 +1565,7 @@ def interpolate_pressure_to_height_levels( - :ref:`/examples/interpolate_pl_to_hl.ipynb` - See also + See Also -------- interpolate_monotonic diff --git a/src/earthkit/meteo/vertical/interpolation.py b/src/earthkit/meteo/vertical/interpolation.py index 27cbf40f..73118314 100644 --- a/src/earthkit/meteo/vertical/interpolation.py +++ b/src/earthkit/meteo/vertical/interpolation.py @@ -1,7 +1,5 @@ import dataclasses as dc -from typing import Any -from typing import Literal -from typing import Sequence +from typing import Any, Literal, Sequence import numpy as np import xarray as xr @@ -178,8 +176,7 @@ def interpolate_to_pressure_levels( target_values = np.array(sorted(target_p)) * target_factor if np.any((target_values < target_p_min) | (target_values > target_p_max)): raise ValueError( - "target coordinate value out of range " - f"(must be in interval [{target_p_min}, {target_p_max}]Pa)" + f"target coordinate value out of range (must be in interval [{target_p_min}, {target_p_max}]Pa)" ) target = TargetCoordinates( type_of_level="isobaricInPa", @@ -340,9 +337,7 @@ def interpolate_sleve_to_theta_levels( # not monotonous wrt to height tc values are stored in K tc_values = np.array(target_theta) * th_tc_unit_conversions[target_t_units] if np.any((tc_values < th_tc_min) | (tc_values > th_tc_max)): - raise ValueError( - "target coordinate value " f"out of range (must be in interval [{th_tc_min}, {th_tc_max}]K)" - ) + raise ValueError(f"target coordinate value out of range (must be in interval [{th_tc_min}, {th_tc_max}]K)") tc = TargetCoordinates( type_of_level="theta", values=tc_values.tolist(), diff --git a/src/earthkit/meteo/wind/array/__init__.py b/src/earthkit/meteo/wind/array/__init__.py index 1d61d5a7..03e7bf8f 100644 --- a/src/earthkit/meteo/wind/array/__init__.py +++ b/src/earthkit/meteo/wind/array/__init__.py @@ -7,8 +7,6 @@ # nor does it submit to any jurisdiction. # -""" -Wind related functions operating on numpy arrays. -""" +"""Wind related functions operating on numpy arrays.""" from .wind import * # noqa diff --git a/tests/extreme/_data.py b/tests/extreme/_data.py index c5573b21..28fc2955 100644 --- a/tests/extreme/_data.py +++ b/tests/extreme/_data.py @@ -1,493 +1,481 @@ import numpy as np -ens = np.array( - [ - 273.56037903, - 273.61122131, - 273.63012695, - 273.59335327, - 273.57876587, - 273.65365601, - 273.61138916, - 273.56483459, - 273.51663208, - 273.58750916, - 273.65254211, - 273.60098267, - 273.58737183, - 273.60316467, - 273.59606934, - 273.60876465, - 273.50160217, - 273.57923889, - 273.56584167, - 273.49942017, - 273.62069702, - 273.58499146, - 273.58781433, - 273.59477234, - 273.68614197, - 273.60530090, - 273.53808594, - 273.54582214, - 273.65734863, - 273.65567017, - 273.57177734, - 273.55506897, - 273.55067444, - 273.56301880, - 273.54679871, - 273.58064270, - 273.56997681, - 273.61926270, - 273.57879639, - 273.60385132, - 273.63429260, - 273.61221313, - 273.40080261, - 273.52342224, - 273.63037109, - 273.54180908, - 273.56486511, - 273.60005188, - 273.49523926, - 273.56845093, - 273.621521, - ] -) +ens = np.array([ + 273.56037903, + 273.61122131, + 273.63012695, + 273.59335327, + 273.57876587, + 273.65365601, + 273.61138916, + 273.56483459, + 273.51663208, + 273.58750916, + 273.65254211, + 273.60098267, + 273.58737183, + 273.60316467, + 273.59606934, + 273.60876465, + 273.50160217, + 273.57923889, + 273.56584167, + 273.49942017, + 273.62069702, + 273.58499146, + 273.58781433, + 273.59477234, + 273.68614197, + 273.60530090, + 273.53808594, + 273.54582214, + 273.65734863, + 273.65567017, + 273.57177734, + 273.55506897, + 273.55067444, + 273.56301880, + 273.54679871, + 273.58064270, + 273.56997681, + 273.61926270, + 273.57879639, + 273.60385132, + 273.63429260, + 273.61221313, + 273.40080261, + 273.52342224, + 273.63037109, + 273.54180908, + 273.56486511, + 273.60005188, + 273.49523926, + 273.56845093, + 273.621521, +]) ens = ens[:, np.newaxis] -ens_eps = np.array( - [ - 0.00458526611328125, - 0.0051116943359375, - 0.0065460205078125, - 0.00421142578125, - 0.0057830810546875, - 0.008014678955078125, - 0.004108428955078125, - 0.00556182861328125, - 0.004291534423828125, - 0.00505828857421875, - 0.002925872802734375, - 0.007587432861328125, - 0.00222015380859375, - 0.007373809814453125, - 0.0028076171875, - 0.005950927734375, - 0.004131317138671875, - 0.004913330078125, - 0.00797271728515625, - 0.005260467529296875, - 0.0038299560546875, - 0.007537841796875, - 0.00521087646484375, - 0.0087127685546875, - 0.007537841796875, - 0.00247955322265625, - 0.002040863037109375, - 0.0065765380859375, - 0.006870269775390625, - 0.0051422119140625, - 0.004852294921875, - 0.00359344482421875, - 0.00231170654296875, - 0.0048675537109375, - 0.0071563720703125, - 0.0062408447265625, - 0.00475311279296875, - 0.00350189208984375, - 0.00746917724609375, - 0.0047607421875, - 0.00444793701171875, - 0.00389862060546875, - 0.0027923583984375, - 0.004535675048828125, - 0.00621795654296875, - 0.00292205810546875, - 0.00658416748046875, - 0.0046234130859375, - 0.00444793701171875, - 0.0044403076171875, - 0.005218505859375, - ] -) +ens_eps = np.array([ + 0.00458526611328125, + 0.0051116943359375, + 0.0065460205078125, + 0.00421142578125, + 0.0057830810546875, + 0.008014678955078125, + 0.004108428955078125, + 0.00556182861328125, + 0.004291534423828125, + 0.00505828857421875, + 0.002925872802734375, + 0.007587432861328125, + 0.00222015380859375, + 0.007373809814453125, + 0.0028076171875, + 0.005950927734375, + 0.004131317138671875, + 0.004913330078125, + 0.00797271728515625, + 0.005260467529296875, + 0.0038299560546875, + 0.007537841796875, + 0.00521087646484375, + 0.0087127685546875, + 0.007537841796875, + 0.00247955322265625, + 0.002040863037109375, + 0.0065765380859375, + 0.006870269775390625, + 0.0051422119140625, + 0.004852294921875, + 0.00359344482421875, + 0.00231170654296875, + 0.0048675537109375, + 0.0071563720703125, + 0.0062408447265625, + 0.00475311279296875, + 0.00350189208984375, + 0.00746917724609375, + 0.0047607421875, + 0.00444793701171875, + 0.00389862060546875, + 0.0027923583984375, + 0.004535675048828125, + 0.00621795654296875, + 0.00292205810546875, + 0.00658416748046875, + 0.0046234130859375, + 0.00444793701171875, + 0.0044403076171875, + 0.005218505859375, +]) ens_eps = ens_eps[:, np.newaxis] -clim = np.array( - [ - 270.31404114, - 271.04376221, - 271.88911438, - 272.17996216, - 272.29629517, - 272.45881653, - 272.55467224, - 272.67176819, - 272.76577759, - 272.89602661, - 272.93855286, - 273.01394653, - 273.08497620, - 273.13676453, - 273.19113159, - 273.22419739, - 273.25590515, - 273.26646423, - 273.30346680, - 273.31695557, - 273.31858826, - 273.33642578, - 273.37065125, - 273.39843750, - 273.40194702, - 273.44013977, - 273.44854736, - 273.47970581, - 273.49264526, - 273.50416565, - 273.52960205, - 273.54267883, - 273.56050110, - 273.55299377, - 273.56401062, - 273.58006287, - 273.58612061, - 273.62063599, - 273.60931396, - 273.60650635, - 273.61801147, - 273.62141418, - 273.65431213, - 273.63867188, - 273.65670776, - 273.67222595, - 273.65425110, - 273.68618774, - 273.68138123, - 273.67578125, - 273.69137573, - 273.68821716, - 273.69772339, - 273.69319153, - 273.70791626, - 273.73648071, - 273.71696472, - 273.72935486, - 273.74089050, - 273.73986816, - 273.76908875, - 273.75910950, - 273.78601074, - 273.78836060, - 273.80046082, - 273.79504395, - 273.81297302, - 273.80552673, - 273.81306458, - 273.83451843, - 273.83099365, - 273.84126282, - 273.84231567, - 273.85221863, - 273.88308716, - 273.87722778, - 273.89910889, - 273.90881348, - 273.88888550, - 273.90223694, - 273.92579651, - 273.93246460, - 273.93844604, - 273.95280457, - 273.96545410, - 273.98947144, - 273.98995972, - 274.02645874, - 274.05247498, - 274.06059265, - 274.05807495, - 274.10276794, - 274.12309265, - 274.15106201, - 274.16433716, - 274.19328308, - 274.20668030, - 274.27967834, - 274.33123779, - 274.42781067, - 274.55859375, - ] -) +clim = np.array([ + 270.31404114, + 271.04376221, + 271.88911438, + 272.17996216, + 272.29629517, + 272.45881653, + 272.55467224, + 272.67176819, + 272.76577759, + 272.89602661, + 272.93855286, + 273.01394653, + 273.08497620, + 273.13676453, + 273.19113159, + 273.22419739, + 273.25590515, + 273.26646423, + 273.30346680, + 273.31695557, + 273.31858826, + 273.33642578, + 273.37065125, + 273.39843750, + 273.40194702, + 273.44013977, + 273.44854736, + 273.47970581, + 273.49264526, + 273.50416565, + 273.52960205, + 273.54267883, + 273.56050110, + 273.55299377, + 273.56401062, + 273.58006287, + 273.58612061, + 273.62063599, + 273.60931396, + 273.60650635, + 273.61801147, + 273.62141418, + 273.65431213, + 273.63867188, + 273.65670776, + 273.67222595, + 273.65425110, + 273.68618774, + 273.68138123, + 273.67578125, + 273.69137573, + 273.68821716, + 273.69772339, + 273.69319153, + 273.70791626, + 273.73648071, + 273.71696472, + 273.72935486, + 273.74089050, + 273.73986816, + 273.76908875, + 273.75910950, + 273.78601074, + 273.78836060, + 273.80046082, + 273.79504395, + 273.81297302, + 273.80552673, + 273.81306458, + 273.83451843, + 273.83099365, + 273.84126282, + 273.84231567, + 273.85221863, + 273.88308716, + 273.87722778, + 273.89910889, + 273.90881348, + 273.88888550, + 273.90223694, + 273.92579651, + 273.93246460, + 273.93844604, + 273.95280457, + 273.96545410, + 273.98947144, + 273.98995972, + 274.02645874, + 274.05247498, + 274.06059265, + 274.05807495, + 274.10276794, + 274.12309265, + 274.15106201, + 274.16433716, + 274.19328308, + 274.20668030, + 274.27967834, + 274.33123779, + 274.42781067, + 274.55859375, +]) clim = clim[:, np.newaxis] -clim_eps = np.array( - [ - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 6.103515625e-05, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.0001220703125, - 0.00018310546875, - 0.000244140625, - 0.0002441406250, - 0.0002441406250, - 0.0002441406250, - 0.0002441406250, - 0.00030517578125, - 0.0003662109375, - 0.0003662109375, - 0.0003662109375, - 0.0003662109375, - 0.00042724609375, - 0.00048828125, - 0.00048828125, - 0.00048828125, - 0.00054931640625, - 0.0006103515625, - 0.0006103515625, - 0.000732421875, - 0.000732421875, - 0.00079345703125, - 0.0008544921875, - 0.0009765625, - 0.0009765625, - 0.0010986328125, - 0.001220703125, - 0.0013427734375, - 0.00146484375, - 0.00152587890625, - 0.001708984375, - 0.0018310546875, - 0.001953125, - 0.0020751953125, - 0.00225830078125, - 0.00244140625, - 0.0025634765625, - 0.0028076171875, - 0.0030517578125, - 0.0032958984375, - 0.00341796875, - 0.003662109375, - 0.0037841796875, - 0.004150390625, - 0.00433349609375, - 0.0045166015625, - 0.0047607421875, - 0.005126953125, - 0.0054931640625, - 0.005859375, - 0.0062255859375, - 0.006591796875, - 0.007080078125, - 0.00732421875, - 0.0078125, - 0.00830078125, - 0.0089111328125, - 0.009765625, - 0.01177978515625, - 0.01336669921875, - 0.015625, - 0.02294921875, - ] -) +clim_eps = np.array([ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 6.103515625e-05, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.0001220703125, + 0.00018310546875, + 0.000244140625, + 0.0002441406250, + 0.0002441406250, + 0.0002441406250, + 0.0002441406250, + 0.00030517578125, + 0.0003662109375, + 0.0003662109375, + 0.0003662109375, + 0.0003662109375, + 0.00042724609375, + 0.00048828125, + 0.00048828125, + 0.00048828125, + 0.00054931640625, + 0.0006103515625, + 0.0006103515625, + 0.000732421875, + 0.000732421875, + 0.00079345703125, + 0.0008544921875, + 0.0009765625, + 0.0009765625, + 0.0010986328125, + 0.001220703125, + 0.0013427734375, + 0.00146484375, + 0.00152587890625, + 0.001708984375, + 0.0018310546875, + 0.001953125, + 0.0020751953125, + 0.00225830078125, + 0.00244140625, + 0.0025634765625, + 0.0028076171875, + 0.0030517578125, + 0.0032958984375, + 0.00341796875, + 0.003662109375, + 0.0037841796875, + 0.004150390625, + 0.00433349609375, + 0.0045166015625, + 0.0047607421875, + 0.005126953125, + 0.0054931640625, + 0.005859375, + 0.0062255859375, + 0.006591796875, + 0.007080078125, + 0.00732421875, + 0.0078125, + 0.00830078125, + 0.0089111328125, + 0.009765625, + 0.01177978515625, + 0.01336669921875, + 0.015625, + 0.02294921875, +]) clim_eps = clim_eps[:, np.newaxis] -ens_eps2 = np.array( - [ - 0.41961669921875, - 0.4482269287109375, - 0.60272216796875, - 0.6380081176757812, - 0.46539306640625, - 0.461578369140625, - 0.06389617919921875, - 0.6380081176757812, - 0.3070831298828125, - 1.1043548583984375, - 0.9107589721679688, - 0.5931854248046875, - 0.48732757568359375, - 0.20503997802734375, - 0.16117095947265625, - 0.30231475830078125, - 0.5283355712890625, - 0.24127960205078125, - 0.823974609375, - 0.8258819580078125, - 0.7534027099609375, - 0.2307891845703125, - 0.09632110595703125, - 0.08106231689453125, - 0.4138946533203125, - 0.5159378051757812, - 1.7681121826171875, - 1.1615753173828125, - 0.8449554443359375, - 0.5655288696289062, - 0.05054473876953125, - 0.29754638671875, - 0.4482269287109375, - 0.186920166015625, - 0.8544921875, - 0.286102294921875, - 0.0514984130859375, - 1.2903213500976562, - 0.7524490356445312, - 0.820159912109375, - 0.13446807861328125, - 0.6704330444335938, - 0.8459091186523438, - 0.45490264892578125, - 0.4444122314453125, - 0.0896453857421875, - 0.03147125244140625, - 0.4062652587890625, - 0.16117095947265625, - 0.8554458618164062, - 1.5535354614257812, - ] -) +ens_eps2 = np.array([ + 0.41961669921875, + 0.4482269287109375, + 0.60272216796875, + 0.6380081176757812, + 0.46539306640625, + 0.461578369140625, + 0.06389617919921875, + 0.6380081176757812, + 0.3070831298828125, + 1.1043548583984375, + 0.9107589721679688, + 0.5931854248046875, + 0.48732757568359375, + 0.20503997802734375, + 0.16117095947265625, + 0.30231475830078125, + 0.5283355712890625, + 0.24127960205078125, + 0.823974609375, + 0.8258819580078125, + 0.7534027099609375, + 0.2307891845703125, + 0.09632110595703125, + 0.08106231689453125, + 0.4138946533203125, + 0.5159378051757812, + 1.7681121826171875, + 1.1615753173828125, + 0.8449554443359375, + 0.5655288696289062, + 0.05054473876953125, + 0.29754638671875, + 0.4482269287109375, + 0.186920166015625, + 0.8544921875, + 0.286102294921875, + 0.0514984130859375, + 1.2903213500976562, + 0.7524490356445312, + 0.820159912109375, + 0.13446807861328125, + 0.6704330444335938, + 0.8459091186523438, + 0.45490264892578125, + 0.4444122314453125, + 0.0896453857421875, + 0.03147125244140625, + 0.4062652587890625, + 0.16117095947265625, + 0.8554458618164062, + 1.5535354614257812, +]) ens_eps2 = ens_eps2[:, np.newaxis] -clim_eps2 = np.array( - [ - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.48828125, - ] -) +clim_eps2 = np.array([ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.48828125, +]) clim_eps2 = clim_eps2[:, np.newaxis] diff --git a/tests/regimes/test_index_array.py b/tests/regimes/test_index_array.py index 248b98ef..fe1d0abe 100644 --- a/tests/regimes/test_index_array.py +++ b/tests/regimes/test_index_array.py @@ -15,7 +15,6 @@ @pytest.fixture def patterns(): class MockPatterns: - _lat = np.linspace(90.0, 0.0, 91) _lon = np.linspace(60.0, -60.0, 121) _dipole = np.cos(np.deg2rad(_lon[None, :])) * np.cos(np.deg2rad(_lat[:, None]) * 2) @@ -90,18 +89,14 @@ def test_project_generates_weights_by_default(patterns): def test_project_with_single_pattern_return(patterns): - proj = array.project( - np.ones((2, *patterns.shape)), patterns, weights=np.ones(patterns.shape), single=True - ) + proj = array.project(np.ones((2, *patterns.shape)), patterns, weights=np.ones(patterns.shape), single=True) # All patterns are the same; monopole has nonzero projection assert proj["monopole"].shape == (2,) assert np.isclose(proj["monopole"][0], proj["monopole"][1]) def test_project_with_multiple_pattern_return(patterns): - proj = array.project( - np.ones((2, *patterns.shape)), patterns, weights=np.ones(patterns.shape), single=False - ) + proj = array.project(np.ones((2, *patterns.shape)), patterns, weights=np.ones(patterns.shape), single=False) # Second pattern has twice the amplitude; monopole has nonzero projection assert proj["monopole"].shape == (2,) assert np.isclose(proj["monopole"][0], 0.5 * proj["monopole"][1]) diff --git a/tests/regimes/test_index_xarray.py b/tests/regimes/test_index_xarray.py index 544ba8a9..a5d0dbc3 100644 --- a/tests/regimes/test_index_xarray.py +++ b/tests/regimes/test_index_xarray.py @@ -9,8 +9,7 @@ import numpy as np import pytest -from earthkit.meteo.regimes.xarray import project -from earthkit.meteo.regimes.xarray import regime_index +from earthkit.meteo.regimes.xarray import project, regime_index xr = pytest.importorskip("xarray") @@ -33,7 +32,6 @@ def data3d(): @pytest.fixture def patterns(): class MockPatterns: - # Necessary properties to mock Patterns shape = (2, 4) size = 2 * 4 diff --git a/tests/regimes/test_patterns.py b/tests/regimes/test_patterns.py index b42cae38..82ababaa 100644 --- a/tests/regimes/test_patterns.py +++ b/tests/regimes/test_patterns.py @@ -13,7 +13,6 @@ class TestConstantPatterns: - lat = np.linspace(90.0, 0.0, 91) lon = np.linspace(60.0, -60.0, 121) dipole = np.cos(np.deg2rad(lon[None, :])) * np.cos(np.deg2rad(lat[:, None]) * 2) @@ -47,7 +46,6 @@ def test_patterns(self, patterns): class TestModulatedPatterns: - lat = np.linspace(90.0, 0.0, 91) lon = np.linspace(60.0, -60.0, 121) dipole = np.cos(np.deg2rad(lon[None, :])) * np.cos(np.deg2rad(lat[:, None]) * 2) diff --git a/tests/score/array/test_score.py b/tests/score/array/test_score.py index 6c11f525..bfe1b90a 100644 --- a/tests/score/array/test_score.py +++ b/tests/score/array/test_score.py @@ -20,7 +20,7 @@ def crps_quaver2(x, y): """Compute Continuous Ranked Probability Score (CRPS) from Quaver - Used for testing + Used for testing. Parameters ---------- @@ -36,7 +36,6 @@ def crps_quaver2(x, y): The method is described in [Hersbach2000]_. """ - n_ens = x.shape[0] anarr = y earr = x @@ -74,9 +73,7 @@ def crps_quaver2(x, y): def _get_crps_data(): here = os.path.dirname(__file__) sys.path.insert(0, here) - from _crps import ens - from _crps import obs - from _crps import v_ref + from _crps import ens, obs, v_ref return obs, ens, v_ref @@ -153,8 +150,7 @@ def test_crps_quaver2(xp, obs, ens, v_ref): def _get_pearson_data(): here = os.path.dirname(__file__) sys.path.insert(0, here) - from _pearson import SAMPLE_X - from _pearson import SAMPLE_Y + from _pearson import SAMPLE_X, SAMPLE_Y rs = np.array([1.0, -1.0, 0.0, 0.42, -0.13, np.nan]) @@ -165,16 +161,14 @@ def _get_pearson_data(): ymiss = SAMPLE_Y.copy() ymiss[12:20] = np.nan x = np.vstack([SAMPLE_X, SAMPLE_X, SAMPLE_X, SAMPLE_Y, SAMPLE_X, SAMPLE_X]) - y = np.vstack( - [ - SAMPLE_X, - -SAMPLE_X, - SAMPLE_Y, - crs[3] * SAMPLE_X + rs[3] * SAMPLE_Y, - rs[4] * SAMPLE_X + crs[4] * SAMPLE_Y, - ymiss, - ] - ) + y = np.vstack([ + SAMPLE_X, + -SAMPLE_X, + SAMPLE_Y, + crs[3] * SAMPLE_X + rs[3] * SAMPLE_Y, + rs[4] * SAMPLE_X + crs[4] * SAMPLE_Y, + ymiss, + ]) return x.tolist(), y.tolist(), rs.tolist() diff --git a/tests/score/xarray/test_deterministic.py b/tests/score/xarray/test_deterministic.py index 38578e89..078be27a 100644 --- a/tests/score/xarray/test_deterministic.py +++ b/tests/score/xarray/test_deterministic.py @@ -4,17 +4,19 @@ import pytest import xarray as xr -from earthkit.meteo.score import abs_error -from earthkit.meteo.score import cosine_similarity -from earthkit.meteo.score import error -from earthkit.meteo.score import kge -from earthkit.meteo.score import mean_abs_error -from earthkit.meteo.score import mean_error -from earthkit.meteo.score import mean_squared_error -from earthkit.meteo.score import pearson_correlation -from earthkit.meteo.score import root_mean_squared_error -from earthkit.meteo.score import squared_error -from earthkit.meteo.score import standard_deviation_of_error +from earthkit.meteo.score import ( + abs_error, + cosine_similarity, + error, + kge, + mean_abs_error, + mean_error, + mean_squared_error, + pearson_correlation, + root_mean_squared_error, + squared_error, + standard_deviation_of_error, +) from earthkit.meteo.utils.testing import NO_SCORES LATITUDES = [40.0, 41.0, 42.0] @@ -712,27 +714,23 @@ def test_pearson_correlation_with_weights(rng): def test_cosine_similarity(): - fcst = make_dataset( - [ - [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0], [6.0, 7.0, 8.0]], - [[9.0, 10.0, 11.0], [12.0, 13.0, 14.0], [15.0, 16.0, 17.0]], - ] - ) + fcst = make_dataset([ + [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0], [6.0, 7.0, 8.0]], + [[9.0, 10.0, 11.0], [12.0, 13.0, 14.0], [15.0, 16.0, 17.0]], + ]) # Add a bias of 10 and some noise on top - obs = fcst + make_dataset( + obs = fcst + make_dataset([ [ - [ - [10.8784503, 9.95007409, 9.81513764], - [9.31907046, 11.22254134, 9.84547052], - [9.57167218, 9.64786645, 10.53230919], - ], - [ - [10.36544406, 10.41273261, 10.430821], - [12.1416476, 9.59358498, 9.48775727], - [9.18622727, 10.61597942, 11.12897229], - ], - ] - ) + [10.8784503, 9.95007409, 9.81513764], + [9.31907046, 11.22254134, 9.84547052], + [9.57167218, 9.64786645, 10.53230919], + ], + [ + [10.36544406, 10.41273261, 10.430821], + [12.1416476, 9.59358498, 9.48775727], + [9.18622727, 10.61597942, 11.12897229], + ], + ]) result = cosine_similarity(fcst, obs, over=["latitude", "longitude"]) expected_values = np.asarray([0.9208026403441584, 0.9954127943028712]) @@ -749,27 +747,23 @@ def test_cosine_similarity(): @pytest.mark.skipif(NO_SCORES, reason="Scores tests disabled") def test_cosine_similarity_with_weights(): - fcst = make_dataset( - [ - [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0], [6.0, 7.0, 8.0]], - [[9.0, 10.0, 11.0], [12.0, 13.0, 14.0], [15.0, 16.0, 17.0]], - ] - ) + fcst = make_dataset([ + [[0.0, 1.0, 2.0], [3.0, 4.0, 5.0], [6.0, 7.0, 8.0]], + [[9.0, 10.0, 11.0], [12.0, 13.0, 14.0], [15.0, 16.0, 17.0]], + ]) # Add a bias of 10 and some noise on top - obs = fcst + make_dataset( + obs = fcst + make_dataset([ [ - [ - [10.8784503, 9.95007409, 9.81513764], - [9.31907046, 11.22254134, 9.84547052], - [9.57167218, 9.64786645, 10.53230919], - ], - [ - [10.36544406, 10.41273261, 10.430821], - [12.1416476, 9.59358498, 9.48775727], - [9.18622727, 10.61597942, 11.12897229], - ], - ] - ) + [10.8784503, 9.95007409, 9.81513764], + [9.31907046, 11.22254134, 9.84547052], + [9.57167218, 9.64786645, 10.53230919], + ], + [ + [10.36544406, 10.41273261, 10.430821], + [12.1416476, 9.59358498, 9.48775727], + [9.18622727, 10.61597942, 11.12897229], + ], + ]) weights = xr.DataArray( np.array([[1, 1, 1], [1, 2, 1], [1, 1, 1]], dtype=float), dims=["latitude", "longitude"], @@ -911,14 +905,12 @@ def test_kge_modified_return_components(rng): name="r", ) - expected = xr.Dataset( - { - "kge": expected_kge, - "gamma": expected_gamma, - "beta": expected_beta, - "rho": expected_rho, - } - ) + expected = xr.Dataset({ + "kge": expected_kge, + "gamma": expected_gamma, + "beta": expected_beta, + "rho": expected_rho, + }) xr.testing.assert_allclose(result, expected) @@ -989,13 +981,11 @@ def test_kge_original_return_components(rng): name="r", ) - expected = xr.Dataset( - { - "kge": expected_kge, - "alpha": expected_alpha, - "beta": expected_beta, - "rho": expected_rho, - } - ) + expected = xr.Dataset({ + "kge": expected_kge, + "alpha": expected_alpha, + "beta": expected_beta, + "rho": expected_rho, + }) xr.testing.assert_allclose(result, expected) diff --git a/tests/score/xarray/test_ensemble.py b/tests/score/xarray/test_ensemble.py index 969208f3..e7ac3983 100644 --- a/tests/score/xarray/test_ensemble.py +++ b/tests/score/xarray/test_ensemble.py @@ -4,11 +4,7 @@ import pytest import xarray as xr -from earthkit.meteo.score import crps_from_cdf -from earthkit.meteo.score import crps_from_ensemble -from earthkit.meteo.score import crps_from_gaussian -from earthkit.meteo.score import quantile_score -from earthkit.meteo.score import spread +from earthkit.meteo.score import crps_from_cdf, crps_from_ensemble, crps_from_gaussian, quantile_score, spread from earthkit.meteo.utils.testing import NO_SCORES LATITUDES = [40.0, 41.0] @@ -124,20 +120,17 @@ def rng(): def test_spread_without_reference(): """Test spread calculation without a reference, using the ensemble mean.""" - # Forecast data: shape (time, lat, lon, number) - fcst_values = np.array( + fcst_values = np.array([ [ - [ - [[10, 12, 14], [20, 22, 24]], - [[30, 32, 34], [40, 42, 44]], - ], - [ - [[11, 13, 15], [21, 23, 25]], - [[31, 33, 35], [41, 43, 45]], - ], - ] - ) + [[10, 12, 14], [20, 22, 24]], + [[30, 32, 34], [40, 42, 44]], + ], + [ + [[11, 13, 15], [21, 23, 25]], + [[31, 33, 35], [41, 43, 45]], + ], + ]) spread_values = np.std(fcst_values, axis=3, ddof=0) fcst = make_ensemble_dataset(fcst_values) @@ -148,24 +141,20 @@ def test_spread_without_reference(): def test_spread_with_reference(): - fcst_values = np.array( + fcst_values = np.array([ [ - [ - [[10, 12, 14], [20, 22, 24]], - [[30, 32, 34], [40, 42, 44]], - ], - [ - [[11, 13, 15], [21, 23, 25]], - [[31, 33, 35], [41, 43, 45]], - ], - ] - ) - reference_values = np.array( + [[10, 12, 14], [20, 22, 24]], + [[30, 32, 34], [40, 42, 44]], + ], [ - [[12, 24], [36, 48]], - [[7, 19], [31, 43]], - ] - ) + [[11, 13, 15], [21, 23, 25]], + [[31, 33, 35], [41, 43, 45]], + ], + ]) + reference_values = np.array([ + [[12, 24], [36, 48]], + [[7, 19], [31, 43]], + ]) spread_values = np.sqrt( np.mean( (fcst_values - reference_values[:, :, :, np.newaxis]) ** 2, @@ -202,30 +191,26 @@ def test_spread_with_reference(): ) def test_quantile_score(tau: float, expected_scores: np.ndarray): # Forecast data: shape (time, lat, lon, number) - fcst_values = np.array( + fcst_values = np.array([ [ - [ - [np.linspace(0.0, 1.0, 10), np.linspace(10.0, 11.0, 10)], - [np.linspace(20.0, 21.0, 10), np.linspace(30.0, 31.0, 10)], - ], - [ - [np.linspace(1.0, 2.0, 10), np.linspace(11.0, 12.0, 10)], - [np.linspace(21.0, 22.0, 10), np.linspace(31.0, 32.0, 10)], - ], - ] - ) - obs_values = np.array( + [np.linspace(0.0, 1.0, 10), np.linspace(10.0, 11.0, 10)], + [np.linspace(20.0, 21.0, 10), np.linspace(30.0, 31.0, 10)], + ], [ - [ - [-1.0, 10.0], - [21.0, 32.0], - ], - [ - [1.5, 11.1], - [-2.0, 33.0], - ], - ] - ) + [np.linspace(1.0, 2.0, 10), np.linspace(11.0, 12.0, 10)], + [np.linspace(21.0, 22.0, 10), np.linspace(31.0, 32.0, 10)], + ], + ]) + obs_values = np.array([ + [ + [-1.0, 10.0], + [21.0, 32.0], + ], + [ + [1.5, 11.1], + [-2.0, 33.0], + ], + ]) fcst = make_ensemble_dataset(fcst_values) obs = make_deterministic_dataset(obs_values) @@ -256,30 +241,22 @@ def test_quantile_score_invalid_tau(tau: float): def test_crps_from_gaussian(): - fcst_mean_values = np.array( - [ - [[10.0, 20.0], [30.0, 40.0]], - [[15.0, 25.0], [35.0, 45.0]], - ] - ) - fcst_stddev_values = np.array( - [ - [[2.0, 3.0], [4.0, 5.0]], - [[2.5, 3.5], [4.5, 5.5]], - ] - ) - obs_values = np.array( - [ - [[12.0, 19.0], [29.0, 42.0]], - [[14.0, 27.0], [36.0, 44.0]], - ] - ) - expected_crps_values = np.array( - [ - [[1.20488272, 0.83284794], [1.03399925, 1.48344045]], - [[0.74172023, 1.26185368], [1.1399182, 1.35765817]], - ] - ) + fcst_mean_values = np.array([ + [[10.0, 20.0], [30.0, 40.0]], + [[15.0, 25.0], [35.0, 45.0]], + ]) + fcst_stddev_values = np.array([ + [[2.0, 3.0], [4.0, 5.0]], + [[2.5, 3.5], [4.5, 5.5]], + ]) + obs_values = np.array([ + [[12.0, 19.0], [29.0, 42.0]], + [[14.0, 27.0], [36.0, 44.0]], + ]) + expected_crps_values = np.array([ + [[1.20488272, 0.83284794], [1.03399925, 1.48344045]], + [[0.74172023, 1.26185368], [1.1399182, 1.35765817]], + ]) # TODO: Settle on stdev vs stddev fcst = make_gaussian_ensemble_dataset(fcst_mean_values, fcst_stddev_values) @@ -291,12 +268,10 @@ def test_crps_from_gaussian(): def test_crps_from_gaussian_invalid_input(): - fcst_mean_values = np.array( - [ - [[10.0, 20.0], [30.0, 40.0]], - [[15.0, 25.0], [35.0, 45.0]], - ] - ) + fcst_mean_values = np.array([ + [[10.0, 20.0], [30.0, 40.0]], + [[15.0, 25.0], [35.0, 45.0]], + ]) # Missing 'stdev' variable fcst = xr.Dataset( { @@ -308,12 +283,10 @@ def test_crps_from_gaussian_invalid_input(): "longitude": LONGITUDES, }, ) - obs_values = np.array( - [ - [[12.0, 19.0], [29.0, 42.0]], - [[14.0, 27.0], [36.0, 44.0]], - ] - ) + obs_values = np.array([ + [[12.0, 19.0], [29.0, 42.0]], + [[14.0, 27.0], [36.0, 44.0]], + ]) obs = make_deterministic_dataset(obs_values)["2t"] with pytest.raises(ValueError, match="Expected fcst to have 'mean' and 'stdev' data variables"): @@ -674,9 +647,7 @@ def test_crps_from_cdf(): xr.testing.assert_allclose(result["crps"], expected_total_da) xr.testing.assert_allclose(result["underforecast_penalty"], expected_under_da) xr.testing.assert_allclose(result["overforecast_penalty"], expected_over_da) - xr.testing.assert_allclose( - result["crps"], result["underforecast_penalty"] + result["overforecast_penalty"] - ) + xr.testing.assert_allclose(result["crps"], result["underforecast_penalty"] + result["overforecast_penalty"]) xr.testing.assert_allclose(total_only, expected_total_da) diff --git a/tests/solar/test_solar.py b/tests/solar/test_solar.py index 07fd055f..0cec334f 100644 --- a/tests/solar/test_solar.py +++ b/tests/solar/test_solar.py @@ -71,15 +71,11 @@ def test_cos_solar_zenith_angle_1(xp, device, date, lat, lon, v_ref): (datetime.datetime(2024, 4, 22), datetime.datetime(2024, 4, 23), 40.0, 18.0, 4, 0.3110738757), ], ) -def test_cos_solar_zenith_angle_integrated( - xp, device, begin_date, end_date, lat, lon, integration_order, v_ref -): +def test_cos_solar_zenith_angle_integrated(xp, device, begin_date, end_date, lat, lon, integration_order, v_ref): lat = xp.asarray(lat, device=device) lon = xp.asarray(lon, device=device) v_ref = xp.asarray(v_ref, device=device) - v = solar.cos_solar_zenith_angle_integrated( - begin_date, end_date, lat, lon, integration_order=integration_order - ) + v = solar.cos_solar_zenith_angle_integrated(begin_date, end_date, lat, lon, integration_order=integration_order) v_ref = xp.asarray(v_ref, dtype=v.dtype) assert xp.allclose(v, v_ref) diff --git a/tests/stats/_quantile.py b/tests/stats/_quantile.py index b7f4671e..06fb3f37 100644 --- a/tests/stats/_quantile.py +++ b/tests/stats/_quantile.py @@ -5,55 +5,53 @@ # a = np.random.rand(6, 6, 6) # a.ravel()[np.random.choice(a.size, 100, replace=False)] = np.nan -q_test_array = np.array( +q_test_array = np.array([ [ - [ - [np.nan, 0.1052613555, 0.9971272030, 0.0965764379, np.nan, 0.7038418626], - [0.7398191101, 0.9197553237, 0.8342123307, 0.5730007883, 0.3646376776, np.nan], - [0.7042756213, 0.8366483142, np.nan, np.nan, 0.2677947082, np.nan], - [np.nan, 0.8885534671, 0.0449239621, 0.1592224605, np.nan, 0.1392665495], - [np.nan, np.nan, np.nan, 0.5943524428, np.nan, np.nan], - [np.nan, 0.0077392384, np.nan, 0.1523712728, np.nan, np.nan], - ], - [ - [np.nan, 0.1636851736, 0.5704691348, np.nan, 0.6271074513, np.nan], - [0.0999832682, 0.2058537465, np.nan, np.nan, 0.7798492835, np.nan], - [0.8760257570, 0.2803770533, 0.1423455547, 0.9669903821, 0.2018830164, np.nan], - [np.nan, np.nan, 0.0336160817, 0.3247922918, np.nan, 0.6234968194], - [0.8725840576, np.nan, np.nan, np.nan, np.nan, 0.1061212098], - [np.nan, np.nan, 0.9168396876, 0.5625501642, 0.1897713163, np.nan], - ], - [ - [np.nan, np.nan, 0.6229428959, np.nan, 0.4999321007, 0.4452414850], - [0.9418003849, 0.6423931383, np.nan, np.nan, np.nan, np.nan], - [0.0104274852, 0.4989510391, 0.6564322127, np.nan, 0.0280848905, np.nan], - [np.nan, 0.2757613311, 0.0742301095, 0.8303453527, 0.1398419485, 0.0885808277], - [0.2193725299, np.nan, np.nan, 0.1443073125, np.nan, np.nan], - [0.7847564471, 0.6636393090, 0.9790714151, 0.9302066806, np.nan, 0.6253438391], - ], - [ - [np.nan, np.nan, np.nan, 0.8932090343, 0.9545972068, 0.1197048845], - [0.9145126879, np.nan, np.nan, np.nan, np.nan, 0.8536921108], - [0.0551577512, 0.6504191445, np.nan, np.nan, 0.9546922542, 0.5769744779], - [np.nan, 0.0739486398, 0.5475077690, np.nan, 0.6458414986, 0.6194371732], - [np.nan, 0.1611261015, 0.2158372926, 0.3069969111, np.nan, 0.6000506409], - [0.4019886119, np.nan, 0.6872770240, np.nan, 0.2411599117, 0.7803608582], - ], - [ - [np.nan, 0.2142083399, np.nan, 0.6622038575, np.nan, 0.2470772890], - [np.nan, 0.8830092190, 0.8802131079, np.nan, np.nan, np.nan], - [np.nan, 0.9056422867, 0.9835925318, 0.6040755487, 0.9071748415, np.nan], - [np.nan, np.nan, 0.8067820800, 0.5242645299, np.nan, np.nan], - [np.nan, np.nan, np.nan, np.nan, 0.8333722558, np.nan], - [np.nan, 0.6342632073, 0.8304079700, 0.7623822007, np.nan, 0.4217113827], - ], - [ - [0.9059394227, 0.7993060705, np.nan, np.nan, np.nan, 0.3341504699], - [np.nan, np.nan, 0.1373453389, 0.4079903533, np.nan, 0.2313296271], - [np.nan, np.nan, np.nan, 0.0456485354, 0.9577664570, 0.0659499117], - [0.5731775044, 0.4913744057, np.nan, np.nan, 0.4034703651, np.nan], - [0.3089540229, np.nan, 0.3395109486, 0.4711076397, 0.2485256612, np.nan], - [0.6474237082, 0.7014179041, np.nan, 0.3952555796, 0.2498164288, np.nan], - ], - ] -) + [np.nan, 0.1052613555, 0.9971272030, 0.0965764379, np.nan, 0.7038418626], + [0.7398191101, 0.9197553237, 0.8342123307, 0.5730007883, 0.3646376776, np.nan], + [0.7042756213, 0.8366483142, np.nan, np.nan, 0.2677947082, np.nan], + [np.nan, 0.8885534671, 0.0449239621, 0.1592224605, np.nan, 0.1392665495], + [np.nan, np.nan, np.nan, 0.5943524428, np.nan, np.nan], + [np.nan, 0.0077392384, np.nan, 0.1523712728, np.nan, np.nan], + ], + [ + [np.nan, 0.1636851736, 0.5704691348, np.nan, 0.6271074513, np.nan], + [0.0999832682, 0.2058537465, np.nan, np.nan, 0.7798492835, np.nan], + [0.8760257570, 0.2803770533, 0.1423455547, 0.9669903821, 0.2018830164, np.nan], + [np.nan, np.nan, 0.0336160817, 0.3247922918, np.nan, 0.6234968194], + [0.8725840576, np.nan, np.nan, np.nan, np.nan, 0.1061212098], + [np.nan, np.nan, 0.9168396876, 0.5625501642, 0.1897713163, np.nan], + ], + [ + [np.nan, np.nan, 0.6229428959, np.nan, 0.4999321007, 0.4452414850], + [0.9418003849, 0.6423931383, np.nan, np.nan, np.nan, np.nan], + [0.0104274852, 0.4989510391, 0.6564322127, np.nan, 0.0280848905, np.nan], + [np.nan, 0.2757613311, 0.0742301095, 0.8303453527, 0.1398419485, 0.0885808277], + [0.2193725299, np.nan, np.nan, 0.1443073125, np.nan, np.nan], + [0.7847564471, 0.6636393090, 0.9790714151, 0.9302066806, np.nan, 0.6253438391], + ], + [ + [np.nan, np.nan, np.nan, 0.8932090343, 0.9545972068, 0.1197048845], + [0.9145126879, np.nan, np.nan, np.nan, np.nan, 0.8536921108], + [0.0551577512, 0.6504191445, np.nan, np.nan, 0.9546922542, 0.5769744779], + [np.nan, 0.0739486398, 0.5475077690, np.nan, 0.6458414986, 0.6194371732], + [np.nan, 0.1611261015, 0.2158372926, 0.3069969111, np.nan, 0.6000506409], + [0.4019886119, np.nan, 0.6872770240, np.nan, 0.2411599117, 0.7803608582], + ], + [ + [np.nan, 0.2142083399, np.nan, 0.6622038575, np.nan, 0.2470772890], + [np.nan, 0.8830092190, 0.8802131079, np.nan, np.nan, np.nan], + [np.nan, 0.9056422867, 0.9835925318, 0.6040755487, 0.9071748415, np.nan], + [np.nan, np.nan, 0.8067820800, 0.5242645299, np.nan, np.nan], + [np.nan, np.nan, np.nan, np.nan, 0.8333722558, np.nan], + [np.nan, 0.6342632073, 0.8304079700, 0.7623822007, np.nan, 0.4217113827], + ], + [ + [0.9059394227, 0.7993060705, np.nan, np.nan, np.nan, 0.3341504699], + [np.nan, np.nan, 0.1373453389, 0.4079903533, np.nan, 0.2313296271], + [np.nan, np.nan, np.nan, 0.0456485354, 0.9577664570, 0.0659499117], + [0.5731775044, 0.4913744057, np.nan, np.nan, 0.4034703651, np.nan], + [0.3089540229, np.nan, 0.3395109486, 0.4711076397, 0.2485256612, np.nan], + [0.6474237082, 0.7014179041, np.nan, 0.3952555796, 0.2498164288, np.nan], + ], +]) diff --git a/tests/stats/test_stats.py b/tests/stats/test_stats.py index a1cae8cf..cd5009b6 100644 --- a/tests/stats/test_stats.py +++ b/tests/stats/test_stats.py @@ -199,9 +199,7 @@ def test_GumbelDistribution_along_axis(): def test_return_period_identity(): dist = stats.GumbelDistribution.fit([6.0, 5.0, 6.0, 7.0, 9.0, 5.0, 6.0, 7.0]) values = np.linspace(4.0, 10.0, 21) - np.testing.assert_allclose( - stats.return_period_to_value(dist, stats.value_to_return_period(dist, values)), values - ) + np.testing.assert_allclose(stats.return_period_to_value(dist, stats.value_to_return_period(dist, values)), values) def test_return_period_with_frequency(): diff --git a/tests/thermo/test_thermo.py b/tests/thermo/test_thermo.py index 857930b1..f0c895e7 100644 --- a/tests/thermo/test_thermo.py +++ b/tests/thermo/test_thermo.py @@ -7,9 +7,7 @@ # nor does it submit to any jurisdiction. # -""" -Tests for the array level thermo functions -""" +"""Tests for the array level thermo functions.""" import os @@ -38,7 +36,7 @@ def read_data_file(path): def save_test_reference(file_name, data): - """Helper function to save test reference data into csv""" + """Helper function to save test reference data into csv.""" np.savetxt( data_file(file_name), np.column_stack(tuple(data.values())), @@ -114,9 +112,7 @@ def test_vapour_pressure_from_specific_humidity(xp, device, q, p, v_ref): @pytest.mark.parametrize("xp, device", NAMESPACE_DEVICES) -@pytest.mark.parametrize( - "mr, p, v_ref", [([0.0080645161, 0.0183299389], [700, 1000], [895.992614, 2862.662152])] -) +@pytest.mark.parametrize("mr, p, v_ref", [([0.0080645161, 0.0183299389], [700, 1000], [895.992614, 2862.662152])]) def test_vapour_pressure_from_mixing_ratio(xp, device, mr, p, v_ref): mr = xp.asarray(mr, device=device) p = xp.asarray(p, device=device) diff --git a/tests/vertical/test_array_interpolate.py b/tests/vertical/test_array_interpolate.py index 12a4e0ce..db4c8c39 100644 --- a/tests/vertical/test_array_interpolate.py +++ b/tests/vertical/test_array_interpolate.py @@ -63,7 +63,7 @@ def _get_data_2(name): def _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device): - """Test to_pressure with scalar value, scalar pres, scalar target""" + """Test to_pressure with scalar value, scalar pres, scalar target.""" data = xp.asarray(data, device=device) coord = xp.asarray(coord, device=device) target_coord = xp.asarray(target_coord, device=device) @@ -80,7 +80,7 @@ def _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected DATA["pressure_s_s_s"], ) def test_array_interpolate_monotonic_s_s_s_1(data, coord, target_coord, mode, expected_data, xp, device): - """Test with scalar data, scalar coord, scalar target_coord""" + """Test with scalar data, scalar coord, scalar target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -90,7 +90,7 @@ def test_array_interpolate_monotonic_s_s_s_1(data, coord, target_coord, mode, ex DATA["height_s_s_s"], ) def test_array_interpolate_monotonic_s_s_s_2(data, coord, target_coord, mode, expected_data, xp, device): - """Test with scalar data, scalar coord, scalar target_coord""" + """Test with scalar data, scalar coord, scalar target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -100,7 +100,7 @@ def test_array_interpolate_monotonic_s_s_s_2(data, coord, target_coord, mode, ex DATA["pressure_a_a_s"], ) def test_array_interpolate_monotonic_a_a_s(data, coord, target_coord, mode, expected_data, xp, device): - """Test with array data, array coord, scalar target_coord""" + """Test with array data, array coord, scalar target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -110,7 +110,7 @@ def test_array_interpolate_monotonic_a_a_s(data, coord, target_coord, mode, expe DATA["pressure_a_s_s"], ) def test_array_interpolate_monotonic_a_s_s(data, coord, target_coord, mode, expected_data, xp, device): - """Test with array data, scalar coord, scalar target_coord""" + """Test with array data, scalar coord, scalar target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -120,7 +120,7 @@ def test_array_interpolate_monotonic_a_s_s(data, coord, target_coord, mode, expe DATA["pressure_a_s_a"], ) def test_array_interpolate_monotonic_a_s_a(data, coord, target_coord, mode, expected_data, xp, device): - """Test with array data, scalar coord, array target_coord""" + """Test with array data, scalar coord, array target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -130,7 +130,7 @@ def test_array_interpolate_monotonic_a_s_a(data, coord, target_coord, mode, expe DATA["pressure_a_a_a"], ) def test_array_interpolate_monotonic_a_a_a_1(data, coord, target_coord, mode, expected_data, xp, device): - """Test with array data, array coord, array target_coord""" + """Test with array data, array coord, array target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -140,7 +140,7 @@ def test_array_interpolate_monotonic_a_a_a_1(data, coord, target_coord, mode, ex DATA["height_a_a_a"], ) def test_array_interpolate_monotonic_a_a_a_2(data, coord, target_coord, mode, expected_data, xp, device): - """Test with array data, array coord, array target_coord""" + """Test with array data, array coord, array target_coord.""" _check_array_interpolate_monotonic(data, coord, target_coord, mode, expected_data, xp, device) @@ -223,10 +223,8 @@ def test_array_interpolate_monotonic_s_a_a(value, pres, target, mode, xp, device ), ], ) -def test_array_interpolate_monotonic_to_pressure_s_s_s_aux( - data, coord, target_coord, mode, expected_data, xp, device -): - """Test interpolation with auxiliary min/max level data""" +def test_array_interpolate_monotonic_to_pressure_s_s_s_aux(data, coord, target_coord, mode, expected_data, xp, device): + """Test interpolation with auxiliary min/max level data.""" data = xp.asarray(data, device=device) coord = xp.asarray(coord, device=device) target_coord = xp.asarray(target_coord, device=device) @@ -257,10 +255,8 @@ def test_array_interpolate_monotonic_to_pressure_s_s_s_aux( ), ], ) -def test_array_interpolate_monotonic_to_height_s_s_s_aux( - data, coord, target_coord, mode, expected_data, xp, device -): - """Test interpolation with auxiliary min/max level data""" +def test_array_interpolate_monotonic_to_height_s_s_s_aux(data, coord, target_coord, mode, expected_data, xp, device): + """Test interpolation with auxiliary min/max level data.""" data = xp.asarray(data, device=device) coord = xp.asarray(coord, device=device) target_coord = xp.asarray(target_coord, device=device) @@ -314,7 +310,7 @@ def test_array_interpolate_monotonic_to_height_s_s_s_aux( def test_array_interpolate_monotonic_to_height_a_a_a_aux( data, coord, target_coord, aux_data, aux_coord, mode, expected_data, xp, device ): - """Test interpolation with auxiliary min/max level data""" + """Test interpolation with auxiliary min/max level data.""" data = xp.asarray(data, device=device) coord = xp.asarray(coord, device=device) target_coord = xp.asarray(target_coord, device=device) @@ -403,9 +399,7 @@ def test_array_interpolate_hybrid_to_pressure_levels(_kwargs, expected_values, p tolerance = Tolerance({64: (1e-8, 1e-6), 32: (10, 1e-6)}) atol, rtol = tolerance.get(dtype=t.dtype) - assert xp.allclose( - r, r_ref, atol=atol, rtol=rtol, equal_nan=True - ), f"max abs diff={xp.max(xp.abs(r - r_ref))}" + assert xp.allclose(r, r_ref, atol=atol, rtol=rtol, equal_nan=True), f"max abs diff={xp.max(xp.abs(r - r_ref))}" @pytest.mark.parametrize("xp, device", [(_NUMPY_NAMESPACE, "cpu")]) @@ -532,9 +526,7 @@ def test_array_interpolate_hybrid_to_height_levels(_kwargs, expected_values, par tolerance = Tolerance({64: (1e-8, 1e-6), 32: (10, 1e-6)}) atol, rtol = tolerance.get(dtype=t.dtype) - assert xp.allclose( - r, r_ref, atol=atol, rtol=rtol, equal_nan=True - ), f"max abs diff={xp.max(xp.abs(r - r_ref))}" + assert xp.allclose(r, r_ref, atol=atol, rtol=rtol, equal_nan=True), f"max abs diff={xp.max(xp.abs(r - r_ref))}" @pytest.mark.parametrize("xp, device", [(_NUMPY_NAMESPACE, "cpu")]) @@ -656,6 +648,4 @@ def test_array_interpolate_pressure_to_height_levels(_kwargs, expected_values, x tolerance = Tolerance({64: (1e-8, 1e-6), 32: (10, 1e-6)}) atol, rtol = tolerance.get(dtype=t.dtype) - assert xp.allclose( - r, r_ref, atol=atol, rtol=rtol, equal_nan=True - ), f"max abs diff={xp.max(xp.abs(r - r_ref))}" + assert xp.allclose(r, r_ref, atol=atol, rtol=rtol, equal_nan=True), f"max abs diff={xp.max(xp.abs(r - r_ref))}" diff --git a/tests/vertical/test_array_vertical.py b/tests/vertical/test_array_vertical.py index b2770b6d..392c15dc 100644 --- a/tests/vertical/test_array_vertical.py +++ b/tests/vertical/test_array_vertical.py @@ -189,14 +189,12 @@ def test_pressure_on_hybrid_levels_core(index, xp, device): # print("delta diff", repr(xp.max(xp.abs(delta - ref_delta)))) # print("alpha diff", repr(xp.max(xp.abs(alpha - ref_alpha)))) - tolerance = Tolerance( - { - "p_full": {64: (1e-8, 1e-6)}, - "p_half": {64: (1e-8, 1e-6)}, - "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, - "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, - } - ) + tolerance = Tolerance({ + "p_full": {64: (1e-8, 1e-6)}, + "p_half": {64: (1e-8, 1e-6)}, + "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, + "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, + }) atol, rtol = tolerance.get(key="p_full", dtype=sp.dtype) assert xp.allclose(p_full, ref_p_full, atol=atol, rtol=rtol) @@ -257,14 +255,12 @@ def test_pressure_on_hybrid_levels_axis(index, xp, device): delta = xp.moveaxis(delta, vertical_axis, 0) alpha = xp.moveaxis(alpha, vertical_axis, 0) - tolerance = Tolerance( - { - "p_full": {64: (1e-8, 1e-6)}, - "p_half": {64: (1e-8, 1e-6)}, - "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, - "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, - } - ) + tolerance = Tolerance({ + "p_full": {64: (1e-8, 1e-6)}, + "p_half": {64: (1e-8, 1e-6)}, + "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, + "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, + }) atol, rtol = tolerance.get(key="p_full", dtype=sp.dtype) assert xp.allclose(p_full, ref_p_full, atol=atol, rtol=rtol) @@ -281,9 +277,7 @@ def test_pressure_on_hybrid_levels_axis(index, xp, device): @pytest.mark.parametrize("xp, device", NAMESPACE_DEVICES) # @pytest.mark.parametrize("xp, device", [(_NUMPY_NAMESPACE, "cpu")]) @pytest.mark.parametrize("index", [(slice(None), slice(None)), (slice(None), 0), (slice(None), 1)]) -@pytest.mark.parametrize( - "levels", [None, list(range(90, 138)), list(range(137, 90, -1)), [1, 2], [2, 1], [1]] -) +@pytest.mark.parametrize("levels", [None, list(range(90, 138)), list(range(137, 90, -1)), [1, 2], [2, 1], [1]]) @pytest.mark.parametrize( "output", [ @@ -307,7 +301,8 @@ def test_pressure_on_hybrid_levels_output(index, levels, output, xp, device): ref_delta = DATA_HYBRID_CORE.delta ref_alpha = DATA_HYBRID_CORE.alpha - # ref_def = {"full": DATA.p_full, "half": DATA.p_half, "delta": DATA_HYBRID_CORE.delta, "alpha": DATA_HYBRID_CORE.alpha} + # ref_def = {"full": DATA.p_full, "half": DATA.p_half, + # "delta": DATA_HYBRID_CORE.delta, "alpha": DATA_HYBRID_CORE.alpha} # ref = { # key: val # for key, val in ref_def.items() @@ -352,31 +347,29 @@ def test_pressure_on_hybrid_levels_output(index, levels, output, xp, device): res = vertical.pressure_on_hybrid_levels(A, B, sp, levels=levels, alpha_top="ifs", output=output) # atol and rtol for different outputs, due to different precisions in backends - tolerance = Tolerance( - { - "full": {64: (1e-8, 1e-6)}, - "half": {64: (1e-8, 1e-6)}, - "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, - "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, - } - ) + tolerance = Tolerance({ + "full": {64: (1e-8, 1e-6)}, + "half": {64: (1e-8, 1e-6)}, + "delta": {64: (1e-8, 1e-6), 32: (1e-6, 1e-5)}, + "alpha": {64: (1e-8, 1e-6), 32: (1e-4, 1e-5)}, + }) if isinstance(output, str) or len(output) == 1: key = output if isinstance(output, str) else output[0] # print(f"{key=}, max abs diff={xp.max(xp.abs(res - ref[key]))}") atol, rtol = tolerance.get(key=key, dtype=sp.dtype) - assert xp.allclose( - res, ref[key], atol=atol, rtol=rtol - ), f"{key=}, max abs diff={xp.max(xp.abs(res - ref[key]))}" + assert xp.allclose(res, ref[key], atol=atol, rtol=rtol), ( + f"{key=}, max abs diff={xp.max(xp.abs(res - ref[key]))}" + ) else: assert isinstance(res, tuple) assert len(res) == len(output) for key, rd in zip(output, res): atol, rtol = tolerance.get(key=key, dtype=sp.dtype) - assert xp.allclose( - rd, ref[key], atol=atol, rtol=rtol - ), f"{key=}, max abs diff={xp.max(xp.abs(rd - ref[key]))}" + assert xp.allclose(rd, ref[key], atol=atol, rtol=rtol), ( + f"{key=}, max abs diff={xp.max(xp.abs(rd - ref[key]))}" + ) @pytest.mark.parametrize("xp, device", NAMESPACE_DEVICES) diff --git a/tests/vertical/test_xr_monotonic.py b/tests/vertical/test_xr_monotonic.py index 3e9ee800..989abb32 100644 --- a/tests/vertical/test_xr_monotonic.py +++ b/tests/vertical/test_xr_monotonic.py @@ -26,7 +26,7 @@ def _make_xr(xp, data, coord, vdim): - """Make xarray Dataset from data and coord arrays""" + """Make xarray Dataset from data and coord arrays.""" import xarray as xr data_is_scalar = xp.ndim(data[0]) == 0 @@ -110,8 +110,7 @@ def _make_input_xr(data, coord, target_coord, expected_data, vdim, xp=np, device @pytest.mark.skipif(NO_XARRAY, reason="Xarray tests disabled") @pytest.mark.parametrize("ds_input,ds_expected,vdim,mode", make_input("pressure_s_s_s")) def test_xr_interpolate_monotonic_s_s_s(ds_input, ds_expected, vdim, mode): - from earthkit.meteo.vertical.interpolation import TargetCoordinates - from earthkit.meteo.vertical.interpolation import interpolate_monotonic + from earthkit.meteo.vertical.interpolation import TargetCoordinates, interpolate_monotonic tc = TargetCoordinates("isobaricInhPa", ds_expected[vdim]) observed = interpolate_monotonic(ds_input.data, ds_input[vdim], tc, mode, vertical_dim=vdim) @@ -121,8 +120,7 @@ def test_xr_interpolate_monotonic_s_s_s(ds_input, ds_expected, vdim, mode): @pytest.mark.skipif(NO_XARRAY, reason="Xarray tests disabled") @pytest.mark.parametrize("ds_input,ds_expected,vdim,mode", make_input("pressure_a_a_s")) def test_xr_interpolate_monotonic_a_a_s(ds_input, ds_expected, vdim, mode): - from earthkit.meteo.vertical.interpolation import TargetCoordinates - from earthkit.meteo.vertical.interpolation import interpolate_monotonic + from earthkit.meteo.vertical.interpolation import TargetCoordinates, interpolate_monotonic tc = TargetCoordinates("isobaricInhPa", ds_expected[vdim].values) observed = interpolate_monotonic(ds_input.data, ds_input.coord, tc, mode, vertical_dim=vdim) @@ -132,8 +130,7 @@ def test_xr_interpolate_monotonic_a_a_s(ds_input, ds_expected, vdim, mode): @pytest.mark.skipif(NO_XARRAY, reason="Xarray tests disabled") @pytest.mark.parametrize("ds_input,ds_expected,vdim,mode", make_input("pressure_a_s_s")) def test_xr_interpolate_monotonic_a_s_s(ds_input, ds_expected, vdim, mode): - from earthkit.meteo.vertical.interpolation import TargetCoordinates - from earthkit.meteo.vertical.interpolation import interpolate_monotonic + from earthkit.meteo.vertical.interpolation import TargetCoordinates, interpolate_monotonic tc = TargetCoordinates("isobaricInhPa", ds_expected[vdim].values) observed = interpolate_monotonic(ds_input.data, ds_input[vdim], tc, mode, vertical_dim=vdim) diff --git a/tests/vertical/test_xr_theta.py b/tests/vertical/test_xr_theta.py index 340f6b64..f82e9809 100644 --- a/tests/vertical/test_xr_theta.py +++ b/tests/vertical/test_xr_theta.py @@ -49,9 +49,7 @@ def test_interpolate_to_pressure(interpolation, input_ds, output_ds): target_p = [40.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1100.0] target_p_units = "hPa" - observed = interpolate_to_pressure_levels( - t, p, target_p, target_p_units, interpolation=interpolation - ).values + observed = interpolate_to_pressure_levels(t, p, target_p, target_p_units, interpolation=interpolation).values expected = output_ds("pressure", interpolation)["T"].values.squeeze() np.testing.assert_allclose(observed, expected, 1e-5, 1e-7) diff --git a/tests/wind/test_wind.py b/tests/wind/test_wind.py index f5656448..ee98b859 100644 --- a/tests/wind/test_wind.py +++ b/tests/wind/test_wind.py @@ -251,21 +251,19 @@ def test_coriolis(lat, v_ref, xp, device): 6, [0, 1, 2, 3, 4], True, - np.array( + np.array([ + [0.0000000000, 0.0000000000, 1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000], + [1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000], + [0.0000000000, 0.0000000000, 1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000], [ - [0.0000000000, 0.0000000000, 1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000], - [1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000], - [0.0000000000, 0.0000000000, 1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000], - [ - 2.0000000000, - 0.0000000000, - 0.0000000000, - 0.0000000000, - 0.0000000000, - 1.0000000000, - ], - ] - ) + 2.0000000000, + 0.0000000000, + 0.0000000000, + 0.0000000000, + 0.0000000000, + 1.0000000000, + ], + ]) * 100 / 11.0, [ @@ -299,9 +297,7 @@ def test_coriolis(lat, v_ref, xp, device): ], ) def test_windrose_1(sp, d, sectors, sp_bins, percent, v_ref, dir_bin_ref, xp, device): - sp, d, sp_bins, v_ref, dir_bin_ref = ( - xp.asarray(x, device=device) for x in [sp, d, sp_bins, v_ref, dir_bin_ref] - ) + sp, d, sp_bins, v_ref, dir_bin_ref = (xp.asarray(x, device=device) for x in [sp, d, sp_bins, v_ref, dir_bin_ref]) dir_bin_ref = xp.astype(dir_bin_ref, sp.dtype) @@ -315,9 +311,7 @@ def test_windrose_1(sp, d, sectors, sp_bins, percent, v_ref, dir_bin_ref, xp, de @pytest.mark.parametrize("xp, device", NAMESPACE_DEVICES) -@pytest.mark.parametrize( - "sp,d,sectors,sp_bins", [(3.4, 90.01, 0, [0, 1]), (3.4, 90.01, 6, [0]), (3.4, 90.01, 6, None)] -) +@pytest.mark.parametrize("sp,d,sectors,sp_bins", [(3.4, 90.01, 0, [0, 1]), (3.4, 90.01, 6, [0]), (3.4, 90.01, 6, None)]) def test_windrose_invalid(sp, d, sectors, sp_bins, xp, device): if sp_bins is not None: sp_bins = xp.asarray(sp_bins)