Coverage for src / CSET / operators / regrid.py: 82%
113 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-01 14:49 +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; also (currently) if the difference between the actual and target
254 points exceed 0.1 degrees.
255 NotImplementedError
256 If the cubes grid, or the method for regridding, is not yet supported.
258 Notes
259 -----
260 The acceptable coordinate names for X and Y coordinates are currently
261 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used
262 coordinate types, though a user can append new ones. Currently rectilinear
263 grids (uniform) are supported. Warnings are raised if the selected gridpoint
264 is within boundary_margin grid lengths of the domain boundary as data here
265 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 # Use different logic for single point obs data.
302 if len(cube.coord(x_coord).points) > 1:
303 # Get boundaries of frame to avoid selecting gridpoint close to domain edge
304 lat_min_bound, lon_min_bound = (
305 lat.points[boundary_margin - 1],
306 lon.points[boundary_margin - 1],
307 )
308 lat_max_bound, lon_max_bound = (
309 lat.points[-boundary_margin],
310 lon.points[-boundary_margin],
311 )
313 # Check to see if selected point is outside the domain
314 if (lat_tr < lat_min) or (lat_tr > lat_max):
315 raise ValueError("Selected point is outside the domain.")
316 else:
317 if (lon_tr < lon_min) or (lon_tr > lon_max):
318 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
319 lon_tr += 360.0
320 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
321 lon_tr -= 360.0
322 else:
323 raise ValueError("Selected point is outside the domain.")
325 # Check to see if selected point is near the domain boundaries
326 if (
327 (lat_tr < lat_min_bound)
328 or (lat_tr > lat_max_bound)
329 or (lon_tr < lon_min_bound)
330 or (lon_tr > lon_max_bound)
331 ):
332 warnings.warn(
333 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
334 category=BoundaryWarning,
335 stacklevel=2,
336 )
338 regrid_method = getattr(iris.analysis, method, None)
339 if not callable(regrid_method):
340 raise NotImplementedError(
341 f"Does not currently support {method} regrid method"
342 )
344 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
345 regridded_cubes.append(cube_rgd)
346 else:
347 if (
348 np.abs((lat_tr - lat.points[0])) > 0.1
349 or np.abs((lon_tr - lon.points[0])) > 0.1
350 ):
351 raise ValueError(
352 "Selected point is too far from the specified coordinates. It should be within 0.1 degrees."
353 )
354 else:
355 print(
356 "*** lat/long diffs",
357 np.abs(lat_tr - lat_pt),
358 np.abs(lon_tr - lon_pt),
359 )
360 regridded_cubes.append(cube)
362 # Preserve returning a cube if only a cube has been supplied to regrid.
363 if len(regridded_cubes) == 1:
364 return regridded_cubes[0]
365 else:
366 return regridded_cubes
369def transform_lat_long_points(lon, lat, cube):
370 """Transform a selected point in longitude and latitude.
372 Transform the coordinates of a point from the real world
373 grid to the corresponding point on the rotated grid of a cube.
375 Parameters
376 ----------
377 cube: Cube
378 An iris cube of data defining the rotated grid to be used in
379 the longitude-latitude transformation.
380 lon: float
381 Selected value of longitude: this should be in the range -180 degrees to
382 180 degrees.
383 lat: float
384 Selected value of latitude: this should be in the range -90 degrees to
385 90 degrees.
387 Returns
388 -------
389 lon_rot, lat_rot: float
390 Coordinates of the selected point on the rotated grid specified within
391 the selected cube.
393 """
394 import cartopy.crs as ccrs
396 rot_pole = cube.coord_system().as_cartopy_crs()
397 true_grid = ccrs.Geodetic()
398 rot_coords = rot_pole.transform_point(lon, lat, true_grid)
399 lon_rot = rot_coords[0]
400 lat_rot = rot_coords[1]
402 return lon_rot, lat_rot
405def interpolate_to_point_cube(
406 fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs
407) -> iris.cube.Cube | iris.cube.CubeList:
408 """Interpolate from a 2D field to a set of points.
410 Interpolate the 2D field in fld to the set of points
411 specified in point_cube.
413 Parameters
414 ----------
415 fld: Cube
416 An iris cube containing a two-dimensional field.
417 point_cube: Cube
418 An iris cube specifying the point to which the data
419 will be interpolated.
421 Returns
422 -------
423 fld_point_cube: Cube
424 An iris cube containing interpolated values at the points
425 specified by the point cube.
427 """
428 #
429 # As a basis, create a copy of the point cube.
430 fld_point_cube = point_cube.copy()
431 # Get indices of positional coordinates. We assume that the
432 # point cube is unrotated.
433 klon = None
434 klat = None
435 for kc in range(len(fld_point_cube.aux_coords)):
436 if fld_point_cube.aux_coords[kc].standard_name == "latitude":
437 klat = kc
438 elif fld_point_cube.aux_coords[kc].standard_name == "longitude":
439 klon = kc
440 #
441 # The input may have a rotated coordinate system.
442 if len(fld.coords("grid_latitude")) > 0:
443 # Interpolate in rotated coordinates.
444 rot_csyst = fld.coords("grid_latitude")[0].coord_system
445 rotpt = iris.analysis.cartography.rotate_pole(
446 fld_point_cube.aux_coords[klon].points,
447 fld_point_cube.aux_coords[klat].points,
448 rot_csyst.grid_north_pole_longitude,
449 rot_csyst.grid_north_pole_latitude,
450 )
451 # Add other interpolation options later.
452 fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator(
453 fld, ["time", "grid_latitude", "grid_longitude"]
454 )
455 for jt in range(len(fld_point_cube.coords("time")[0].points)):
456 fld_point_cube.data[jt, :] = np.ma.masked_invalid(
457 [
458 fld_interpolator(
459 [
460 fld_point_cube.coord("time").points[jt],
461 rotpt[1][k],
462 rotpt[0][k],
463 ]
464 ).data
465 if ~point_cube.data.mask[jt][k]
466 else np.nan
467 for k in range(len(rotpt[0]))
468 ]
469 )
470 else:
471 # Add other interpolation options later.
472 fld_interpolator = iris.analysis.Linear(extrapolation_mode="mask").interpolator(
473 fld, ["time", "latitude", "longitude"]
474 )
475 for jt in range(len(fld_point_cube.coords("time")[0].points)):
476 fld_point_cube.data[jt, :] = np.ma.masked_invalid(
477 [
478 fld_interpolator(
479 [
480 fld_point_cube.coords("time")[0].points[jt],
481 fld_point_cube.coord("latitude").points[k],
482 fld_point_cube.coord("longitude").points[k],
483 ]
484 ).data
485 if ~point_cube.data.mask[jt][k]
486 else np.nan
487 for k in range(fld_point_cube.coord("latitude").points)
488 ]
489 )
490 return fld_point_cube