From 8f9bdaa704b33ec8d21585ba7e0d81fe57c9125b Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 23 Aug 2023 15:27:38 +0100 Subject: [PATCH 01/10] allow extract point to support other grids #2177 --- esmvalcore/preprocessor/_regrid.py | 49 +++++++++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 1008cf67ad..5ec3a9c47f 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -18,6 +18,7 @@ from geopy.geocoders import Nominatim from iris.analysis import AreaWeighted, Linear, Nearest, UnstructuredNearest from iris.util import broadcast_to_shape +import cartopy.crs as ccrs from ..cmor._fixes.shared import add_altitude_from_plev, add_plev_from_altitude from ..cmor.table import CMOR_TABLES @@ -460,9 +461,55 @@ def extract_point(cube, latitude, longitude, scheme): raise ValueError(msg) point = [('latitude', latitude), ('longitude', longitude)] - cube = cube.interpolate(point, scheme=scheme) + cube = get_pt(cube, point, scheme) return cube +def get_pt(cube, point, scheme): + """ + Extract data at a single point from cubes with any coordinate + system or none (if lat-lon only) + + Parameters + ---------- + cube : cube + The source cube to extract a point from. + + latitude, longitude : float, or array of float + The latitude and longitude of the point. + + scheme : str + The interpolation scheme. 'linear' or 'nearest'. No default. + + Returns + ------- + :py:class:`~iris.cube.Cube` + Returns a cube with the extracted point(s) + + Raises + ------ + ValueError + if cube doesn't have a coordinate system and isn't on a lat-lon grid + """ + x_coord = cube.coord(axis='X', dim_coords=True) + y_coord = cube.coord(axis='Y', dim_coords=True) + + # some lon-lat cubes don't have a coordinate system + # check if it is lon-lat and if so just use interpolate + # as it is + if cube.coord_system() is None: + # if coords are equal to lat and lon then + # assume we have a lat-lon grid + + if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': + return cube.interpolate(point, scheme=scheme) + else: + raise ValueError('Can only interpolate lat-lon grids if no coordinate system on cube') + else: + # convert the target point(s) to lat lon and do the interpolation + ll = ccrs.Geodetic() + trpoints = cube.coord_system().as_cartopy_crs().transform_point(point[1][1], point[0][1], ll) + trpoints = [(y_coord.name(), trpoints[1]),(x_coord.name(), trpoints[0])] + return cube.interpolate(trpoints, scheme=scheme) def is_dataset(dataset): """Test if something is an `esmvalcore.dataset.Dataset`.""" From cca63ba94ca6e2c9c2cdbd14cb542912ab146b15 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 30 Aug 2023 11:14:10 +0100 Subject: [PATCH 02/10] extra fixes to work with lists of points - shapes still wrong. #2177 --- esmvalcore/preprocessor/_regrid.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 5ec3a9c47f..85681fbddc 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -9,6 +9,7 @@ from copy import deepcopy from decimal import Decimal from pathlib import Path +from turtle import pd from typing import Dict import dask.array as da @@ -503,12 +504,24 @@ def get_pt(cube, point, scheme): if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': return cube.interpolate(point, scheme=scheme) else: - raise ValueError('Can only interpolate lat-lon grids if no coordinate system on cube') + raise ValueError('If no coordinate system on cube then can only interpolate lat-lon grids') else: # convert the target point(s) to lat lon and do the interpolation ll = ccrs.Geodetic() - trpoints = cube.coord_system().as_cartopy_crs().transform_point(point[1][1], point[0][1], ll) - trpoints = [(y_coord.name(), trpoints[1]),(x_coord.name(), trpoints[0])] + myccrs = cube.coord_system().as_cartopy_crs() + xpoints = np.array(point[1][1]) + ypoints = np.array(point[0][1]) + # repeat length 1 arrays to be the right shape + if xpoints.size == 1 and ypoints.size > 1: + xpoints = np.repeat(xpoints, ypoints.size) + if ypoints.size == 1 and xpoints.size > 1: + ypoints = np.repeat(ypoints, xpoints.size) + + try: + trpoints = myccrs.transform_points(ll, xpoints, ypoints) + except: + import pdb;pdb.set_trace() + trpoints = [(y_coord.name(), trpoints[:,1]),(x_coord.name(), trpoints[:,0])] return cube.interpolate(trpoints, scheme=scheme) def is_dataset(dataset): From 5739878c952e8c6e1877753fa6f2e680dc601a9b Mon Sep 17 00:00:00 2001 From: Emma Hogan Date: Wed, 30 Aug 2023 13:45:46 +0100 Subject: [PATCH 03/10] #2177: Update extract point tests to use a more realistic cubes that has coordinates with an associated coordinate system --- .../_regrid/test_extract_point.py | 53 +++++++++---------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/tests/unit/preprocessor/_regrid/test_extract_point.py b/tests/unit/preprocessor/_regrid/test_extract_point.py index 5c19530ff7..304e071b36 100644 --- a/tests/unit/preprocessor/_regrid/test_extract_point.py +++ b/tests/unit/preprocessor/_regrid/test_extract_point.py @@ -7,7 +7,7 @@ import unittest from unittest import mock -import iris +from iris.tests.stock import lat_lon_cube import tests from esmvalcore.preprocessor import extract_point @@ -15,24 +15,13 @@ class Test(tests.Test): + def setUp(self): - self.coord_system = mock.Mock(return_value=None) - self.coord = mock.sentinel.coord - self.coords = mock.Mock(return_value=[self.coord]) - self.remove_coord = mock.Mock() - self.point_cube = mock.sentinel.point_cube - self.interpolate = mock.Mock(return_value=self.point_cube) - self.src_cube = mock.Mock( - spec=iris.cube.Cube, - coord_system=self.coord_system, - coords=self.coords, - remove_coord=self.remove_coord, - interpolate=self.interpolate) - self.schemes = ['linear', 'nearest'] - - self.mocks = [ - self.coord_system, self.coords, self.interpolate, self.src_cube - ] + # Use an Iris test cube with coordinates that have a coordinate + # system, see the following issue for more details: + # https://github.com/ESMValGroup/ESMValCore/issues/2177. + self.src_cube = lat_lon_cube() + self.schemes = ["linear", "nearest"] def test_invalid_scheme__unknown(self): dummy = mock.sentinel.dummy @@ -41,20 +30,30 @@ def test_invalid_scheme__unknown(self): extract_point(dummy, dummy, dummy, 'non-existent') def test_interpolation_schemes(self): - self.assertEqual( - set(POINT_INTERPOLATION_SCHEMES.keys()), set(self.schemes)) + self.assertEqual(set(POINT_INTERPOLATION_SCHEMES.keys()), + set(self.schemes)) def test_extract_point_interpolation_schemes(self): - dummy = mock.sentinel.dummy + latitude = -90. + longitude = 0. for scheme in self.schemes: - result = extract_point(self.src_cube, dummy, dummy, scheme) - self.assertEqual(result, self.point_cube) + result = extract_point(self.src_cube, latitude, longitude, scheme) + self._assert_coords(result, latitude, longitude) def test_extract_point(self): - dummy = mock.sentinel.dummy - scheme = 'linear' - result = extract_point(self.src_cube, dummy, dummy, scheme) - self.assertEqual(result, self.point_cube) + latitude = 90. + longitude = -180. + for scheme in self.schemes: + result = extract_point(self.src_cube, latitude, longitude, scheme) + self._assert_coords(result, latitude, longitude) + + def _assert_coords(self, cube, ref_lat, ref_lon): + lat_points = cube.coord("latitude").points + lon_points = cube.coord("longitude").points + self.assertEqual(len(lat_points), 1) + self.assertEqual(len(lon_points), 1) + self.assertEqual(lat_points[0], ref_lat) + self.assertEqual(lon_points[0], ref_lon) if __name__ == '__main__': From 07c2e22abe09103ff31fb989803504deb92ff16c Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 30 Aug 2023 13:52:54 +0100 Subject: [PATCH 04/10] remove import added spuriously --- esmvalcore/preprocessor/_regrid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 85681fbddc..57abd4a849 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -9,7 +9,6 @@ from copy import deepcopy from decimal import Decimal from pathlib import Path -from turtle import pd from typing import Dict import dask.array as da From 3f543d534bac5ebde5e042cd83c9969062be4d7d Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 30 Aug 2023 14:12:05 +0100 Subject: [PATCH 05/10] remove debugging statement --- esmvalcore/preprocessor/_regrid.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 57abd4a849..70608293a4 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -516,10 +516,7 @@ def get_pt(cube, point, scheme): if ypoints.size == 1 and xpoints.size > 1: ypoints = np.repeat(ypoints, xpoints.size) - try: - trpoints = myccrs.transform_points(ll, xpoints, ypoints) - except: - import pdb;pdb.set_trace() + trpoints = myccrs.transform_points(ll, xpoints, ypoints) trpoints = [(y_coord.name(), trpoints[:,1]),(x_coord.name(), trpoints[:,0])] return cube.interpolate(trpoints, scheme=scheme) From 8f3c47eb5ac814c082b8dae003e057c35c7804d5 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 30 Aug 2023 15:20:26 +0100 Subject: [PATCH 06/10] fix issues flagged by pre-commit --- esmvalcore/preprocessor/_regrid.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index 70608293a4..c0a68d5cae 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -11,6 +11,7 @@ from pathlib import Path from typing import Dict +import cartopy.crs as ccrs import dask.array as da import iris import numpy as np @@ -18,7 +19,6 @@ from geopy.geocoders import Nominatim from iris.analysis import AreaWeighted, Linear, Nearest, UnstructuredNearest from iris.util import broadcast_to_shape -import cartopy.crs as ccrs from ..cmor._fixes.shared import add_altitude_from_plev, add_plev_from_altitude from ..cmor.table import CMOR_TABLES @@ -464,10 +464,10 @@ def extract_point(cube, latitude, longitude, scheme): cube = get_pt(cube, point, scheme) return cube + def get_pt(cube, point, scheme): - """ - Extract data at a single point from cubes with any coordinate - system or none (if lat-lon only) + """Extract data at a single point from cubes with any coordinate system or + none (if lat-lon only) Parameters ---------- @@ -492,7 +492,7 @@ def get_pt(cube, point, scheme): """ x_coord = cube.coord(axis='X', dim_coords=True) y_coord = cube.coord(axis='Y', dim_coords=True) - + # some lon-lat cubes don't have a coordinate system # check if it is lon-lat and if so just use interpolate # as it is @@ -503,7 +503,8 @@ def get_pt(cube, point, scheme): if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': return cube.interpolate(point, scheme=scheme) else: - raise ValueError('If no coordinate system on cube then can only interpolate lat-lon grids') + raise ValueError('If no coordinate system on cube then ' + + 'can only interpolate lat-lon grids') else: # convert the target point(s) to lat lon and do the interpolation ll = ccrs.Geodetic() @@ -515,11 +516,13 @@ def get_pt(cube, point, scheme): xpoints = np.repeat(xpoints, ypoints.size) if ypoints.size == 1 and xpoints.size > 1: ypoints = np.repeat(ypoints, xpoints.size) - - trpoints = myccrs.transform_points(ll, xpoints, ypoints) - trpoints = [(y_coord.name(), trpoints[:,1]),(x_coord.name(), trpoints[:,0])] + + trpoints = myccrs.transform_points(ll, xpoints, ypoints) + trpoints = [(y_coord.name(), trpoints[:, 1]), + (x_coord.name(), trpoints[:, 0])] return cube.interpolate(trpoints, scheme=scheme) + def is_dataset(dataset): """Test if something is an `esmvalcore.dataset.Dataset`.""" # Use this function to avoid circular imports @@ -617,7 +620,6 @@ def regrid(cube, target_grid, scheme, lat_offset=True, lon_offset=True): target: 1x1 scheme: reference: esmf_regrid.schemes:ESMFAreaWeighted - """ if is_dataset(target_grid): target_grid = target_grid.copy() From 7f0eec3d7e89824ab9377f1f726ace34bcd688b7 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 30 Aug 2023 15:52:55 +0100 Subject: [PATCH 07/10] fix issue with size of returned arrays by refactoring --- esmvalcore/preprocessor/_regrid.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index c0a68d5cae..da3e7e19c4 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -493,19 +493,16 @@ def get_pt(cube, point, scheme): x_coord = cube.coord(axis='X', dim_coords=True) y_coord = cube.coord(axis='Y', dim_coords=True) - # some lon-lat cubes don't have a coordinate system # check if it is lon-lat and if so just use interpolate - # as it is - if cube.coord_system() is None: - # if coords are equal to lat and lon then - # assume we have a lat-lon grid - - if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': - return cube.interpolate(point, scheme=scheme) - else: + # as it is (reproduce previous code) + if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': + return cube.interpolate(point, scheme=scheme) + + else: + if cube.coord_system() is None: raise ValueError('If no coordinate system on cube then ' + 'can only interpolate lat-lon grids') - else: + # convert the target point(s) to lat lon and do the interpolation ll = ccrs.Geodetic() myccrs = cube.coord_system().as_cartopy_crs() From 847ab1129419c5e5f775263acef27b005d92a38f Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Thu, 31 Aug 2023 10:27:27 +0100 Subject: [PATCH 08/10] make changes to pass tests in prospector --- esmvalcore/preprocessor/_regrid.py | 48 ++++++++++++++++-------------- 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py index da3e7e19c4..7caecf2c27 100644 --- a/esmvalcore/preprocessor/_regrid.py +++ b/esmvalcore/preprocessor/_regrid.py @@ -466,7 +466,9 @@ def extract_point(cube, latitude, longitude, scheme): def get_pt(cube, point, scheme): - """Extract data at a single point from cubes with any coordinate system or + """Extract data at a single point from cubes. + + Works for cubes with any coordinate system or none (if lat-lon only) Parameters @@ -474,8 +476,9 @@ def get_pt(cube, point, scheme): cube : cube The source cube to extract a point from. - latitude, longitude : float, or array of float - The latitude and longitude of the point. + point : list containing two tuples - one for + latitude and one for longitude e.g + [('latitude', latitude), ('longitude', longitude)] scheme : str The interpolation scheme. 'linear' or 'nearest'. No default. @@ -498,26 +501,25 @@ def get_pt(cube, point, scheme): if x_coord.name() == 'longitude' and y_coord.name() == 'latitude': return cube.interpolate(point, scheme=scheme) - else: - if cube.coord_system() is None: - raise ValueError('If no coordinate system on cube then ' + - 'can only interpolate lat-lon grids') - - # convert the target point(s) to lat lon and do the interpolation - ll = ccrs.Geodetic() - myccrs = cube.coord_system().as_cartopy_crs() - xpoints = np.array(point[1][1]) - ypoints = np.array(point[0][1]) - # repeat length 1 arrays to be the right shape - if xpoints.size == 1 and ypoints.size > 1: - xpoints = np.repeat(xpoints, ypoints.size) - if ypoints.size == 1 and xpoints.size > 1: - ypoints = np.repeat(ypoints, xpoints.size) - - trpoints = myccrs.transform_points(ll, xpoints, ypoints) - trpoints = [(y_coord.name(), trpoints[:, 1]), - (x_coord.name(), trpoints[:, 0])] - return cube.interpolate(trpoints, scheme=scheme) + if cube.coord_system() is None: + raise ValueError('If no coordinate system on cube then ' + + 'can only interpolate lat-lon grids') + + # convert the target point(s) to lat lon and do the interpolation + ll_cs = ccrs.Geodetic() + myccrs = cube.coord_system().as_cartopy_crs() + xpoints = np.array(point[1][1]) + ypoints = np.array(point[0][1]) + # repeat length 1 arrays to be the right shape + if xpoints.size == 1 and ypoints.size > 1: + xpoints = np.repeat(xpoints, ypoints.size) + if ypoints.size == 1 and xpoints.size > 1: + ypoints = np.repeat(ypoints, xpoints.size) + + trpoints = myccrs.transform_points(ll_cs, xpoints, ypoints) + trpoints = [(y_coord.name(), trpoints[:, 1]), + (x_coord.name(), trpoints[:, 0])] + return cube.interpolate(trpoints, scheme=scheme) def is_dataset(dataset): From 64e5e76b622c5d88fb15f347b2374ed1975334f2 Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Fri, 15 Sep 2023 14:34:40 +0100 Subject: [PATCH 09/10] add name and ORC ID to CITATIN.cff --- CITATION.cff | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CITATION.cff b/CITATION.cff index b8437333e0..f5bebc80b4 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -194,6 +194,11 @@ authors: affiliation: "Forschungszentrum Juelich (FZJ), Germany" family-names: Benke given-names: Joerg + - + affiliation: "MetOffice, UK" + family-names: Savage + given-names: Nicholas Henry + orcid: "https://orcid.org/0000-0001-9391-5100" cff-version: 1.2.0 date-released: 2023-07-04 From 5f254f74f9c7f8d714fe6395207216fe8fa2fb1b Mon Sep 17 00:00:00 2001 From: Nick Savage Date: Wed, 27 Sep 2023 18:18:16 +0100 Subject: [PATCH 10/10] have good framework for tests but results for r pole wrong. also Lambert just dummy --- .../_regrid/test_extract_point.py | 152 ++++++++++++++++++ 1 file changed, 152 insertions(+) diff --git a/tests/unit/preprocessor/_regrid/test_extract_point.py b/tests/unit/preprocessor/_regrid/test_extract_point.py index 304e071b36..3d744beb46 100644 --- a/tests/unit/preprocessor/_regrid/test_extract_point.py +++ b/tests/unit/preprocessor/_regrid/test_extract_point.py @@ -7,12 +7,20 @@ import unittest from unittest import mock +import iris.coords +import iris.cube +import numpy as np +from iris.analysis.cartography import get_xy_grids, unrotate_pole +from iris.coord_systems import GeogCS, RotatedGeogCS from iris.tests.stock import lat_lon_cube import tests from esmvalcore.preprocessor import extract_point from esmvalcore.preprocessor._regrid import POINT_INTERPOLATION_SCHEMES +# TODO: +# use these to test extract point + class Test(tests.Test): @@ -21,6 +29,10 @@ def setUp(self): # system, see the following issue for more details: # https://github.com/ESMValGroup/ESMValCore/issues/2177. self.src_cube = lat_lon_cube() + self.rpole_cube = _stock_rpole_2d() + self.lambert_cube = _stock_lambert_2d() + self.lambert_cube_nocs = _stock_lambert_2d_no_cs() + self.schemes = ["linear", "nearest"] def test_invalid_scheme__unknown(self): @@ -29,6 +41,16 @@ def test_invalid_scheme__unknown(self): with self.assertRaisesRegex(ValueError, emsg): extract_point(dummy, dummy, dummy, 'non-existent') + def test_invalid_coord_sys(self): + latitude = -90. + longitude = 0. + emsg = 'If no coordinate system on cube then ' + \ + 'can only interpolate lat-lon grids' + + with self.assertRaisesRegex(ValueError, emsg): + extract_point(self.lambert_cube_nocs, latitude, longitude, + self.schemes[0]) + def test_interpolation_schemes(self): self.assertEqual(set(POINT_INTERPOLATION_SCHEMES.keys()), set(self.schemes)) @@ -47,7 +69,25 @@ def test_extract_point(self): result = extract_point(self.src_cube, latitude, longitude, scheme) self._assert_coords(result, latitude, longitude) + def test_extract_point_rpole(self): + latitude = 90. + longitude = -180. + for scheme in self.schemes: + result = extract_point(self.rpole_cube, latitude, longitude, + scheme) + self._assert_coords_rpole(result, latitude, longitude) + + def test_extract_point_lambert(self): + latitude = 90. + longitude = -180. + for scheme in self.schemes: + result = extract_point(self.lambert_cube, latitude, longitude, + scheme) + self._assert_coords_lambert(result, latitude, longitude) + def _assert_coords(self, cube, ref_lat, ref_lon): + """For a 1D cube with a lat-lon coord system check that a 1x1 cube is + returned and the points are at the correct location.""" lat_points = cube.coord("latitude").points lon_points = cube.coord("longitude").points self.assertEqual(len(lat_points), 1) @@ -55,6 +95,118 @@ def _assert_coords(self, cube, ref_lat, ref_lon): self.assertEqual(lat_points[0], ref_lat) self.assertEqual(lon_points[0], ref_lon) + def _assert_coords_rpole(self, cube, ref_lat, ref_lon): + """For a 1D cube with a rotated coord system check that a 1x1 cube is + returned and the points are at the correct location Will need to + generate the lat and lon points from the grid using unrotate pole and + test with approx equal.""" + + pole_lat = cube.coord_system().grid_north_pole_latitude + pole_lon = cube.coord_system().grid_north_pole_longitude + rotated_lons, rotated_lats = get_xy_grids(cube) + tlons, tlats = unrotate_pole(rotated_lons, rotated_lats, pole_lon, + pole_lat) + self.assertEqual(len(tlats), 1) + self.assertEqual(len(tlons), 1) + self.assertEqual(tlats[0], ref_lat) + self.assertEqual(tlons[0], ref_lon) + + def _assert_coords_lambert(self, cube, ref_lat, ref_lon): + """For a 1D cube with a Lambert coord system check that a 1x1 cube is + returned and the points are at the correct location Will need to + generate the lat and lon points from the grid.""" + + pass + + +def _stock_rpole_2d(): + """Returns a realistic rotated pole 2d cube.""" + data = np.arange(9 * 11).reshape((9, 11)) + lat_pts = np.linspace(-4, 4, 9) + lon_pts = np.linspace(-5, 5, 11) + ll_cs = RotatedGeogCS(37.5, 177.5, ellipsoid=GeogCS(6371229.0)) + + lat = iris.coords.DimCoord( + lat_pts, + standard_name="grid_latitude", + units="degrees", + coord_system=ll_cs, + ) + lon = iris.coords.DimCoord( + lon_pts, + standard_name="grid_longitude", + units="degrees", + coord_system=ll_cs, + ) + cube = iris.cube.Cube( + data, + standard_name="air_potential_temperature", + units="K", + dim_coords_and_dims=[(lat, 0), (lon, 1)], + attributes={"source": "test cube"}, + ) + return cube + + +def _stock_lambert_2d(): + """Returns a realistic lambert conformal 2d cube.""" + data = np.arange(9 * 11).reshape((9, 11)) + y_pts = np.linspace(-96000., 96000., 9) + x_pts = np.linspace(0., 120000., 11) + lam_cs = iris.coord_systems.LambertConformal(central_lat=48.0, + central_lon=9.75, + false_easting=-6000.0, + false_northing=-6000.0, + secant_latitudes=(30.0, 65.0), + ellipsoid=GeogCS(6371229.0)) + + ydim = iris.coords.DimCoord( + y_pts, + standard_name="projection_y_coordinate", + units="m", + coord_system=lam_cs, + ) + xdim = iris.coords.DimCoord( + x_pts, + standard_name="projection_x_coordinate", + units="m", + coord_system=lam_cs, + ) + cube = iris.cube.Cube( + data, + standard_name="air_potential_temperature", + units="K", + dim_coords_and_dims=[(ydim, 0), (xdim, 1)], + attributes={"source": "test cube"}, + ) + return cube + + +def _stock_lambert_2d_no_cs(): + """Returns a realistic lambert conformal 2d cube.""" + data = np.arange(9 * 11).reshape((9, 11)) + y_pts = np.linspace(-96000., 96000., 9) + x_pts = np.linspace(0., 120000., 11) + + ydim = iris.coords.DimCoord( + y_pts, + standard_name="projection_y_coordinate", + units="m", + ) + xdim = iris.coords.DimCoord( + x_pts, + standard_name="projection_x_coordinate", + units="m", + ) + cube = iris.cube.Cube( + data, + standard_name="air_potential_temperature", + units="K", + dim_coords_and_dims=[(ydim, 0), (xdim, 1)], + attributes={"source": "test cube"}, + ) + return cube + if __name__ == '__main__': unittest.main()