Coverage for src / CSET / operators / regrid.py: 71%
135 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-08 17:20 +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 logging
18import warnings
20import iris
21import iris.coord_systems
22import iris.cube
23import numpy as np
25from CSET._common import iter_maybe
26from CSET.operators._utils import get_cube_yxcoordname
29class BoundaryWarning(UserWarning):
30 """Selected gridpoint is close to the domain edge.
32 In many cases gridpoints near the domain boundary contain non-physical
33 values, so caution is advised when interpreting them.
34 """
37def regrid_onto_cube(
38 toregrid: iris.cube.Cube | iris.cube.CubeList,
39 target: iris.cube.Cube,
40 method: str,
41 **kwargs,
42) -> iris.cube.Cube | iris.cube.CubeList:
43 """Regrid a cube or CubeList, projecting onto a target cube.
45 All cubes must have at least 2 spatial (map) dimensions.
47 Arguments
48 ----------
49 toregrid: iris.cube | iris.cube.CubeList
50 An iris Cube of data to regrid, or multiple cubes to regrid in a
51 CubeList. A minimum requirement is that the cube(s) need to be 2D with a
52 latitude, longitude coordinates.
53 target: Cube
54 An iris cube of the data to regrid onto. It needs to be 2D with a
55 latitude, longitude coordinate.
56 method: str
57 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
59 Returns
60 -------
61 iris.cube | iris.cube.CubeList
62 An iris cube of the data that has been regridded, or a CubeList of the
63 cubes that have been regridded in the same order they were passed in
64 toregrid.
66 Raises
67 ------
68 ValueError
69 If a unique x/y coordinate cannot be found
70 NotImplementedError
71 If the cubes grid, or the method for regridding, is not yet supported.
73 Notes
74 -----
75 Currently rectlinear grids (uniform) are supported.
76 """
77 # To store regridded cubes.
78 regridded_cubes = iris.cube.CubeList()
80 # Iterate over all cubes and regrid.
81 for cube in iter_maybe(toregrid):
82 # Get y,x coord names
83 y_coord, x_coord = get_cube_yxcoordname(cube)
85 # List of supported grids - check if it is compatible
86 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
87 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
88 raise NotImplementedError(
89 f"Does not currently support {cube.coord(x_coord).coord_system} coordinate system"
90 )
91 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
92 raise NotImplementedError(
93 f"Does not currently support {cube.coord(y_coord).coord_system} coordinate system"
94 )
96 regrid_method = getattr(iris.analysis, method, None)
97 if callable(regrid_method):
98 regridded_cubes.append(cube.regrid(target, regrid_method()))
99 else:
100 raise NotImplementedError(
101 f"Does not currently support {method} regrid method"
102 )
104 # Preserve returning a cube if only a cube has been supplied to regrid.
105 if len(regridded_cubes) == 1:
106 return regridded_cubes[0]
107 else:
108 return regridded_cubes
111def regrid_onto_xyspacing(
112 toregrid: iris.cube.Cube | iris.cube.CubeList,
113 xspacing: float,
114 yspacing: float,
115 method: str,
116 **kwargs,
117) -> iris.cube.Cube | iris.cube.CubeList:
118 """Regrid cube or cubelist onto a set x,y spacing.
120 Regrid cube(s) using specified x,y spacing, which is performed linearly.
122 Parameters
123 ----------
124 toregrid: iris.cube | iris.cube.CubeList
125 An iris cube of the data to regrid, or multiple cubes to regrid in a
126 cubelist. A minimum requirement is that the cube(s) need to be 2D with a
127 latitude, longitude coordinates.
128 xspacing: float
129 Spacing of points in longitude direction (could be degrees, meters etc.)
130 yspacing: float
131 Spacing of points in latitude direction (could be degrees, meters etc.)
132 method: str
133 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
135 Returns
136 -------
137 iris.cube | iris.cube.CubeList
138 An iris cube of the data that has been regridded, or a cubelist of the
139 cubes that have been regridded in the same order they were passed in
140 toregrid.
142 Raises
143 ------
144 ValueError
145 If a unique x/y coordinate cannot be found
146 NotImplementedError
147 If the cubes grid, or the method for regridding, is not yet supported.
149 Notes
150 -----
151 Currently rectlinear grids (uniform) are supported.
153 """
154 # To store regridded cubes.
155 regridded_cubes = iris.cube.CubeList()
157 # Iterate over all cubes and regrid.
158 for cube in iter_maybe(toregrid):
159 # Get x,y coord names
160 y_coord, x_coord = get_cube_yxcoordname(cube)
162 # List of supported grids - check if it is compatible
163 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
164 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
165 raise NotImplementedError(
166 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
167 )
168 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
169 raise NotImplementedError(
170 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
171 )
173 # Get axis
174 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
176 # Get bounds
177 lat_min, lon_min = lat.points.min(), lon.points.min()
178 lat_max, lon_max = lat.points.max(), lon.points.max()
180 # Generate new mesh
181 latout = np.arange(lat_min, lat_max, yspacing)
182 lonout = np.arange(lon_min, lon_max, xspacing)
184 regrid_method = getattr(iris.analysis, method, None)
185 if callable(regrid_method):
186 regridded_cubes.append(
187 cube.interpolate(
188 [(y_coord, latout), (x_coord, lonout)], regrid_method()
189 )
190 )
191 else:
192 raise NotImplementedError(
193 f"Does not currently support {method} regrid method"
194 )
196 # Preserve returning a cube if only a cube has been supplied to regrid.
197 if len(regridded_cubes) == 1:
198 return regridded_cubes[0]
199 else:
200 return regridded_cubes
203def regrid_to_single_point(
204 cubes: iris.cube.Cube | iris.cube.CubeList,
205 lat_pt: float,
206 lon_pt: float,
207 latlon_in_type: str = "rotated",
208 method: str = "Nearest",
209 boundary_margin: int = 8,
210 **kwargs,
211) -> iris.cube.Cube:
212 """Select data at a single point by longitude and latitude.
214 Selection of model grid point is performed by a regrid function, either
215 selecting the nearest gridpoint to the selected longitude and latitude
216 values or using linear interpolation across the surrounding points.
218 Parameters
219 ----------
220 cubes: Cube | CubeList
221 An iris cube or CubeList of the data to regrid. As a minimum, it needs
222 to be 2D with latitude, longitude coordinates.
223 lon_pt: float
224 Selected value of longitude: this should be in the range -180 degrees to
225 180 degrees.
226 lat_pt: float
227 Selected value of latitude: this should be in the range -90 degrees to
228 90 degrees.
229 latlon_in_type: str, optional
230 Specify whether the input longitude and latitude point is in standard
231 geographic realworld coordinates ("realworld") or on the rotated grid
232 of the cube ("rotated"). Default is "rotated".
233 method: str
234 Method used to determine the values at the selected longitude and
235 latitude. The recommended approach is to use iris.analysis.Nearest(),
236 which selects the nearest gridpoint. An alternative is
237 iris.analysis.Linear(), which obtains the values at the selected
238 longitude and latitude by linear interpolation.
239 boundary_margin: int, optional
240 Number of grid points from the domain boundary considered "unreliable".
241 Defaults to 8.
243 Returns
244 -------
245 regridded_cubes: Cube | CubeList
246 An iris cube or CubeList of the data at the specified point (this may
247 have time and/or height dimensions).
249 Raises
250 ------
251 ValueError
252 If a unique x/y coordinate cannot be found; also if, for selecting a
253 single gridpoint, the chosen longitude and latitude point is outside the
254 domain; also (currently) if the difference between the actual and target
255 points exceed 0.1 degrees.
256 NotImplementedError
257 If the cubes grid, or the method for regridding, is not yet supported.
259 Notes
260 -----
261 The acceptable coordinate names for X and Y coordinates are currently
262 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used
263 coordinate types, though a user can append new ones. Currently rectilinear
264 grids (uniform) are supported. Warnings are raised if the selected gridpoint
265 is within boundary_margin grid lengths of the domain boundary as data here
266 is potentially unreliable.
267 """
268 # To store regridded cubes.
269 regridded_cubes = iris.cube.CubeList()
271 # Iterate over all cubes and regrid.
272 for cube in iter_maybe(cubes):
273 # Get x and y coordinate names.
274 y_coord, x_coord = get_cube_yxcoordname(cube)
276 # List of supported grids - check if it is compatible
277 # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --
278 # this may need to be added in other places in these Operators.
279 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
280 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
281 raise NotImplementedError(
282 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
283 )
284 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
285 raise NotImplementedError(
286 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
287 )
289 # Transform input coordinates onto rotated grid if requested
290 if latlon_in_type == "realworld":
291 lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)
292 elif latlon_in_type == "rotated": 292 ↛ 296line 292 didn't jump to line 296 because the condition on line 292 was always true
293 lon_tr, lat_tr = lon_pt, lat_pt
295 # Get axis
296 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
298 # Get bounds
299 lat_min, lon_min = lat.points.min(), lon.points.min()
300 lat_max, lon_max = lat.points.max(), lon.points.max()
302 # Use different logic for single point obs data.
303 if len(cube.coord(x_coord).points) > 1:
304 # Get boundaries of frame to avoid selecting gridpoint close to domain edge
305 lat_min_bound, lon_min_bound = (
306 lat.points[boundary_margin - 1],
307 lon.points[boundary_margin - 1],
308 )
309 lat_max_bound, lon_max_bound = (
310 lat.points[-boundary_margin],
311 lon.points[-boundary_margin],
312 )
314 # Check to see if selected point is outside the domain
315 if (lat_tr < lat_min) or (lat_tr > lat_max):
316 raise ValueError("Selected point is outside the domain.")
317 else:
318 if (lon_tr < lon_min) or (lon_tr > lon_max):
319 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
320 lon_tr += 360.0
321 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
322 lon_tr -= 360.0
323 else:
324 raise ValueError("Selected point is outside the domain.")
326 # Check to see if selected point is near the domain boundaries
327 if (
328 (lat_tr < lat_min_bound)
329 or (lat_tr > lat_max_bound)
330 or (lon_tr < lon_min_bound)
331 or (lon_tr > lon_max_bound)
332 ):
333 warnings.warn(
334 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
335 category=BoundaryWarning,
336 stacklevel=2,
337 )
339 regrid_method = getattr(iris.analysis, method, None)
340 if not callable(regrid_method):
341 raise NotImplementedError(
342 f"Does not currently support {method} regrid method"
343 )
345 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
346 regridded_cubes.append(cube_rgd)
347 else:
348 if (
349 np.abs((lat_tr - lat.points[0])) > 0.1
350 or np.abs((lon_tr - lon.points[0])) > 0.1
351 ):
352 raise ValueError(
353 "Selected point is too far from the specified coordinates. It should be within 0.1 degrees."
354 )
355 else:
356 print(
357 "*** lat/long diffs",
358 np.abs(lat_tr - lat_pt),
359 np.abs(lon_tr - lon_pt),
360 )
361 regridded_cubes.append(cube)
363 # Preserve returning a cube if only a cube has been supplied to regrid.
364 if len(regridded_cubes) == 1:
365 return regridded_cubes[0]
366 else:
367 return regridded_cubes
370def transform_lat_long_points(lon, lat, cube):
371 """Transform a selected point in longitude and latitude.
373 Transform the coordinates of a point from the real world
374 grid to the corresponding point on the rotated grid of a cube.
376 Parameters
377 ----------
378 cube: Cube
379 An iris cube of data defining the rotated grid to be used in
380 the longitude-latitude transformation.
381 lon: float
382 Selected value of longitude: this should be in the range -180 degrees to
383 180 degrees.
384 lat: float
385 Selected value of latitude: this should be in the range -90 degrees to
386 90 degrees.
388 Returns
389 -------
390 lon_rot, lat_rot: float
391 Coordinates of the selected point on the rotated grid specified within
392 the selected cube.
394 """
395 import cartopy.crs as ccrs
397 rot_pole = cube.coord_system().as_cartopy_crs()
398 true_grid = ccrs.Geodetic()
399 rot_coords = rot_pole.transform_point(lon, lat, true_grid)
400 lon_rot = rot_coords[0]
401 lat_rot = rot_coords[1]
403 return lon_rot, lat_rot
406def interpolate_to_point_cube(
407 fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs
408) -> iris.cube.Cube | iris.cube.CubeList:
409 """Interpolate a 2D field in cube or CubeList to a set of points.
411 Regrid cube(s) to set of sample points specified by point_cube.
412 Ensures only matching times between input and point_cube are included.
414 Parameters
415 ----------
416 fld: Cube | CubeList
417 An iris cube or CubeList containing a two-dimensional field(s).
418 point_cube: Cube
419 An iris cube specifying the coord point(s) to which the data
420 will be interpolated.
422 Returns
423 -------
424 fld_point_cube: Cube | CubeList
425 An iris cube or CubeList containing interpolated values at the
426 points specified by the point cube for matching times common to
427 both fld and point_cube.
428 """
429 # Empty CubeList To store regridded cubes.
430 regridded_cubes = iris.cube.CubeList()
432 # Iterate over all cubes and regrid.
433 for cube in iter_maybe(fld):
434 # Ensure matching times in fld cube and point_cube
435 base_time_coord = point_cube.coord("time")
436 other_time_coord = cube.coord("time")
437 base_times = base_time_coord.units.num2date(base_time_coord.points)
438 other_times = other_time_coord.units.num2date(other_time_coord.points)
439 shared_times = set.intersection(set(base_times), set(other_times))
440 logging.debug("Shared times: %s", shared_times)
441 time_constraint = iris.Constraint(
442 coord_values={
443 "time": lambda cell, shared_times=shared_times: (
444 cell.point in shared_times
445 )
446 }
447 )
449 # Extract points matching the shared times.
450 cube = cube.extract(time_constraint)
451 point_cube = point_cube.extract(time_constraint)
452 if cube is None or point_cube is None:
453 raise ValueError("No common time points found!")
455 # Generate array of point cube lat and lon points.
456 lat_name, lon_name = get_cube_yxcoordname(point_cube)
457 lats = point_cube.coord(lat_name).points
458 lons = point_cube.coord(lon_name).points
459 if lat_name == "latitude":
460 sample_points = [
461 ("latitude", np.array(lats)),
462 ("longitude", np.array(lons)),
463 ]
464 else:
465 sample_points = [("grid_latitude", lats), ("grid_longitude", lons)]
467 # Interpolate fld cube to required sample points
468 fld_point_cube = cube.interpolate(
469 sample_points, iris.analysis.Linear(extrapolation_mode="mask")
470 )
472 # Retain only diagonal elements of 2D interpolated cube to vector points
473 od_index = point_cube.coord_dims("station")[0]
474 fv_cube = iris.cube.Cube(
475 fld_point_cube.data.diagonal(axis1=od_index, axis2=od_index + 1),
476 standard_name=cube.standard_name,
477 long_name=cube.long_name,
478 units=cube.units,
479 )
480 # Copy all non-lat/lon coordinates and cube attributes
481 if "time" in [coord.name() for coord in fld_point_cube.coords(dim_coords=True)]:
482 fv_cube.add_dim_coord(point_cube.coord("time"), 0)
483 else:
484 fv_cube.add_aux_coord(point_cube.coord("time"), od_index)
485 for coord in cube.coords():
486 if coord.name() not in [
487 "time",
488 "latitude",
489 "longitude",
490 "grid_latitude",
491 "grid_longitude",
492 ]:
493 fv_cube.add_aux_coord(coord.copy(), cube.coord_dims(coord))
494 for coord in point_cube.coords():
495 if coord.name() not in ["time", "realization", "station"]:
496 fv_cube.add_aux_coord(coord.copy(), point_cube.coord_dims(coord))
497 fv_cube.add_dim_coord(point_cube.coord("station"), od_index)
498 fv_cube.attributes = cube.attributes.copy()
499 fv_cube.cell_methods = cube.cell_methods
500 fv_cube.units = cube.units
501 regridded_cubes.append(fv_cube)
503 # Preserve returning a cube if only a cube has been supplied to regrid.
504 if len(regridded_cubes) == 1:
505 return regridded_cubes[0]
506 else:
507 return regridded_cubes