Coverage for src/CSET/operators/regrid.py: 98%
89 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-10-30 15:17 +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 method: str
229 Method used to determine the values at the selected longitude and
230 latitude. The recommended approach is to use iris.analysis.Nearest(),
231 which selects the nearest gridpoint. An alternative is
232 iris.analysis.Linear(), which obtains the values at the selected
233 longitude and latitude by linear interpolation.
234 boundary_margin: int, optional
235 Number of grid points from the domain boundary considered "unreliable".
236 Defaults to 8.
238 Returns
239 -------
240 regridded_cubes: Cube | CubeList
241 An iris cube or CubeList of the data at the specified point (this may
242 have time and/or height dimensions).
244 Raises
245 ------
246 ValueError
247 If a unique x/y coordinate cannot be found; also if, for selecting a
248 single gridpoint, the chosen longitude and latitude point is outside the
249 domain.
250 NotImplementedError
251 If the cubes grid, or the method for regridding, is not yet supported.
253 Notes
254 -----
255 The acceptable coordinate names for X and Y coordinates are currently
256 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used
257 coordinate types, though a user can append new ones. Currently rectilinear
258 grids (uniform) are supported. Warnings are raised if the selected gridpoint
259 is within boundary_margin grid lengths of the domain boundary as data here
260 is potentially unreliable.
262 """
263 # To store regridded cubes.
264 regridded_cubes = iris.cube.CubeList()
266 # Iterate over all cubes and regrid.
267 for cube in iter_maybe(cubes):
268 # Get x and y coordinate names.
269 y_coord, x_coord = get_cube_yxcoordname(cube)
271 # List of supported grids - check if it is compatible
272 # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --
273 # this may need to be added in other places in these Operators.
274 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
275 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
276 raise NotImplementedError(
277 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
278 )
279 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
280 raise NotImplementedError(
281 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
282 )
284 # Transform input coordinates onto rotated grid if requested
285 if latlon_in_type == "realworld":
286 lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)
287 elif latlon_in_type == "rotated": 287 ↛ 291line 287 didn't jump to line 291 because the condition on line 287 was always true
288 lon_tr, lat_tr = lon_pt, lat_pt
290 # Get axis
291 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
293 # Get bounds
294 lat_min, lon_min = lat.points.min(), lon.points.min()
295 lat_max, lon_max = lat.points.max(), lon.points.max()
297 # Get boundaries of frame to avoid selecting gridpoint close to domain edge
298 lat_min_bound, lon_min_bound = (
299 lat.points[boundary_margin - 1],
300 lon.points[boundary_margin - 1],
301 )
302 lat_max_bound, lon_max_bound = (
303 lat.points[-boundary_margin],
304 lon.points[-boundary_margin],
305 )
307 # Check to see if selected point is outside the domain
308 if (lat_tr < lat_min) or (lat_tr > lat_max):
309 raise ValueError("Selected point is outside the domain.")
310 else:
311 if (lon_tr < lon_min) or (lon_tr > lon_max):
312 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
313 lon_tr += 360.0
314 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
315 lon_tr -= 360.0
316 else:
317 raise ValueError("Selected point is outside the domain.")
319 # Check to see if selected point is near the domain boundaries
320 if (
321 (lat_tr < lat_min_bound)
322 or (lat_tr > lat_max_bound)
323 or (lon_tr < lon_min_bound)
324 or (lon_tr > lon_max_bound)
325 ):
326 warnings.warn(
327 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
328 category=BoundaryWarning,
329 stacklevel=2,
330 )
332 regrid_method = getattr(iris.analysis, method, None)
333 if not callable(regrid_method):
334 raise NotImplementedError(
335 f"Does not currently support {method} regrid method"
336 )
337 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
338 regridded_cubes.append(cube_rgd)
339 # Preserve returning a cube if only a cube has been supplied to regrid.
340 if len(regridded_cubes) == 1: 340 ↛ 343line 340 didn't jump to line 343 because the condition on line 340 was always true
341 return regridded_cubes[0]
342 else:
343 return regridded_cubes
346def transform_lat_long_points(lon, lat, cube):
347 """Transform a selected point in longitude and latitude.
349 Transform the coordinates of a point from the real world
350 grid to the corresponding point on the rotated grid of a cube.
352 Parameters
353 ----------
354 cube: Cube
355 An iris cube of data defining the rotated grid to be used in
356 the longitude-latitude transformation.
357 lon: float
358 Selected value of longitude: this should be in the range -180 degrees to
359 180 degrees.
360 lat: float
361 Selected value of latitude: this should be in the range -90 degrees to
362 90 degrees.
364 Returns
365 -------
366 lon_rot, lat_rot: float
367 Coordinates of the selected point on the rotated grid specified within
368 the selected cube.
370 """
371 import cartopy.crs as ccrs
373 rot_pole = cube.coord_system().as_cartopy_crs()
374 true_grid = ccrs.Geodetic()
375 rot_coords = rot_pole.transform_point(lon, lat, true_grid)
376 lon_rot = rot_coords[0]
377 lat_rot = rot_coords[1]
379 return lon_rot, lat_rot