Coverage for src/CSET/operators/regrid.py: 97%

1# © Crown copyright, Met Office (2022-2024) and CSET contributors.

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

7# http://www.apache.org/licenses/LICENSE-2.0

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

24from iris.analysis.cartography import rotate_pole

26from CSET._common import iter_maybe

27from CSET.operators._utils import get_cube_yxcoordname

30class BoundaryWarning(UserWarning):

31 """Selected gridpoint is close to the domain edge.

33 In many cases gridpoints near the domain boundary contain non-physical

34 values, so caution is advised when interpreting them.

35 """

38def regrid_onto_cube(

39 toregrid: iris.cube.Cube | iris.cube.CubeList,

40 target: iris.cube.Cube,

41 method: str,

42 **kwargs,

43) -> iris.cube.Cube | iris.cube.CubeList:

44 """Regrid a cube or CubeList, projecting onto a target cube.

46 All cubes must have at least 2 spatial (map) dimensions.

48 Arguments

49 ----------

50 toregrid: iris.cube | iris.cube.CubeList

51 An iris Cube of data to regrid, or multiple cubes to regrid in a

52 CubeList. A minimum requirement is that the cube(s) need to be 2D with a

53 latitude, longitude coordinates.

54 target: Cube

55 An iris cube of the data to regrid onto. It needs to be 2D with a

56 latitude, longitude coordinate.

57 method: str

58 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()

60 Returns

61 -------

62 iris.cube | iris.cube.CubeList

63 An iris cube of the data that has been regridded, or a CubeList of the

64 cubes that have been regridded in the same order they were passed in

65 toregrid.

67 Raises

68 ------

69 ValueError

70 If a unique x/y coordinate cannot be found

71 NotImplementedError

72 If the cubes grid, or the method for regridding, is not yet supported.

74 Notes

75 -----

76 Currently rectlinear grids (uniform) are supported.

77 """

78 # To store regridded cubes.

79 regridded_cubes = iris.cube.CubeList()

81 # Iterate over all cubes and regrid.

82 for cube in iter_maybe(toregrid):

83 # Get y,x coord names

84 y_coord, x_coord = get_cube_yxcoordname(cube)

86 # List of supported grids - check if it is compatible

87 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)

88 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):

89 raise NotImplementedError(

90 f"Does not currently support {cube.coord(x_coord).coord_system} coordinate system"

91 )

92 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):

93 raise NotImplementedError(

94 f"Does not currently support {cube.coord(y_coord).coord_system} coordinate system"

95 )

97 regrid_method = getattr(iris.analysis, method, None)

98 if callable(regrid_method):

99 regridded_cubes.append(cube.regrid(target, regrid_method()))

100 else:

101 raise NotImplementedError(

102 f"Does not currently support {method} regrid method"

103 )

104

105 # Preserve returning a cube if only a cube has been supplied to regrid.

106 if len(regridded_cubes) == 1:

107 return regridded_cubes[0]

108 else:

109 return regridded_cubes

110

111

112def regrid_onto_xyspacing(

113 toregrid: iris.cube.Cube | iris.cube.CubeList,

114 xspacing: float,

115 yspacing: float,

116 method: str,

117 **kwargs,

118) -> iris.cube.Cube | iris.cube.CubeList:

119 """Regrid cube or cubelist onto a set x,y spacing.

120

121 Regrid cube(s) using specified x,y spacing, which is performed linearly.

122

123 Parameters

124 ----------

125 toregrid: iris.cube | iris.cube.CubeList

126 An iris cube of the data to regrid, or multiple cubes to regrid in a

127 cubelist. A minimum requirement is that the cube(s) need to be 2D with a

128 latitude, longitude coordinates.

129 xspacing: float

130 Spacing of points in longitude direction (could be degrees, meters etc.)

131 yspacing: float

132 Spacing of points in latitude direction (could be degrees, meters etc.)

133 method: str

134 Method used to regrid onto, etc. Linear will use iris.analysis.Linear()

135

136 Returns

137 -------

138 iris.cube | iris.cube.CubeList

139 An iris cube of the data that has been regridded, or a cubelist of the

140 cubes that have been regridded in the same order they were passed in

141 toregrid.

142

143 Raises

144 ------

145 ValueError

146 If a unique x/y coordinate cannot be found

147 NotImplementedError

148 If the cubes grid, or the method for regridding, is not yet supported.

149

150 Notes

151 -----

152 Currently rectlinear grids (uniform) are supported.

153

154 """

155 # To store regridded cubes.

156 regridded_cubes = iris.cube.CubeList()

157

158 # Iterate over all cubes and regrid.

159 for cube in iter_maybe(toregrid):

160 # Get x,y coord names

161 y_coord, x_coord = get_cube_yxcoordname(cube)

162

163 # List of supported grids - check if it is compatible

164 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)

165 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):

166 raise NotImplementedError(

167 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"

168 )

169 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):

170 raise NotImplementedError(

171 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"

172 )

173

174 # Get axis

175 lat, lon = cube.coord(y_coord), cube.coord(x_coord)

176

177 # Get bounds

178 lat_min, lon_min = lat.points.min(), lon.points.min()

179 lat_max, lon_max = lat.points.max(), lon.points.max()

180

181 # Generate new mesh

182 latout = np.arange(lat_min, lat_max, yspacing)

183 lonout = np.arange(lon_min, lon_max, xspacing)

184

185 regrid_method = getattr(iris.analysis, method, None)

186 if callable(regrid_method):

187 regridded_cubes.append(

188 cube.interpolate(

189 [(y_coord, latout), (x_coord, lonout)], regrid_method()

190 )

191 )

192 else:

193 raise NotImplementedError(

194 f"Does not currently support {method} regrid method"

195 )

196

197 # Preserve returning a cube if only a cube has been supplied to regrid.

198 if len(regridded_cubes) == 1:

199 return regridded_cubes[0]

200 else:

201 return regridded_cubes

202

203

204def regrid_to_single_point(

205 cubes: iris.cube.Cube | iris.cube.CubeList,

206 lat_pt: float,

207 lon_pt: float,

208 latlon_in_type: str = "rotated",

209 method: str = "Nearest",

210 boundary_margin: int = 8,

211 **kwargs,

212) -> iris.cube.Cube:

213 """Select data at a single point by longitude and latitude.

214

215 Selection of model grid point is performed by a regrid function, either

216 selecting the nearest gridpoint to the selected longitude and latitude

217 values or using linear interpolation across the surrounding points.

218

219 Parameters

220 ----------

221 cubes: Cube | CubeList

222 An iris cube or CubeList of the data to regrid. As a minimum, it needs

223 to be 2D with latitude, longitude coordinates.

224 lon_pt: float

225 Selected value of longitude: this should be in the range -180 degrees to

226 180 degrees.

227 lat_pt: float

228 Selected value of latitude: this should be in the range -90 degrees to

229 90 degrees.

230 latlon_in_type: str, optional

231 Specify whether the input longitude and latitude point is in standard

232 geographic realworld coordinates ("realworld") or on the rotated grid

233 of the cube ("rotated"). Default is "rotated".

234 method: str

235 Method used to determine the values at the selected longitude and

236 latitude. The recommended approach is to use iris.analysis.Nearest(),

237 which selects the nearest gridpoint. An alternative is

238 iris.analysis.Linear(), which obtains the values at the selected

239 longitude and latitude by linear interpolation.

240 boundary_margin: int, optional

241 Number of grid points from the domain boundary considered "unreliable".

242 Defaults to 8.

243

244 Returns

245 -------

246 regridded_cubes: Cube | CubeList

247 An iris cube or CubeList of the data at the specified point (this may

248 have time and/or height dimensions).

249

250 Raises

251 ------

252 ValueError

253 If a unique x/y coordinate cannot be found; also if, for selecting a

254 single gridpoint, the chosen longitude and latitude point is outside the

255 domain; also (currently) if the difference between the actual and target

256 points exceed 0.1 degrees.

257 NotImplementedError

258 If the cubes grid, or the method for regridding, is not yet supported.

259

260 Notes

261 -----

262 The acceptable coordinate names for X and Y coordinates are currently

263 described in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used

264 coordinate types, though a user can append new ones. Currently rectilinear

265 grids (uniform) are supported. Warnings are raised if the selected gridpoint

266 is within boundary_margin grid lengths of the domain boundary as data here

267 is potentially unreliable.

268 """

269 # To store regridded cubes.

270 regridded_cubes = iris.cube.CubeList()

271

272 # Iterate over all cubes and regrid.

273 for cube in iter_maybe(cubes):

274 # Get x and y coordinate names.

275 y_coord, x_coord = get_cube_yxcoordname(cube)

276

277 # List of supported grids - check if it is compatible

278 # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --

279 # this may need to be added in other places in these Operators.

280 supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)

281 if not isinstance(cube.coord(x_coord).coord_system, supported_grids):

282 raise NotImplementedError(

283 f"Does not currently support {cube.coord(x_coord).coord_system} regrid method"

284 )

285 if not isinstance(cube.coord(y_coord).coord_system, supported_grids):

286 raise NotImplementedError(

287 f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"

288 )

289

290 # Transform input coordinates onto rotated grid if requested

291 if latlon_in_type == "realworld":

292 lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)

293 elif latlon_in_type == "rotated": 293 ↛ 297line 293 didn't jump to line 297 because the condition on line 293 was always true

294 lon_tr, lat_tr = lon_pt, lat_pt

295

296 # Get axis

297 lat, lon = cube.coord(y_coord), cube.coord(x_coord)

298

299 # Get bounds

300 lat_min, lon_min = lat.points.min(), lon.points.min()

301 lat_max, lon_max = lat.points.max(), lon.points.max()

302

303 # Use different logic for single point obs data.

304 if len(cube.coord(x_coord).points) > 1:

305 # Get boundaries of frame to avoid selecting gridpoint close to domain edge

306 lat_min_bound, lon_min_bound = (

307 lat.points[boundary_margin - 1],

308 lon.points[boundary_margin - 1],

309 )

310 lat_max_bound, lon_max_bound = (

311 lat.points[-boundary_margin],

312 lon.points[-boundary_margin],

313 )

314

315 # Check to see if selected point is outside the domain

316 if (lat_tr < lat_min) or (lat_tr > lat_max):

317 raise ValueError("Selected point is outside the domain.")

318 else:

319 if (lon_tr < lon_min) or (lon_tr > lon_max):

320 if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):

321 lon_tr += 360.0

322 elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):

323 lon_tr -= 360.0

324 else:

325 raise ValueError("Selected point is outside the domain.")

326

327 # Check to see if selected point is near the domain boundaries

328 if (

329 (lat_tr < lat_min_bound)

330 or (lat_tr > lat_max_bound)

331 or (lon_tr < lon_min_bound)

332 or (lon_tr > lon_max_bound)

333 ):

334 warnings.warn(

335 f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",

336 category=BoundaryWarning,

337 stacklevel=2,

338 )

339

340 regrid_method = getattr(iris.analysis, method, None)

341 if not callable(regrid_method):

342 raise NotImplementedError(

343 f"Does not currently support {method} regrid method"

344 )

345

346 cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())

347 regridded_cubes.append(cube_rgd)

348 else:

349 if (

350 np.abs((lat_tr - lat.points[0])) > 0.1

351 or np.abs((lon_tr - lon.points[0])) > 0.1

352 ):

353 raise ValueError(

354 "Selected point is too far from the specified coordinates. It should be within 0.1 degrees."

355 )

356 else:

357 print(

358 "*** lat/long diffs",

359 np.abs(lat_tr - lat_pt),

360 np.abs(lon_tr - lon_pt),

361 )

362 regridded_cubes.append(cube)

363

364 # Preserve returning a cube if only a cube has been supplied to regrid.

365 if len(regridded_cubes) == 1:

366 return regridded_cubes[0]

367 else:

368 return regridded_cubes

369

370

371def transform_lat_long_points(lon, lat, cube):

372 """Transform a selected point in longitude and latitude.

373

374 Transform the coordinates of a point from the real world

375 grid to the corresponding point on the rotated grid of a cube.

376

377 Parameters

378 ----------

379 cube: Cube

380 An iris cube of data defining the rotated grid to be used in

381 the longitude-latitude transformation.

382 lon: float

383 Selected value of longitude: this should be in the range -180 degrees to

384 180 degrees.

385 lat: float

386 Selected value of latitude: this should be in the range -90 degrees to

387 90 degrees.

388

389 Returns

390 -------

391 lon_rot, lat_rot: float

392 Coordinates of the selected point on the rotated grid specified within

393 the selected cube.

394

395 """

396 import cartopy.crs as ccrs

397

398 rot_pole = cube.coord_system().as_cartopy_crs()

399 true_grid = ccrs.Geodetic()

400 rot_coords = rot_pole.transform_point(lon, lat, true_grid)

401 lon_rot = rot_coords[0]

402 lat_rot = rot_coords[1]

403

404 return lon_rot, lat_rot

405

406

407def interpolate_to_point_cube(

408 fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs

409) -> iris.cube.Cube | iris.cube.CubeList:

410 """Interpolate a 2D field in cube or CubeList to a set of points.

411

412 Regrid cube(s) to set of sample points specified by point_cube.

413 Ensures only matching times between input and point_cube are included.

414

415 Parameters

416 ----------

417 fld: Cube | CubeList

418 An iris cube or CubeList containing a two-dimensional field(s).

419 point_cube: Cube

420 An iris cube specifying the coord point(s) to which the data

421 will be interpolated.

422

423 Returns

424 -------

425 fld_point_cube: Cube | CubeList

426 An iris cube or CubeList containing interpolated values at the

427 points specified by the point cube for matching times common to

428 both fld and point_cube.

429 """

430 # Empty CubeList To store regridded cubes.

431 regridded_cubes = iris.cube.CubeList()

432

433 # Iterate over all cubes and regrid.

434 for cube in iter_maybe(fld):

435 # Ensure matching times in fld cube and point_cube

436 base_time_coord = point_cube.coord("time")

437 other_time_coord = cube.coord("time")

438 base_times = base_time_coord.units.num2date(base_time_coord.points)

439 other_times = other_time_coord.units.num2date(other_time_coord.points)

440 shared_times = set.intersection(set(base_times), set(other_times))

441 logging.debug("Shared times: %s", shared_times)

442 time_constraint = iris.Constraint(

443 coord_values={

444 "time": lambda cell, shared_times=shared_times: (

445 cell.point in shared_times

446 )

447 }

448 )

449

450 # Extract points matching the shared times.

451 cube = cube.extract(time_constraint)

452 point_cube = point_cube.extract(time_constraint)

453 if cube is None or point_cube is None: 453 ↛ 454line 453 didn't jump to line 454 because the condition on line 453 was never true

454 raise ValueError("No common time points found!")

455

456 # Generate array of point cube lat and lon points.

457 lat_name, lon_name = get_cube_yxcoordname(point_cube)

458 lats = point_cube.coord(lat_name).points

459 lons = point_cube.coord(lon_name).points

460

461 y_coord, x_coord = get_cube_yxcoordname(cube)

462 if isinstance( 462 ↛ 473line 462 didn't jump to line 473 because the condition on line 462 was always true

463 cube.coord(x_coord).coord_system, iris.coord_systems.RotatedGeogCS

464 ):

465 lons_rp, lats_rp = rotate_pole(

466 lons,

467 lats,

468 pole_lon=cube.coord(x_coord).coord_system.grid_north_pole_longitude,

469 pole_lat=cube.coord(x_coord).coord_system.grid_north_pole_latitude,

470 )

471 sample_points = [("grid_latitude", lats_rp), ("grid_longitude", lons_rp)]

472 else:

473 sample_points = [

474 ("latitude", np.array(lats)),

475 ("longitude", np.array(lons)),

476 ]

477

478 # Interpolate fld cube to required sample points

479 fld_point_cube = cube.interpolate(

480 sample_points, iris.analysis.Linear(extrapolation_mode="mask")

481 )

482

483 # Retain only diagonal elements of 2D interpolated cube to vector points

484 od_index = point_cube.coord_dims("station")[0]

485 fv_cube = iris.cube.Cube(

486 fld_point_cube.data.diagonal(axis1=od_index, axis2=od_index + 1),

487 standard_name=cube.standard_name,

488 long_name=cube.long_name,

489 units=cube.units,

490 )

491 # Copy all non-lat/lon coordinates and cube attributes

492 if "time" in [coord.name() for coord in fld_point_cube.coords(dim_coords=True)]: 492 ↛ 496line 492 didn't jump to line 496 because the condition on line 492 was always true

493 fv_cube.add_dim_coord(

494 point_cube.coord("time"), point_cube.coord_dims("time")[0]

495 )

496 for coord in cube.coords():

497 if coord.name() not in [

498 "time",

499 "latitude",

500 "longitude",

501 "grid_latitude",

502 "grid_longitude",

503 ] and coord.name() not in [coord.name() for coord in fv_cube.coords()]:

504 fv_cube.add_aux_coord(coord.copy(), cube.coord_dims(coord))

505 for coord in point_cube.coords():

506 if coord.name() not in [

507 "time",

508 "forecast_period",

509 "forecast_reference_time",

510 "realization",

511 "station",

512 ] and coord.name() not in [coord.name() for coord in fv_cube.coords()]:

513 fv_cube.add_aux_coord(coord.copy(), point_cube.coord_dims(coord))

514 fv_cube.add_dim_coord(point_cube.coord("station"), od_index)

515 fv_cube.attributes = cube.attributes.copy()

516 fv_cube.cell_methods = cube.cell_methods

517 fv_cube.units = cube.units

518 regridded_cubes.append(fv_cube)

519

520 # Preserve returning a cube if only a cube has been supplied to regrid.

521 if len(regridded_cubes) == 1:

522 return regridded_cubes[0]

523 else:

524 return regridded_cubes

525

526

527def vertical_interpolation(

528 cubes: iris.cube.Cube | iris.cube.CubeList,

529 coordinate: str,

530 target: iris.cube.Cube | iris.cube.CubeList,

531) -> iris.cube.Cube | iris.cube.CubeList:

532 """Vertical interpolation of a cube to match that off a different cube.

533

534 Acts as a wrapper around the `cube.interpolate` functionality and uses

535 linear interpolation as the method.

536

537 Parameters

538 ----------

539 cubes: iris.cube.Cube | iris.cube.CubeList

540 An iris cube or cubelist of data defining field that should be

541 vertically interpolated.

542 coordinate: str

543 The coordinate the interpolation occurs over.

544 target: iris.cube.Cube | iris.cube.CubeList

545 The target cube or cubelist that provides the vertical coordinate

546 information. It will use `cube.coord(coordinate).points` to provide

547 the vertical target. The number of target cubes should match the number

548 of cubes used as input.

549

550 Returns

551 -------

552 interpolated_cubes: iris.cube.Cube | iris.cube.CubeList

553 Coordinates of the selected point on the rotated grid specified within

554 the selected cube.

555 """

556 interpolated_cubes = iris.cube.CubeList([])

557 for cube, cube_t in zip(iter_maybe(cubes), iter_maybe(target), strict=True):

558 target_levels = cube_t.coord(coordinate).points

559 new_cube = cube.interpolate(

560 [(coordinate, target_levels)], iris.analysis.Linear()

561 )

562 interpolated_cubes.append(new_cube)

563 if len(interpolated_cubes) == 1:

564 return interpolated_cubes[0]

565 else:

566 return interpolated_cubes