Coverage for src/CSET/operators/regrid.py: 99%
83 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-05 21:08 +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.cube
21import numpy as np
23from CSET._common import iter_maybe
24from CSET.operators._utils import get_cube_yxcoordname
27class BoundaryWarning(UserWarning):
28 """Selected gridpoint is close to the domain edge.
30 In many cases gridpoints near the domain boundary contain non-physical
31 values, so caution is advised when interpreting them.
32 """
35def regrid_onto_cube(
36 toregrid: iris.cube.Cube | iris.cube.CubeList,
37 target: iris.cube.Cube,
38 method: str,
39 **kwargs,
40) -> iris.cube.Cube | iris.cube.CubeList:
41 """Regrid a cube or CubeList, projecting onto a target cube.
43 All cubes must have at least 2 spatial (map) dimensions.
45 Arguments
46 ----------
47 toregrid: iris.cube | iris.cube.CubeList
48 An iris Cube of data to regrid, or multiple cubes to regrid in a
49 CubeList. A minimum requirement is that the cube(s) need to be 2D with a
50 latitude, longitude coordinates.
51 target: Cube
52 An iris cube of the data to regrid onto. It needs to be 2D with a
53 latitude, longitude coordinate.
54 method: str
55 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
57 Returns
58 -------
59 iris.cube | iris.cube.CubeList
60 An iris cube of the data that has been regridded, or a CubeList of the
61 cubes that have been regridded in the same order they were passed in
62 toregrid.
64 Raises
65 ------
66 ValueError
67 If a unique x/y coordinate cannot be found
68 NotImplementedError
69 If the cubes grid, or the method for regridding, is not yet supported.
71 Notes
72 -----
73 Currently rectlinear grids (uniform) are supported.
74 """
75 # To store regridded cubes.
76 regridded_cubes = iris.cube.CubeList()
78 # Iterate over all cubes and regrid.
79 for cube in iter_maybe(toregrid):
80 # Get y,x coord names
81 y_coord, x_coord = get_cube_yxcoordname(cube)
83 # List of supported grids - check if it is compatible
84 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
85 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
86 raise NotImplementedError(
87 f"Does not currently support {cube.coord(x_coord).coord_system} coordinate system"
88 )
89 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
90 raise NotImplementedError(
91 f"Does not currently support {cube.coord(y_coord).coord_system} coordinate system"
92 )
94 regrid_method = getattr(iris.analysis, method, None)
95 if callable(regrid_method):
96 regridded_cubes.append(cube.regrid(target, regrid_method()))
97 else:
98 raise NotImplementedError(
99 f"Does not currently support {method} regrid method"
100 )
102 # Preserve returning a cube if only a cube has been supplied to regrid.
103 if len(regridded_cubes) == 1:
104 return regridded_cubes[0]
105 else:
106 return regridded_cubes
109def regrid_onto_xyspacing(
110 toregrid: iris.cube.Cube | iris.cube.CubeList,
111 xspacing: float,
112 yspacing: float,
113 method: str,
114 **kwargs,
115) -> iris.cube.Cube | iris.cube.CubeList:
116 """Regrid cube or cubelist onto a set x,y spacing.
118 Regrid cube(s) using specified x,y spacing, which is performed linearly.
120 Parameters
121 ----------
122 toregrid: iris.cube | iris.cube.CubeList
123 An iris cube of the data to regrid, or multiple cubes to regrid in a
124 cubelist. A minimum requirement is that the cube(s) need to be 2D with a
125 latitude, longitude coordinates.
126 xspacing: float
127 Spacing of points in longitude direction (could be degrees, meters etc.)
128 yspacing: float
129 Spacing of points in latitude direction (could be degrees, meters etc.)
130 method: str
131 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()
133 Returns
134 -------
135 iris.cube | iris.cube.CubeList
136 An iris cube of the data that has been regridded, or a cubelist of the
137 cubes that have been regridded in the same order they were passed in
138 toregrid.
140 Raises
141 ------
142 ValueError
143 If a unique x/y coordinate cannot be found
144 NotImplementedError
145 If the cubes grid, or the method for regridding, is not yet supported.
147 Notes
148 -----
149 Currently rectlinear grids (uniform) are supported.
151 """
152 # To store regridded cubes.
153 regridded_cubes = iris.cube.CubeList()
155 # Iterate over all cubes and regrid.
156 for cube in iter_maybe(toregrid):
157 # Get x,y coord names
158 y_coord, x_coord = get_cube_yxcoordname(cube)
160 # List of supported grids - check if it is compatible
161 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
162 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
163 raise NotImplementedError(
164 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
165 )
166 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
167 raise NotImplementedError(
168 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
169 )
171 # Get axis
172 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
174 # Get bounds
175 lat_min, lon_min = lat.points.min(), lon.points.min()
176 lat_max, lon_max = lat.points.max(), lon.points.max()
178 # Generate new mesh
179 latout = np.arange(lat_min, lat_max, yspacing)
180 lonout = np.arange(lon_min, lon_max, xspacing)
182 regrid_method = getattr(iris.analysis, method, None)
183 if callable(regrid_method):
184 regridded_cubes.append(
185 cube.interpolate(
186 [(y_coord, latout), (x_coord, lonout)], regrid_method()
187 )
188 )
189 else:
190 raise NotImplementedError(
191 f"Does not currently support {method} regrid method"
192 )
194 # Preserve returning a cube if only a cube has been supplied to regrid.
195 if len(regridded_cubes) == 1:
196 return regridded_cubes[0]
197 else:
198 return regridded_cubes
201def regrid_to_single_point(
202 cube: iris.cube.Cube,
203 lat_pt: float,
204 lon_pt: float,
205 latlon_in_type: str = "rotated",
206 method: str = "Nearest",
207 boundary_margin: int = 8,
208 **kwargs,
209) -> iris.cube.Cube:
210 """Select data at a single point by longitude and latitude.
212 Selection of model grid point is performed by a regrid function, either
213 selecting the nearest gridpoint to the selected longitude and latitude
214 values or using linear interpolation across the surrounding points.
216 Parameters
217 ----------
218 cube: Cube
219 An iris cube of the data to regrid. As a minimum, it needs to be 2D with
220 latitude, longitude coordinates.
221 lon_pt: float
222 Selected value of longitude: this should be in the range -180 degrees to
223 180 degrees.
224 lat_pt: float
225 Selected value of latitude: this should be in the range -90 degrees to
226 90 degrees.
227 method: str
228 Method used to determine the values at the selected longitude and
229 latitude. The recommended approach is to use iris.analysis.Nearest(),
230 which selects the nearest gridpoint. An alternative is
231 iris.analysis.Linear(), which obtains the values at the selected
232 longitude and latitude by linear interpolation.
233 boundary_margin: int, optional
234 Number of grid points from the domain boundary considered "unreliable".
235 Defaults to 8.
237 Returns
238 -------
239 cube_rgd: Cube
240 An iris cube of the data at the specified point (this may have time
241 and/or height dimensions).
243 Raises
244 ------
245 ValueError
246 If a unique x/y coordinate cannot be found; also if, for selecting a
247 single gridpoint, the chosen longitude and latitude point is outside the
248 domain.
249 NotImplementedError
250 If the cubes grid, or the method for regridding, is not yet supported.
252 Notes
253 -----
254 The acceptable coordinate names for X and Y coordinates are currently
255 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used
256 coordinate types, though a user can append new ones. Currently rectilinear
257 grids (uniform) are supported. Warnings are raised if the selected gridpoint
258 is within boundary_margin grid lengths of the domain boundary as data here
259 is potentially unreliable.
261 """
262 # Get x and y coordinate names.
263 y_coord, x_coord = get_cube_yxcoordname(cube)
265 # List of supported grids - check if it is compatible
266 # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --
267 # this may need to be added in other places in these Operators.
268 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
269 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):
270 raise NotImplementedError(
271 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"
272 )
273 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):
274 raise NotImplementedError(
275 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
276 )
278 # Transform input coordinates onto rotated grid if requested
279 if latlon_in_type == "realworld":
280 lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)
281 elif latlon_in_type == "rotated": 281 ↛ 285line 281 didn't jump to line 285 because the condition on line 281 was always true
282 lon_tr, lat_tr = lon_pt, lat_pt
284 # Get axis
285 lat, lon = cube.coord(y_coord), cube.coord(x_coord)
287 # Get bounds
288 lat_min, lon_min = lat.points.min(), lon.points.min()
289 lat_max, lon_max = lat.points.max(), lon.points.max()
291 # Get boundaries of frame to avoid selecting gridpoint close to domain edge
292 lat_min_bound, lon_min_bound = (
293 lat.points[boundary_margin - 1],
294 lon.points[boundary_margin - 1],
295 )
296 lat_max_bound, lon_max_bound = (
297 lat.points[-boundary_margin],
298 lon.points[-boundary_margin],
299 )
301 # Check to see if selected point is outside the domain
302 if (lat_tr < lat_min) or (lat_tr > lat_max):
303 raise ValueError("Selected point is outside the domain.")
304 else:
305 if (lon_tr < lon_min) or (lon_tr > lon_max):
306 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
307 lon_tr += 360.0
308 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
309 lon_tr -= 360.0
310 else:
311 raise ValueError("Selected point is outside the domain.")
313 # Check to see if selected point is near the domain boundaries
314 if (
315 (lat_tr < lat_min_bound)
316 or (lat_tr > lat_max_bound)
317 or (lon_tr < lon_min_bound)
318 or (lon_tr > lon_max_bound)
319 ):
320 warnings.warn(
321 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
322 category=BoundaryWarning,
323 stacklevel=2,
324 )
326 regrid_method = getattr(iris.analysis, method, None)
327 if not callable(regrid_method):
328 raise NotImplementedError(f"Does not currently support {method} regrid method")
329 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
330 return cube_rgd
333def transform_lat_long_points(lon, lat, cube):
334 """Transform a selected point in longitude and latitude.
336 Transform the coordinates of a point from the real world
337 grid to the corresponding point on the rotated grid of a cube.
339 Parameters
340 ----------
341 cube: Cube
342 An iris cube of data defining the rotated grid to be used in
343 the longitude-latitude transformation.
344 lon: float
345 Selected value of longitude: this should be in the range -180 degrees to
346 180 degrees.
347 lat: float
348 Selected value of latitude: this should be in the range -90 degrees to
349 90 degrees.
351 Returns
352 -------
353 lon_rot, lat_rot: float
354 Coordinates of the selected point on the rotated grid specified within
355 the selected cube.
357 """
358 import cartopy.crs as ccrs
360 rot_pole = cube.coord_system().as_cartopy_crs()
361 true_grid = ccrs.Geodetic()
362 rot_coords = rot_pole.transform_point(lon, lat, true_grid)
363 lon_rot = rot_coords[0]
364 lat_rot = rot_coords[1]
366 return lon_rot, lat_rot