From e3414f44128ae476d8fbc7c1a23006b2575d5d2e Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Sun, 24 Mar 2024 16:58:14 +0000 Subject: [PATCH 1/4] Add generate_latlon method --- earthkit/geo/__init__.py | 2 + earthkit/geo/latlon.py | 63 ++++++++++++++++++++++++++++++ tests/grid/test_generate_latlon.py | 41 +++++++++++++++++++ 3 files changed, 106 insertions(+) create mode 100644 earthkit/geo/latlon.py create mode 100644 tests/grid/test_generate_latlon.py diff --git a/earthkit/geo/__init__.py b/earthkit/geo/__init__.py index 2641196..a4fff31 100644 --- a/earthkit/geo/__init__.py +++ b/earthkit/geo/__init__.py @@ -22,9 +22,11 @@ nearest_point_haversine, nearest_point_kdtree, ) +from earthkit.geo.latlon import generate_latlon __all__ = [ "GeoKDTree", + "generate_latlon", "haversine_distance", "nearest_point_haversine", "nearest_point_kdtree", diff --git a/earthkit/geo/latlon.py b/earthkit/geo/latlon.py new file mode 100644 index 0000000..c2d6ece --- /dev/null +++ b/earthkit/geo/latlon.py @@ -0,0 +1,63 @@ +# (C) Copyright 2024 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +from . import constants + + +def generate_latlon(dx, dy): + """Generate latitude and longitude arrays for a global regular latitude-longitude grid. + + Parameters + ---------- + dx: number + Increment in longitudes (degrees) + dy: number + Increment in latitudes (degrees) + + Returns + ------- + ndarray + Latitudes (degrees) with the shape of the grid + ndarray + Longitudes (degrees) with the shape of the grid + + + The global regular latitude-longitude grid is formed as follows: + + - the points start at the north pole and go down to the south pole from west to east in each latitude band + - the shape of the grid is (ny, nx) where ny is the number of latitudes and nx is the number of longitudes + + Examples + -------- + >>> import earthkit.geo + >>> lat, lon = earthkit.geo.generate_latlon(120, 45) + >>> lat + array([[ 90., 90., 90.], + [ 45., 45., 45.], + [ 0., 0., 0.], + [-45., -45., -45.], + [-90., -90., -90.]]) + >>> lon + array([[ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.]]) + + """ + import numpy as np + + lat_v = np.linspace( + constants.north, + constants.south, + int((constants.north - constants.south) / dy) + 1, + ) + lon_v = np.linspace(0, constants.full_circle - dx, int(constants.full_circle / dx)) + lon, lat = np.meshgrid(lon_v, lat_v) + return lat, lon diff --git a/tests/grid/test_generate_latlon.py b/tests/grid/test_generate_latlon.py new file mode 100644 index 0000000..be1553e --- /dev/null +++ b/tests/grid/test_generate_latlon.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 + +# (C) Copyright 2024 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + +import numpy as np +import pytest + +from earthkit.geo import generate_latlon + + +@pytest.mark.parametrize( + "dx,dy,shape,vals_x,vals_y", + [ + (20, 20, (10, 18), np.arange(0, 360, 20), np.arange(90, -90 - 20, -20)), + (30, 20, (10, 12), np.arange(0, 360, 30), np.arange(90, -90 - 20, -20)), + ], +) +def test_generate_latlon(dx, dy, shape, vals_x, vals_y): + lat, lon = generate_latlon(dx, dy) + assert lat.shape == shape + assert lon.shape == shape + + ny = lat.shape[0] + nx = lat.shape[1] + for i in range(ny): + first = lat[i, 0] + assert np.allclose(lat[i, :], np.ones(nx) * first) + + for i in range(nx): + first = lon[0, i] + assert np.allclose(lon[:, i], np.ones(ny) * first) + + assert np.allclose(lat[:, 0], vals_y) + assert np.allclose(lon[0, :], vals_x) From aa8f4d105b565b0f76e8a52e0aa4dacafd53ff76 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Thu, 28 Mar 2024 09:57:26 +0000 Subject: [PATCH 2/4] Add to_latlon method --- earthkit/geo/__init__.py | 2 - earthkit/geo/constants.py | 6 + earthkit/geo/grid.py | 119 +++++++++ earthkit/geo/latlon.py | 63 ----- tests/data/global_ll_ref.json | 402 +++++++++++++++++++++++++++++ tests/grid/test_generate_latlon.py | 71 ++++- 6 files changed, 593 insertions(+), 70 deletions(-) create mode 100644 earthkit/geo/grid.py delete mode 100644 earthkit/geo/latlon.py create mode 100644 tests/data/global_ll_ref.json diff --git a/earthkit/geo/__init__.py b/earthkit/geo/__init__.py index a4fff31..2641196 100644 --- a/earthkit/geo/__init__.py +++ b/earthkit/geo/__init__.py @@ -22,11 +22,9 @@ nearest_point_haversine, nearest_point_kdtree, ) -from earthkit.geo.latlon import generate_latlon __all__ = [ "GeoKDTree", - "generate_latlon", "haversine_distance", "nearest_point_haversine", "nearest_point_kdtree", diff --git a/earthkit/geo/constants.py b/earthkit/geo/constants.py index e0f90a9..ae18e6e 100644 --- a/earthkit/geo/constants.py +++ b/earthkit/geo/constants.py @@ -18,6 +18,12 @@ full_circle = 360 r"""Full circle in degrees""" +half_circle = 180 +r"""Half circle in degrees""" + +quarter_circle = 90 +r"""Quarter circle in degrees""" + north = 90 r"""Latitude of the north pole in degrees""" diff --git a/earthkit/geo/grid.py b/earthkit/geo/grid.py new file mode 100644 index 0000000..795b226 --- /dev/null +++ b/earthkit/geo/grid.py @@ -0,0 +1,119 @@ +# (C) Copyright 2024 ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. +# + + +from . import constants + + +def to_latlon(grid): + """Generate latitude and longitude arrays for a given ``grid``. + + Parameters + ---------- + grid: dict + Grid specification. At the moment, only global regular latitude-longitude + grids are supported. The required format is as follows:: + + {"grid": [dlon, dlat]} + + Returns + ------- + ndarray + Latitudes (degrees) with the shape of the grid + ndarray + Longitudes (degrees) with the shape of the grid + + + The global regular latitude-longitude grid is defined as follows: + + - the points start at the north pole and go down to the south pole from + west (0°) to east in each latitude band + - the **shape** of the grid is (ny, nx) where ny is the number of latitudes + and nx is the number of longitudes + + Examples + -------- + >>> import earthkit.geo + >>> lat, lon = earthkit.geo.grid.to_latlon({"grid": [120, 45]}) + >>> lat + array([[ 90., 90., 90.], + [ 45., 45., 45.], + [ 0., 0., 0.], + [-45., -45., -45.], + [-90., -90., -90.]]) + >>> lon + array([[ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.], + [ 0., 120., 240.]]) + + """ + if not isinstance(grid, dict): + raise TypeError(f"Invalid grid type {type(grid)}, must be a dict.") + + grid = dict(grid) + dlon, dlat = grid.pop("grid") + if grid: + raise KeyError(f"Unexpected keys in grid = {list(grid.keys())}") + + if dlat <= 0: + raise ValueError(f"Increment in latitude must be positive: {dlat=}") + if dlon <= 0: + raise ValueError(f"Increment in longitude must be positive: {dlon=}") + + import numpy as np + + eps_lat = dlat / 1e2 + eps_lon = dlon / 1e2 + + north = constants.north + south = constants.south + west = 0.0 + east = constants.full_circle + + # north-south + ny2 = round(constants.quarter_circle / dlat) + half_lat_range = ny2 * dlat + delta = constants.quarter_circle - half_lat_range + if delta < -eps_lat: + ny2 -= 1 + half_lat_range = ny2 * dlat + delta = constants.quarter_circle - half_lat_range + + ny = ny2 * 2 + 1 + + if abs(delta) < eps_lat: + north = constants.north + south = constants.south + else: + north = half_lat_range + south = -north + + assert ny % 2 == 1 + lats = np.array([north - i * dlat for i in range(ny)]) + lats[-1] = south + + # west-east + west = 0.0 + east = constants.full_circle + + nx = round(constants.full_circle / dlon) + lon_range = nx * dlon + delta = constants.full_circle - lon_range + if abs(delta) < eps_lon or delta < 0: + nx = nx - 1 + nx = nx + 1 + east = (nx - 1) * dlon + + lons = np.array([west + i * dlon for i in range(nx)]) + lons[-1] = east + + lon, lat = np.meshgrid(lons, lats) + return lat, lon diff --git a/earthkit/geo/latlon.py b/earthkit/geo/latlon.py deleted file mode 100644 index c2d6ece..0000000 --- a/earthkit/geo/latlon.py +++ /dev/null @@ -1,63 +0,0 @@ -# (C) Copyright 2024 ECMWF. -# -# This software is licensed under the terms of the Apache Licence Version 2.0 -# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -# In applying this licence, ECMWF does not waive the privileges and immunities -# granted to it by virtue of its status as an intergovernmental organisation -# nor does it submit to any jurisdiction. -# - -from . import constants - - -def generate_latlon(dx, dy): - """Generate latitude and longitude arrays for a global regular latitude-longitude grid. - - Parameters - ---------- - dx: number - Increment in longitudes (degrees) - dy: number - Increment in latitudes (degrees) - - Returns - ------- - ndarray - Latitudes (degrees) with the shape of the grid - ndarray - Longitudes (degrees) with the shape of the grid - - - The global regular latitude-longitude grid is formed as follows: - - - the points start at the north pole and go down to the south pole from west to east in each latitude band - - the shape of the grid is (ny, nx) where ny is the number of latitudes and nx is the number of longitudes - - Examples - -------- - >>> import earthkit.geo - >>> lat, lon = earthkit.geo.generate_latlon(120, 45) - >>> lat - array([[ 90., 90., 90.], - [ 45., 45., 45.], - [ 0., 0., 0.], - [-45., -45., -45.], - [-90., -90., -90.]]) - >>> lon - array([[ 0., 120., 240.], - [ 0., 120., 240.], - [ 0., 120., 240.], - [ 0., 120., 240.], - [ 0., 120., 240.]]) - - """ - import numpy as np - - lat_v = np.linspace( - constants.north, - constants.south, - int((constants.north - constants.south) / dy) + 1, - ) - lon_v = np.linspace(0, constants.full_circle - dx, int(constants.full_circle / dx)) - lon, lat = np.meshgrid(lon_v, lat_v) - return lat, lon diff --git a/tests/data/global_ll_ref.json b/tests/data/global_ll_ref.json new file mode 100644 index 0000000..572bd17 --- /dev/null +++ b/tests/data/global_ll_ref.json @@ -0,0 +1,402 @@ +{ + "0.1": { + "grid": [ + 0.1, + 0.1 + ], + "shape": [ + 1801, + 3600 + ], + "area": [ + 90, + 0, + -90, + 359.9 + ] + }, + "0.125": { + "grid": [ + 0.125, + 0.125 + ], + "shape": [ + 1441, + 2880 + ], + "area": [ + 90, + 0, + -90, + 359.875 + ] + }, + "0.15": { + "grid": [ + 0.15, + 0.15 + ], + "shape": [ + 1201, + 2400 + ], + "area": [ + 90, + 0, + -90, + 359.85 + ] + }, + "0.2": { + "grid": [ + 0.2, + 0.2 + ], + "shape": [ + 901, + 1800 + ], + "area": [ + 90, + 0, + -90, + 359.8 + ] + }, + "0.25": { + "grid": [ + 0.25, + 0.25 + ], + "shape": [ + 721, + 1440 + ], + "area": [ + 90, + 0, + -90, + 359.75 + ] + }, + "0.3": { + "grid": [ + 0.3, + 0.3 + ], + "shape": [ + 601, + 1200 + ], + "area": [ + 90, + 0, + -90, + 359.7 + ] + }, + "0.4": { + "grid": [ + 0.4, + 0.4 + ], + "shape": [ + 451, + 900 + ], + "area": [ + 90, + 0, + -90, + 359.6 + ] + }, + "0.5": { + "grid": [ + 0.5, + 0.5 + ], + "shape": [ + 361, + 720 + ], + "area": [ + 90, + 0, + -90, + 359.5 + ] + }, + "0.6": { + "grid": [ + 0.6, + 0.6 + ], + "shape": [ + 301, + 600 + ], + "area": [ + 90, + 0, + -90, + 359.4 + ] + }, + "0.7": { + "grid": [ + 0.7, + 0.7 + ], + "shape": [ + 257, + 515 + ], + "area": [ + 89.6, + 0, + -89.6, + 359.8 + ] + }, + "0.75": { + "grid": [ + 0.75, + 0.75 + ], + "shape": [ + 241, + 480 + ], + "area": [ + 90, + 0, + -90, + 359.25 + ] + }, + "0.8": { + "grid": [ + 0.8, + 0.8 + ], + "shape": [ + 225, + 450 + ], + "area": [ + 89.6, + 0, + -89.6, + 359.2 + ] + }, + "0.9": { + "grid": [ + 0.9, + 0.9 + ], + "shape": [ + 201, + 400 + ], + "area": [ + 90, + 0, + -90, + 359.1 + ] + }, + "1": { + "grid": [ + 1, + 1 + ], + "shape": [ + 181, + 360 + ], + "area": [ + 90, + 0, + -90, + 359 + ] + }, + "1.2": { + "grid": [ + 1.2, + 1.2 + ], + "shape": [ + 151, + 300 + ], + "area": [ + 90, + 0, + -90, + 358.8 + ] + }, + "1.25": { + "grid": [ + 1.25, + 1.25 + ], + "shape": [ + 145, + 288 + ], + "area": [ + 90, + 0, + -90, + 358.75 + ] + }, + "1.4": { + "grid": [ + 1.4, + 1.4 + ], + "shape": [ + 129, + 258 + ], + "area": [ + 89.6, + 0, + -89.6, + 359.8 + ] + }, + "1.5": { + "grid": [ + 1.5, + 1.5 + ], + "shape": [ + 121, + 240 + ], + "area": [ + 90, + 0, + -90, + 358.5 + ] + }, + "1.6": { + "grid": [ + 1.6, + 1.6 + ], + "shape": [ + 113, + 225 + ], + "area": [ + 89.6, + 0, + -89.6, + 358.4 + ] + }, + "1.8": { + "grid": [ + 1.8, + 1.8 + ], + "shape": [ + 101, + 200 + ], + "area": [ + 90, + 0, + -90, + 358.2 + ] + }, + "2": { + "grid": [ + 2, + 2 + ], + "shape": [ + 91, + 180 + ], + "area": [ + 90, + 0, + -90, + 358 + ] + }, + "2.5": { + "grid": [ + 2.5, + 2.5 + ], + "shape": [ + 73, + 144 + ], + "area": [ + 90, + 0, + -90, + 357.5 + ] + }, + "5": { + "grid": [ + 5, + 5 + ], + "shape": [ + 37, + 72 + ], + "area": [ + 90, + 0, + -90, + 355 + ] + }, + "10": { + "grid": [ + 10, + 10 + ], + "shape": [ + 19, + 36 + ], + "area": [ + 90, + 0, + -90, + 350 + ] + }, + "20": { + "grid": [ + 20, + 20 + ], + "shape": [ + 9, + 19 + ], + "area": [ + 90, + 0, + -90, + 340 + ] + } +} diff --git a/tests/grid/test_generate_latlon.py b/tests/grid/test_generate_latlon.py index be1553e..71ffc1c 100644 --- a/tests/grid/test_generate_latlon.py +++ b/tests/grid/test_generate_latlon.py @@ -9,21 +9,68 @@ # nor does it submit to any jurisdiction. # +import os + import numpy as np import pytest -from earthkit.geo import generate_latlon +from earthkit.geo.grid import to_latlon + +DATA_PATH = os.path.join(os.path.dirname(__file__), "data") + + +def file_in_testdir(filename): + return os.path.join(DATA_PATH, filename) + + +def _global_grids(): + import json + + from earthkit.data.utils.paths import earthkit_conf_file + + with open(earthkit_conf_file("global_grids.json"), "r") as f: + d = json.load(f) + + for _, v in d.items(): + yield v + + +# these are all the cases for the earthkit-regrid target latlon grids +@pytest.mark.parametrize("grid", _global_grids()) +def test_to_latlon_1(grid): + dx, dy = grid["grid"] + lat, lon = to_latlon({"grid": [dx, dy]}) + + assert list(lat.shape) == grid["shape"] + assert list(lon.shape) == grid["shape"] + + # bounding box + north = lat[0, 0] + south = lat[-1, 0] + west = lon[0, 0] + east = lon[0, -1] + ref_north, ref_west, ref_south, ref_east = grid["area"] + assert np.isclose(north, ref_north) + assert np.isclose(south, ref_south) + assert np.isclose(west, ref_west) + assert np.isclose(east, ref_east) + + # equator + ny = lat.shape[0] + assert ny % 2 == 1 + assert np.isclose(lat[int((ny - 1) / 2), 0], 0.0) @pytest.mark.parametrize( "dx,dy,shape,vals_x,vals_y", [ - (20, 20, (10, 18), np.arange(0, 360, 20), np.arange(90, -90 - 20, -20)), - (30, 20, (10, 12), np.arange(0, 360, 30), np.arange(90, -90 - 20, -20)), + (10, 10, (19, 36), np.arange(0, 360, 10), np.arange(90, -90 - 10, -10)), + (30, 20, (9, 12), np.arange(0, 360, 30), np.arange(80, -80 - 20, -20)), ], ) -def test_generate_latlon(dx, dy, shape, vals_x, vals_y): - lat, lon = generate_latlon(dx, dy) +def test_to_latlon_2(dx, dy, shape, vals_x, vals_y): + gs = {"grid": [dx, dy]} + lat, lon = to_latlon(gs) assert lat.shape == shape assert lon.shape == shape @@ -39,3 +86,17 @@ def test_generate_latlon(dx, dy, shape, vals_x, vals_y): assert np.allclose(lat[:, 0], vals_y) assert np.allclose(lon[0, :], vals_x) + + +@pytest.mark.parametrize( + "grid,error", + [ + (1, TypeError), + ({"grid": 1}, TypeError), + ({"grid": [-1, 1]}, ValueError), + ({"grid": [1, 1], "area": [70, 0, 40, 25]}, KeyError), + ], +) +def test_to_latlon_invalid(grid, error): + with pytest.raises(error): + to_latlon(grid) From 8a7e142612568a0b785260a9e7d73373c0f40334 Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Thu, 28 Mar 2024 10:03:58 +0000 Subject: [PATCH 3/4] Add to_latlon method --- tests/data/global_ll_ref.json | 6 +++--- tests/grid/test_generate_latlon.py | 6 ++---- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/tests/data/global_ll_ref.json b/tests/data/global_ll_ref.json index 572bd17..1952c71 100644 --- a/tests/data/global_ll_ref.json +++ b/tests/data/global_ll_ref.json @@ -390,12 +390,12 @@ ], "shape": [ 9, - 19 + 18 ], "area": [ - 90, + 80, 0, - -90, + -80, 340 ] } diff --git a/tests/grid/test_generate_latlon.py b/tests/grid/test_generate_latlon.py index 71ffc1c..1fad8d4 100644 --- a/tests/grid/test_generate_latlon.py +++ b/tests/grid/test_generate_latlon.py @@ -16,7 +16,7 @@ from earthkit.geo.grid import to_latlon -DATA_PATH = os.path.join(os.path.dirname(__file__), "data") +DATA_PATH = os.path.join(os.path.dirname(os.path.dirname(__file__)), "data") def file_in_testdir(filename): @@ -26,9 +26,7 @@ def file_in_testdir(filename): def _global_grids(): import json - from earthkit.data.utils.paths import earthkit_conf_file - - with open(earthkit_conf_file("global_grids.json"), "r") as f: + with open(file_in_testdir("global_ll_ref.json"), "r") as f: d = json.load(f) for _, v in d.items(): From 7587a3046f9c41c5e7d27791e36b1dcb07410b1f Mon Sep 17 00:00:00 2001 From: Sandor Kertesz Date: Wed, 8 May 2024 17:04:54 +0100 Subject: [PATCH 4/4] Merge from develop --- src/earthkit/geo/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/earthkit/geo/constants.py b/src/earthkit/geo/constants.py index 822cdb8..10a8085 100644 --- a/src/earthkit/geo/constants.py +++ b/src/earthkit/geo/constants.py @@ -21,4 +21,4 @@ r"""Full angle in degrees""" RIGHT_ANGLE = 90.0 -r"""Right angle in degrees""" \ No newline at end of file +r"""Right angle in degrees"""