Coverage for src / CSET / operators / regrid.py: 99%
89 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-11-27 11:58 +0000
1# © Crown copyright, Met Office (2022-2024) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to regrid cubes."""
17import warnings
19import iris
20import iris.coord_systems
21import iris.cube
22import numpy as np
24from CSET._common import iter_maybe
25from CSET.operators._utils import get_cube_yxcoordname
28class BoundaryWarning(UserWarning):
29 """Selected gridpoint is close to the domain edge.
31 In many cases gridpoints near the domain boundary contain non-physical
32 values, so caution is advised when interpreting them.
33 """
36def regrid_onto_cube(
37 toregrid: iris.cube.Cube | iris.cube.CubeList,
38 target: iris.cube.Cube,
39 method: str,
40 **kwargs,
41) -> iris.cube.Cube | iris.cube.CubeList:
42 """Regrid a cube or CubeList, projecting onto a target cube.
44 All cubes must have at least 2 spatial (map) dimensions.
46 Arguments
47 ----------
48 toregrid: iris.cube | iris.cube.CubeList
49 An iris Cube of data to regrid, or multiple cubes to regrid in a
50 CubeList. A minimum requirement is that the cube(s) need to be 2D with a
51 latitude, longitude coordinates.
52 target: Cube
53 An iris cube of the data to regrid onto. It needs to be 2D with a
54 latitude, longitude coordinate.
55 method: str
56 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
58 Returns
59 -------
60 iris.cube | iris.cube.CubeList
61 An iris cube of the data that has been regridded, or a CubeList of the
62 cubes that have been regridded in the same order they were passed in
63 toregrid.
65 Raises
66 ------
67 ValueError
68 If a unique x/y coordinate cannot be found
69 NotImplementedError
70 If the cubes grid, or the method for regridding, is not yet supported.
72 Notes
73 -----
74 Currently rectlinear grids (uniform) are supported.
75 """
76 # To store regridded cubes.
77 regridded_cubes = iris.cube.CubeList()
79 # Iterate over all cubes and regrid.
80 for cube in iter_maybe(toregrid):
81 # Get y,x coord names
82 y_coord, x_coord = get_cube_yxcoordname(cube)
84 # List of supported grids - check if it is compatible
85 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
86 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
87 raise NotImplementedError(
88 f"Does not currently support {cube.coord(x_coord).coord_system} coordinate system"
89 )
90 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
91 raise NotImplementedError(
92 f"Does not currently support {cube.coord(y_coord).coord_system} coordinate system"
93 )
95 regrid_method = getattr(iris.analysis, method, None)
96 if callable(regrid_method):
97 regridded_cubes.append(cube.regrid(target, regrid_method()))
98 else:
99 raise NotImplementedError(
100 f"Does not currently support {method} regrid method"
101 )
103 # Preserve returning a cube if only a cube has been supplied to regrid.
104 if len(regridded_cubes) == 1:
105 return regridded_cubes[0]
106 else:
107 return regridded_cubes
110def regrid_onto_xyspacing(
111 toregrid: iris.cube.Cube | iris.cube.CubeList,
112 xspacing: float,
113 yspacing: float,
114 method: str,
115 **kwargs,
116) -> iris.cube.Cube | iris.cube.CubeList:
117 """Regrid cube or cubelist onto a set x,y spacing.
119 Regrid cube(s) using specified x,y spacing, which is performed linearly.
121 Parameters
122 ----------
123 toregrid: iris.cube | iris.cube.CubeList
124 An iris cube of the data to regrid, or multiple cubes to regrid in a
125 cubelist. A minimum requirement is that the cube(s) need to be 2D with a
126 latitude, longitude coordinates.
127 xspacing: float
128 Spacing of points in longitude direction (could be degrees, meters etc.)
129 yspacing: float
130 Spacing of points in latitude direction (could be degrees, meters etc.)
131 method: str
132 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
134 Returns
135 -------
136 iris.cube | iris.cube.CubeList
137 An iris cube of the data that has been regridded, or a cubelist of the
138 cubes that have been regridded in the same order they were passed in
139 toregrid.
141 Raises
142 ------
143 ValueError
144 If a unique x/y coordinate cannot be found
145 NotImplementedError
146 If the cubes grid, or the method for regridding, is not yet supported.
148 Notes
149 -----
150 Currently rectlinear grids (uniform) are supported.
152 """
153 # To store regridded cubes.
154 regridded_cubes = iris.cube.CubeList()
156 # Iterate over all cubes and regrid.
157 for cube in iter_maybe(toregrid):
158 # Get x,y coord names
159 y_coord, x_coord = get_cube_yxcoordname(cube)
161 # List of supported grids - check if it is compatible
162 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
163 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
164 raise NotImplementedError(
165 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
166 )
167 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
168 raise NotImplementedError(
169 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
170 )
172 # Get axis
173 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
175 # Get bounds
176 lat_min, lon_min = lat.points.min(), lon.points.min()
177 lat_max, lon_max = lat.points.max(), lon.points.max()
179 # Generate new mesh
180 latout = np.arange(lat_min, lat_max, yspacing)
181 lonout = np.arange(lon_min, lon_max, xspacing)
183 regrid_method = getattr(iris.analysis, method, None)
184 if callable(regrid_method):
185 regridded_cubes.append(
186 cube.interpolate(
187 [(y_coord, latout), (x_coord, lonout)], regrid_method()
188 )
189 )
190 else:
191 raise NotImplementedError(
192 f"Does not currently support {method} regrid method"
193 )
195 # Preserve returning a cube if only a cube has been supplied to regrid.
196 if len(regridded_cubes) == 1:
197 return regridded_cubes[0]
198 else:
199 return regridded_cubes
202def regrid_to_single_point(
203 cubes: iris.cube.Cube | iris.cube.CubeList,
204 lat_pt: float,
205 lon_pt: float,
206 latlon_in_type: str = "rotated",
207 method: str = "Nearest",
208 boundary_margin: int = 8,
209 **kwargs,
210) -> iris.cube.Cube:
211 """Select data at a single point by longitude and latitude.
213 Selection of model grid point is performed by a regrid function, either
214 selecting the nearest gridpoint to the selected longitude and latitude
215 values or using linear interpolation across the surrounding points.
217 Parameters
218 ----------
219 cubes: Cube | CubeList
220 An iris cube or CubeList of the data to regrid. As a minimum, it needs
221 to be 2D with latitude, longitude coordinates.
222 lon_pt: float
223 Selected value of longitude: this should be in the range -180 degrees to
224 180 degrees.
225 lat_pt: float
226 Selected value of latitude: this should be in the range -90 degrees to
227 90 degrees.
228 latlon_in_type: str, optional
229 Specify whether the input longitude and latitude point is in standard
230 geographic realworld coordinates ("realworld") or on the rotated grid
231 of the cube ("rotated"). Default is "rotated".
232 method: str
233 Method used to determine the values at the selected longitude and
234 latitude. The recommended approach is to use iris.analysis.Nearest(),
235 which selects the nearest gridpoint. An alternative is
236 iris.analysis.Linear(), which obtains the values at the selected
237 longitude and latitude by linear interpolation.
238 boundary_margin: int, optional
239 Number of grid points from the domain boundary considered "unreliable".
240 Defaults to 8.
242 Returns
243 -------
244 regridded_cubes: Cube | CubeList
245 An iris cube or CubeList of the data at the specified point (this may
246 have time and/or height dimensions).
248 Raises
249 ------
250 ValueError
251 If a unique x/y coordinate cannot be found; also if, for selecting a
252 single gridpoint, the chosen longitude and latitude point is outside the
253 domain.
254 NotImplementedError
255 If the cubes grid, or the method for regridding, is not yet supported.
257 Notes
258 -----
259 The acceptable coordinate names for X and Y coordinates are currently
260 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used
261 coordinate types, though a user can append new ones. Currently rectilinear
262 grids (uniform) are supported. Warnings are raised if the selected gridpoint
263 is within boundary_margin grid lengths of the domain boundary as data here
264 is potentially unreliable.
266 """
267 # To store regridded cubes.
268 regridded_cubes = iris.cube.CubeList()
270 # Iterate over all cubes and regrid.
271 for cube in iter_maybe(cubes):
272 # Get x and y coordinate names.
273 y_coord, x_coord = get_cube_yxcoordname(cube)
275 # List of supported grids - check if it is compatible
276 # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --
277 # this may need to be added in other places in these Operators.
278 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
279 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
280 raise NotImplementedError(
281 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
282 )
283 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
284 raise NotImplementedError(
285 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
286 )
288 # Transform input coordinates onto rotated grid if requested
289 if latlon_in_type == "realworld":
290 lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)
291 elif latlon_in_type == "rotated": 291 ↛ 295line 291 didn't jump to line 295 because the condition on line 291 was always true
292 lon_tr, lat_tr = lon_pt, lat_pt
294 # Get axis
295 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
297 # Get bounds
298 lat_min, lon_min = lat.points.min(), lon.points.min()
299 lat_max, lon_max = lat.points.max(), lon.points.max()
301 # Get boundaries of frame to avoid selecting gridpoint close to domain edge
302 lat_min_bound, lon_min_bound = (
303 lat.points[boundary_margin - 1],
304 lon.points[boundary_margin - 1],
305 )
306 lat_max_bound, lon_max_bound = (
307 lat.points[-boundary_margin],
308 lon.points[-boundary_margin],
309 )
311 # Check to see if selected point is outside the domain
312 if (lat_tr < lat_min) or (lat_tr > lat_max):
313 raise ValueError("Selected point is outside the domain.")
314 else:
315 if (lon_tr < lon_min) or (lon_tr > lon_max):
316 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
317 lon_tr += 360.0
318 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
319 lon_tr -= 360.0
320 else:
321 raise ValueError("Selected point is outside the domain.")
323 # Check to see if selected point is near the domain boundaries
324 if (
325 (lat_tr < lat_min_bound)
326 or (lat_tr > lat_max_bound)
327 or (lon_tr < lon_min_bound)
328 or (lon_tr > lon_max_bound)
329 ):
330 warnings.warn(
331 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
332 category=BoundaryWarning,
333 stacklevel=2,
334 )
336 regrid_method = getattr(iris.analysis, method, None)
337 if not callable(regrid_method):
338 raise NotImplementedError(
339 f"Does not currently support {method} regrid method"
340 )
341 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
342 regridded_cubes.append(cube_rgd)
343 # Preserve returning a cube if only a cube has been supplied to regrid.
344 if len(regridded_cubes) == 1:
345 return regridded_cubes[0]
346 else:
347 return regridded_cubes
350def transform_lat_long_points(lon, lat, cube):
351 """Transform a selected point in longitude and latitude.
353 Transform the coordinates of a point from the real world
354 grid to the corresponding point on the rotated grid of a cube.
356 Parameters
357 ----------
358 cube: Cube
359 An iris cube of data defining the rotated grid to be used in
360 the longitude-latitude transformation.
361 lon: float
362 Selected value of longitude: this should be in the range -180 degrees to
363 180 degrees.
364 lat: float
365 Selected value of latitude: this should be in the range -90 degrees to
366 90 degrees.
368 Returns
369 -------
370 lon_rot, lat_rot: float
371 Coordinates of the selected point on the rotated grid specified within
372 the selected cube.
374 """
375 import cartopy.crs as ccrs
377 rot_pole = cube.coord_system().as_cartopy_crs()
378 true_grid = ccrs.Geodetic()
379 rot_coords = rot_pole.transform_point(lon, lat, true_grid)
380 lon_rot = rot_coords[0]
381 lat_rot = rot_coords[1]
383 return lon_rot, lat_rot