Coverage for src/CSET/operators/plot.py: 81%

1# © Crown copyright, Met Office (2022-2025) 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 produce various kinds of plots."""

17import fcntl

18import functools

19import importlib.resources

20import itertools

21import json

22import logging

23import math

24import os

25from typing import Literal

27import cartopy.crs as ccrs

28import cartopy.feature as cfeature

29import iris

30import iris.coords

31import iris.cube

32import iris.exceptions

33import iris.plot as iplt

34import matplotlib as mpl

35import matplotlib.colors as mcolors

36import matplotlib.pyplot as plt

37import numpy as np

38from iris.cube import Cube

39from markdown_it import MarkdownIt

41from CSET._common import (

42 combine_dicts,

43 filename_slugify,

44 get_recipe_metadata,

45 iter_maybe,

46 render_file,

47 slugify,

48)

49from CSET.operators._utils import (

50 fully_equalise_attributes,

51 get_cube_yxcoordname,

52 is_transect,

53 slice_over_maybe,

54)

55from CSET.operators.collapse import collapse

56from CSET.operators.misc import _extract_common_time_points

57from CSET.operators.regrid import regrid_onto_cube

59# Suppress matplotlib debug output

60# logging.getLogger('matplotlib.ticker').setLevel(logging.WARNING)

62# Use a non-interactive plotting backend.

63mpl.use("agg")

65DEFAULT_DISCRETE_COLORS = mpl.colormaps["tab10"].colors + mpl.colormaps["Accent"].colors

67############################

68# Private helper functions #

69############################

72def _append_to_plot_index(plot_index: list) -> list:

73 """Add plots into the plot index, returning the complete plot index."""

74 with open("meta.json", "r+t", encoding="UTF-8") as fp:

75 fcntl.flock(fp, fcntl.LOCK_EX)

76 fp.seek(0)

77 meta = json.load(fp)

78 complete_plot_index = meta.get("plots", [])

79 complete_plot_index = complete_plot_index + plot_index

80 meta["plots"] = complete_plot_index

81 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(

82 os.getenv("DO_CASE_AGGREGATION")

83 ):

84 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")

85 fp.seek(0)

86 fp.truncate()

87 json.dump(meta, fp, indent=2)

88 return complete_plot_index

91def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:

92 """Ensure a single cube is given.

94 If a CubeList of length one is given that the contained cube is returned,

95 otherwise an error is raised.

97 Parameters

98 ----------

99 cube: Cube | CubeList

100 The cube to check.

101

102 Returns

103 -------

104 cube: Cube

105 The checked cube.

106

107 Raises

108 ------

109 TypeError

110 If the input cube is not a Cube or CubeList of a single Cube.

111 """

112 if isinstance(cube, iris.cube.Cube):

113 return cube

114 if isinstance(cube, iris.cube.CubeList):

115 if len(cube) == 1:

116 return cube[0]

117 raise TypeError("Must have a single cube", cube)

118

119

120def _make_plot_html_page(plots: list):

121 """Create a HTML page to display a plot image."""

122 # Debug check that plots actually contains some strings.

123 assert isinstance(plots[0], str)

124

125 # Load HTML template file.

126 operator_files = importlib.resources.files()

127 template_file = operator_files.joinpath("_plot_page_template.html")

128

129 # Get some metadata.

130 meta = get_recipe_metadata()

131 title = meta.get("title", "Untitled")

132 description = MarkdownIt().render(meta.get("description", "*No description.*"))

133

134 # Prepare template variables.

135 variables = {

136 "title": title,

137 "description": description,

138 "initial_plot": plots[0],

139 "plots": plots,

140 "title_slug": slugify(title),

141 }

142

143 # Render template.

144 html = render_file(template_file, **variables)

145

146 # Save completed HTML.

147 with open("index.html", "wt", encoding="UTF-8") as fp:

148 fp.write(html)

149

150

151@functools.cache

152def _load_colorbar_map(user_colorbar_file: str = None) -> dict:

153 """Load the colorbar definitions from a file.

154

155 This is a separate function to make it cacheable.

156 """

157 colorbar_file = importlib.resources.files().joinpath("_colorbar_definition.json")

158 with open(colorbar_file, "rt", encoding="UTF-8") as fp:

159 colorbar = json.load(fp)

160

161 logging.debug("User colour bar file: %s", user_colorbar_file)

162 override_colorbar = {}

163 if user_colorbar_file:

164 try:

165 with open(user_colorbar_file, "rt", encoding="UTF-8") as fp:

166 override_colorbar = json.load(fp)

167 except FileNotFoundError:

168 logging.warning("Colorbar file does not exist. Using default values.")

169

170 # Overwrite values with the user supplied colorbar definition.

171 colorbar = combine_dicts(colorbar, override_colorbar)

172 return colorbar

173

174

175def _get_model_colors_map(cubes: iris.cube.CubeList | iris.cube.Cube) -> dict:

176 """Get an appropriate colors for model lines in line plots.

177

178 For each model in the list of cubes colors either from user provided

179 color definition file (so-called style file) or from default colors are mapped

180 to model_name attribute.

181

182 Parameters

183 ----------

184 cubes: CubeList or Cube

185 Cubes with model_name attribute

186

187 Returns

188 -------

189 model_colors_map:

190 Dictionary mapping model_name attribute to colors

191 """

192 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)

193 colorbar = _load_colorbar_map(user_colorbar_file)

194 model_names = sorted(

195 filter(

196 lambda x: x is not None,

197 (cube.attributes.get("model_name", None) for cube in iter_maybe(cubes)),

198 )

199 )

200 if not model_names:

201 return {}

202 use_user_colors = all(mname in colorbar.keys() for mname in model_names)

203 if use_user_colors: 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true

204 return {mname: colorbar[mname] for mname in model_names}

205

206 color_list = itertools.cycle(DEFAULT_DISCRETE_COLORS)

207 return {mname: color for mname, color in zip(model_names, color_list, strict=False)}

208

209

210def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):

211 """Get an appropriate colorbar for the given cube.

212

213 For the given variable the appropriate colorbar is looked up from a

214 combination of the built-in CSET colorbar definitions, and any user supplied

215 definitions. As well as varying on variables, these definitions may also

216 exist for specific pressure levels to account for variables with

217 significantly different ranges at different heights. The colorbars also exist

218 for masks and mask differences for considering variable presence diagnostics.

219 Specific variable ranges can be separately set in user-supplied definition

220 for x- or y-axis limits, or indicate where automated range preferred.

221

222 Parameters

223 ----------

224 cube: Cube

225 Cube of variable for which the colorbar information is desired.

226 axis: "x", "y", optional

227 Select the levels for just this axis of a line plot. The min and max

228 can be set by xmin/xmax or ymin/ymax respectively. For variables where

229 setting a universal range is not desirable (e.g. temperature), users

230 can set ymin/ymax values to "auto" in the colorbar definitions file.

231 Where no additional xmin/xmax or ymin/ymax values are provided, the

232 axis bounds default to use the vmin/vmax values provided.

233

234 Returns

235 -------

236 cmap:

237 Matplotlib colormap.

238 levels:

239 List of levels to use for plotting. For continuous plots the min and max

240 should be taken as the range.

241 norm:

242 BoundaryNorm information.

243 """

244 # Grab the colorbar file from the recipe global metadata.

245 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)

246 colorbar = _load_colorbar_map(user_colorbar_file)

247 cmap = None

248

249 try:

250 # We assume that pressure is a scalar coordinate here.

251 pressure_level_raw = cube.coord("pressure").points[0]

252 # Ensure pressure_level is a string, as it is used as a JSON key.

253 pressure_level = str(int(pressure_level_raw))

254 except iris.exceptions.CoordinateNotFoundError:

255 pressure_level = None

256

257 # First try long name, then standard name, then var name. This order is used

258 # as long name is the one we correct between models, so it most likely to be

259 # consistent.

260 varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name]))

261 for varname in varnames:

262 # Get the colormap for this variable.

263 try:

264 var_colorbar = colorbar[varname]

265 cmap = plt.get_cmap(colorbar[varname]["cmap"], 51)

266 varname_key = varname

267 break

268 except KeyError:

269 logging.debug("Cube name %s has no colorbar definition.", varname)

270

271 # Get colormap if it is a mask.

272 if any("mask_for_" in name for name in varnames):

273 cmap, levels, norm = _custom_colormap_mask(cube, axis=axis)

274 return cmap, levels, norm

275 # If winds on Beaufort Scale use custom colorbar and levels

276 if any("Beaufort_Scale" in name for name in varnames):

277 cmap, levels, norm = _custom_beaufort_scale(cube, axis=axis)

278 return cmap, levels, norm

279 # If probability is plotted use custom colorbar and levels

280 if any("probability_of_" in name for name in varnames):

281 cmap, levels, norm = _custom_colormap_probability(cube, axis=axis)

282 return cmap, levels, norm

283 # If aviation colour state use custom colorbar and levels

284 if any("aviation_colour_state" in name for name in varnames):

285 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)

286 return cmap, levels, norm

287

288 # If no valid colormap has been defined, use defaults and return.

289 if not cmap:

290 logging.warning("No colorbar definition exists for %s.", cube.name())

291 cmap, levels, norm = mpl.colormaps["viridis"], None, None

292 return cmap, levels, norm

293

294 # Test if pressure-level specific settings are provided for cube.

295 if pressure_level:

296 try:

297 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]

298 except KeyError:

299 logging.debug(

300 "%s has no colorbar definition for pressure level %s.",

301 varname,

302 pressure_level,

303 )

304

305 # Check for availability of x-axis or y-axis user-specific overrides

306 # for setting level bounds for line plot types and return just levels.

307 # Line plots do not need a colormap, and just use the data range.

308 if axis:

309 if axis == "x":

310 try:

311 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]

312 except KeyError:

313 vmin, vmax = var_colorbar["min"], var_colorbar["max"]

314 if axis == "y":

315 try:

316 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]

317 except KeyError:

318 vmin, vmax = var_colorbar["min"], var_colorbar["max"]

319 # Check if user-specified auto-scaling for this variable

320 if vmin == "auto" or vmax == "auto":

321 levels = None

322 else:

323 levels = [vmin, vmax]

324 return None, levels, None

325 # Get and use the colorbar levels for this variable if spatial or histogram.

326 else:

327 try:

328 levels = var_colorbar["levels"]

329 # Use discrete bins when levels are specified, rather

330 # than a smooth range.

331 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)

332 logging.debug("Using levels for %s colorbar.", varname)

333 logging.info("Using levels: %s", levels)

334 except KeyError:

335 # Get the range for this variable.

336 vmin, vmax = var_colorbar["min"], var_colorbar["max"]

337 logging.debug("Using min and max for %s colorbar.", varname)

338 # Calculate levels from range.

339 levels = np.linspace(vmin, vmax, 101)

340 norm = None

341

342 # Overwrite cmap, levels and norm for specific variables that

343 # require custom colorbar_map as these can not be defined in the

344 # JSON file.

345 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)

346 cmap, levels, norm = _custom_colourmap_visibility_in_air(

347 cube, cmap, levels, norm

348 )

349 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)

350 return cmap, levels, norm

351

352

353def _setup_spatial_map(

354 cube: iris.cube.Cube,

355 figure,

356 cmap,

357 grid_size: int | None = None,

358 subplot: int | None = None,

359):

360 """Define map projections, extent and add coastlines and borderlines for spatial plots.

361

362 For spatial map plots, a relevant map projection for rotated or non-rotated inputs

363 is specified, and map extent defined based on the input data.

364

365 Parameters

366 ----------

367 cube: Cube

368 2 dimensional (lat and lon) Cube of the data to plot.

369 figure:

370 Matplotlib Figure object holding all plot elements.

371 cmap:

372 Matplotlib colormap.

373 grid_size: int, optional

374 Size of grid for subplots if multiple spatial subplots in figure.

375 subplot: int, optional

376 Subplot index if multiple spatial subplots in figure.

377

378 Returns

379 -------

380 axes:

381 Matplotlib GeoAxes definition.

382 """

383 # Identify min/max plot bounds.

384 try:

385 lat_axis, lon_axis = get_cube_yxcoordname(cube)

386 x1 = np.min(cube.coord(lon_axis).points)

387 x2 = np.max(cube.coord(lon_axis).points)

388 y1 = np.min(cube.coord(lat_axis).points)

389 y2 = np.max(cube.coord(lat_axis).points)

390

391 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.

392 if np.abs(x2 - x1) > 180.0:

393 x1 = x1 - 180.0

394 x2 = x2 - 180.0

395 logging.debug("Adjusting plot bounds to fit global extent.")

396

397 # Consider map projection orientation.

398 # Adapting orientation enables plotting across international dateline.

399 # Users can adapt the default central_longitude if alternative projections views.

400 if x2 > 180.0 or x1 < -180.0:

401 central_longitude = 180.0

402 else:

403 central_longitude = 0.0

404

405 # Define spatial map projection.

406 coord_system = cube.coord(lat_axis).coord_system

407 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):

408 # Define rotated pole map projection for rotated pole inputs.

409 projection = ccrs.RotatedPole(

410 pole_longitude=coord_system.grid_north_pole_longitude,

411 pole_latitude=coord_system.grid_north_pole_latitude,

412 central_rotated_longitude=central_longitude,

413 )

414 crs = projection

415 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 415 ↛ 417line 415 didn't jump to line 417 because the condition on line 415 was never true

416 # Define Transverse Mercator projection for TM inputs.

417 projection = ccrs.TransverseMercator(

418 central_longitude=coord_system.longitude_of_central_meridian,

419 central_latitude=coord_system.latitude_of_projection_origin,

420 false_easting=coord_system.false_easting,

421 false_northing=coord_system.false_northing,

422 scale_factor=coord_system.scale_factor_at_central_meridian,

423 )

424 crs = projection

425 else:

426 # Define regular map projection for non-rotated pole inputs.

427 # Alternatives might include e.g. for global model outputs:

428 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)

429 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.

430 projection = ccrs.PlateCarree(central_longitude=central_longitude)

431 crs = ccrs.PlateCarree()

432

433 # Define axes for plot (or subplot) with required map projection.

434 if subplot is not None:

435 axes = figure.add_subplot(

436 grid_size, grid_size, subplot, projection=projection

437 )

438 else:

439 axes = figure.add_subplot(projection=projection)

440

441 # Add coastlines and borderlines if cube contains x and y map coordinates.

442 if cmap.name in ["viridis", "Greys"]:

443 coastcol = "magenta"

444 else:

445 coastcol = "black"

446 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)

447 axes.coastlines(resolution="10m", color=coastcol)

448 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)

449

450 # If is lat/lon spatial map, fix extent to keep plot tight.

451 # Specifying crs within set_extent helps ensure only data region is shown.

452 if isinstance(coord_system, iris.coord_systems.GeogCS):

453 axes.set_extent([x1, x2, y1, y2], crs=crs)

454

455 except ValueError:

456 # Skip if not both x and y map coordinates.

457 axes = figure.gca()

458 pass

459

460 return axes

461

462

463def _get_plot_resolution() -> int:

464 """Get resolution of rasterised plots in pixels per inch."""

465 return get_recipe_metadata().get("plot_resolution", 100)

466

467

468def _set_title_and_filename(

469 seq_coord: iris.coords.Coord,

470 nplot: int,

471 recipe_title: str,

472 filename: str,

473):

474 """Set plot title and filename based on cube coordinate.

475

476 Parameters

477 ----------

478 sequence_coordinate: iris.coords.Coord

479 Coordinate about which to make a plot sequence.

480 nplot: int

481 Number of output plots to generate - controls title/naming.

482 recipe_title: str

483 Default plot title, potentially to update.

484 filename: str

485 Input plot filename, potentially to update.

486

487 Returns

488 -------

489 plot_title: str

490 Output formatted plot title string, based on plotted data.

491 plot_filename: str

492 Output formatted plot filename string.

493 """

494 ndim = seq_coord.ndim

495 npoints = np.size(seq_coord.points)

496 sequence_title = ""

497 sequence_fname = ""

498

499 # Account for case with multi-dimension sequence input (e.g. aggregation)

500 if ndim > 1:

501 sequence_title = f"\n [{ndim} cases]"

502 sequence_fname = f"_{ndim}cases"

503

504 else:

505 if npoints == 1:

506 if nplot > 1:

507 # Set default labels for sequence inputs

508 sequence_value = seq_coord.units.title(seq_coord.points[0])

509 sequence_title = f"\n [{sequence_value}]"

510 sequence_fname = f"_{filename_slugify(sequence_value)}"

511 elif seq_coord.has_bounds():

512 ncase = np.size(seq_coord.bounds)

513 sequence_title = f"\n [{ncase} cases]"

514 sequence_fname = f"_{ncase}cases"

515 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs

516 # Take title endpoints from coord points where series input (e.g. timeseries)

517 if npoints > 1:

518 if not seq_coord.has_bounds():

519 startstring = seq_coord.units.title(seq_coord.points[0])

520 endstring = seq_coord.units.title(seq_coord.points[-1])

521 else:

522 # Take title endpoint from coord bounds where single input (e.g. map)

523 startstring = seq_coord.units.title(seq_coord.bounds.flatten()[0])

524 endstring = seq_coord.units.title(seq_coord.bounds.flatten()[-1])

525 sequence_title = f"\n [{startstring} to {endstring}]"

526 sequence_fname = (

527 f"_{filename_slugify(startstring)}_{filename_slugify(endstring)}"

528 )

529

530 # Set plot title and filename

531 plot_title = f"{recipe_title}{sequence_title}"

532

533 # Set plot filename, defaulting to user input if provided.

534 if filename is None:

535 filename = slugify(recipe_title)

536 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"

537 else:

538 if nplot > 1:

539 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"

540 else:

541 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"

542

543 return plot_title, plot_filename

544

545

546def select_series_coord(cube, series_coordinate):

547 """Determine the grid coordinates to use to calculate grid spacing."""

548 try:

549 # Try the requested coordinate first

550 return cube.coord(series_coordinate)

551

552 except iris.exceptions.CoordinateNotFoundError:

553 # Fallback logic

554 if series_coordinate == "frequency":

555 fallbacks = ("physical_wavenumber", "wavelength")

556

557 elif series_coordinate == "physical_wavenumber":

558 fallbacks = ("frequency", "wavelength")

559

560 elif series_coordinate == "wavelength": 560 ↛ 565line 560 didn't jump to line 565 because the condition on line 560 was always true

561 fallbacks = ("frequency", "physical_wavenumber")

562

563 else:

564 # Unknown coordinate type -> re-raise the original exception

565 raise

566

567 # Try the fallbacks in order

568 for fb in fallbacks: 568 ↛ 576line 568 didn't jump to line 576 because the loop on line 568 didn't complete

569 try:

570 return cube.coord(fb)

571 except iris.exceptions.CoordinateNotFoundError:

572 continue

573

574 # If we get here, none of the fallback options were found

575

576 raise iris.exceptions.CoordinateNotFoundError(

577 f"No valid coordinate found for '{series_coordinate}' "

578 f"or fallback options {fallbacks}"

579 ) from None

580

581

582def _plot_and_save_spatial_plot(

583 cube: iris.cube.Cube,

584 filename: str,

585 title: str,

586 method: Literal["contourf", "pcolormesh"],

587 overlay_cube: iris.cube.Cube | None = None,

588 contour_cube: iris.cube.Cube | None = None,

589 **kwargs,

590):

591 """Plot and save a spatial plot.

592

593 Parameters

594 ----------

595 cube: Cube

596 2 dimensional (lat and lon) Cube of the data to plot.

597 filename: str

598 Filename of the plot to write.

599 title: str

600 Plot title.

601 method: "contourf" | "pcolormesh"

602 The plotting method to use.

603 overlay_cube: Cube, optional

604 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube

605 contour_cube: Cube, optional

606 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube

607 """

608 # Setup plot details, size, resolution, etc.

609 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

610

611 # Specify the color bar

612 cmap, levels, norm = _colorbar_map_levels(cube)

613

614 # If overplotting, set required colorbars

615 if overlay_cube:

616 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)

617 if contour_cube:

618 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)

619

620 # Setup plot map projection, extent and coastlines and borderlines.

621 axes = _setup_spatial_map(cube, fig, cmap)

622

623 # Plot the field.

624 if method == "contourf":

625 # Filled contour plot of the field.

626 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)

627 elif method == "pcolormesh":

628 try:

629 vmin = min(levels)

630 vmax = max(levels)

631 except TypeError:

632 vmin, vmax = None, None

633 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax

634 # if levels are defined.

635 if norm is not None:

636 vmin = None

637 vmax = None

638 logging.debug("Plotting using defined levels.")

639 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)

640 else:

641 raise ValueError(f"Unknown plotting method: {method}")

642

643 # Overplot overlay field, if required

644 if overlay_cube:

645 try:

646 over_vmin = min(over_levels)

647 over_vmax = max(over_levels)

648 except TypeError:

649 over_vmin, over_vmax = None, None

650 if over_norm is not None: 650 ↛ 651line 650 didn't jump to line 651 because the condition on line 650 was never true

651 over_vmin = None

652 over_vmax = None

653 overlay = iplt.pcolormesh(

654 overlay_cube,

655 cmap=over_cmap,

656 norm=over_norm,

657 alpha=0.8,

658 vmin=over_vmin,

659 vmax=over_vmax,

660 )

661 # Overplot contour field, if required, with contour labelling.

662 if contour_cube:

663 contour = iplt.contour(

664 contour_cube,

665 colors="darkgray",

666 levels=cntr_levels,

667 norm=cntr_norm,

668 alpha=0.5,

669 linestyles="--",

670 linewidths=1,

671 )

672 plt.clabel(contour)

673

674 # Check to see if transect, and if so, adjust y axis.

675 if is_transect(cube):

676 if "pressure" in [coord.name() for coord in cube.coords()]:

677 axes.invert_yaxis()

678 axes.set_yscale("log")

679 axes.set_ylim(1100, 100)

680 # If both model_level_number and level_height exists, iplt can construct

681 # plot as a function of height above orography (NOT sea level).

682 elif {"model_level_number", "level_height"}.issubset( 682 ↛ 687line 682 didn't jump to line 687 because the condition on line 682 was always true

683 {coord.name() for coord in cube.coords()}

684 ):

685 axes.set_yscale("log")

686

687 axes.set_title(

688 f"{title}\n"

689 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"

690 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"

691 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"

692 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",

693 fontsize=16,

694 )

695

696 else:

697 # Add title.

698 axes.set_title(title, fontsize=16)

699

700 # Add watermark with min/max/mean. Currently not user togglable.

701 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.

702 axes.annotate(

703 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",

704 xy=(1, -0.05),

705 xycoords="axes fraction",

706 xytext=(-5, 5),

707 textcoords="offset points",

708 ha="right",

709 va="bottom",

710 size=11,

711 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),

712 )

713

714 # Add secondary colour bar for overlay_cube field if required.

715 if overlay_cube:

716 cbarB = fig.colorbar(

717 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7

718 )

719 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)

720 # add ticks and tick_labels for every levels if less than 20 levels exist

721 if over_levels is not None and len(over_levels) < 20: 721 ↛ 722line 721 didn't jump to line 722 because the condition on line 721 was never true

722 cbarB.set_ticks(over_levels)

723 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])

724 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():

725 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])

726 logging.debug("Set secondary colorbar ticks and labels.")

727

728 # Add main colour bar.

729 cbar = fig.colorbar(

730 plot, orientation="horizontal", location="bottom", pad=0.042, shrink=0.7

731 )

732 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)

733 # add ticks and tick_labels for every levels if less than 20 levels exist

734 if levels is not None and len(levels) < 20:

735 cbar.set_ticks(levels)

736 cbar.set_ticklabels([f"{level:.2f}" for level in levels])

737 if "rainfall" or "snowfall" or "visibility" in cube.name(): 737 ↛ 739line 737 didn't jump to line 739 because the condition on line 737 was always true

738 cbar.set_ticklabels([f"{level:.3g}" for level in levels])

739 logging.debug("Set colorbar ticks and labels.")

740

741 # Save plot.

742 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

743 logging.info("Saved spatial plot to %s", filename)

744 plt.close(fig)

745

746

747def _plot_and_save_postage_stamp_spatial_plot(

748 cube: iris.cube.Cube,

749 filename: str,

750 stamp_coordinate: str,

751 title: str,

752 method: Literal["contourf", "pcolormesh"],

753 overlay_cube: iris.cube.Cube | None = None,

754 contour_cube: iris.cube.Cube | None = None,

755 **kwargs,

756):

757 """Plot postage stamp spatial plots from an ensemble.

758

759 Parameters

760 ----------

761 cube: Cube

762 Iris cube of data to be plotted. It must have the stamp coordinate.

763 filename: str

764 Filename of the plot to write.

765 stamp_coordinate: str

766 Coordinate that becomes different plots.

767 method: "contourf" | "pcolormesh"

768 The plotting method to use.

769 overlay_cube: Cube, optional

770 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube

771 contour_cube: Cube, optional

772 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube

773

774 Raises

775 ------

776 ValueError

777 If the cube doesn't have the right dimensions.

778 """

779 # Use the smallest square grid that will fit the members.

780 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))

781

782 fig = plt.figure(figsize=(10, 10))

783

784 # Specify the color bar

785 cmap, levels, norm = _colorbar_map_levels(cube)

786 # If overplotting, set required colorbars

787 if overlay_cube: 787 ↛ 788line 787 didn't jump to line 788 because the condition on line 787 was never true

788 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)

789 if contour_cube: 789 ↛ 790line 789 didn't jump to line 790 because the condition on line 789 was never true

790 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)

791

792 # Make a subplot for each member.

793 for member, subplot in zip(

794 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False

795 ):

796 # Setup subplot map projection, extent and coastlines and borderlines.

797 axes = _setup_spatial_map(

798 member, fig, cmap, grid_size=grid_size, subplot=subplot

799 )

800 if method == "contourf":

801 # Filled contour plot of the field.

802 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)

803 elif method == "pcolormesh":

804 if levels is not None:

805 vmin = min(levels)

806 vmax = max(levels)

807 else:

808 raise TypeError("Unknown vmin and vmax range.")

809 vmin, vmax = None, None

810 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax

811 # if levels are defined.

812 if norm is not None: 812 ↛ 813line 812 didn't jump to line 813 because the condition on line 812 was never true

813 vmin = None

814 vmax = None

815 # pcolormesh plot of the field.

816 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)

817 else:

818 raise ValueError(f"Unknown plotting method: {method}")

819

820 # Overplot overlay field, if required

821 if overlay_cube: 821 ↛ 822line 821 didn't jump to line 822 because the condition on line 821 was never true

822 try:

823 over_vmin = min(over_levels)

824 over_vmax = max(over_levels)

825 except TypeError:

826 over_vmin, over_vmax = None, None

827 if over_norm is not None:

828 over_vmin = None

829 over_vmax = None

830 iplt.pcolormesh(

831 overlay_cube[member.coord(stamp_coordinate).points[0]],

832 cmap=over_cmap,

833 norm=over_norm,

834 alpha=0.6,

835 vmin=over_vmin,

836 vmax=over_vmax,

837 )

838 # Overplot contour field, if required

839 if contour_cube: 839 ↛ 840line 839 didn't jump to line 840 because the condition on line 839 was never true

840 iplt.contour(

841 contour_cube[member.coord(stamp_coordinate).points[0]],

842 colors="darkgray",

843 levels=cntr_levels,

844 norm=cntr_norm,

845 alpha=0.6,

846 linestyles="--",

847 linewidths=1,

848 )

849 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")

850 axes.set_axis_off()

851

852 # Put the shared colorbar in its own axes.

853 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])

854 colorbar = fig.colorbar(

855 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7

856 )

857 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)

858

859 # Overall figure title.

860 fig.suptitle(title, fontsize=16)

861

862 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

863 logging.info("Saved contour postage stamp plot to %s", filename)

864 plt.close(fig)

865

866

867def _plot_and_save_line_series(

868 cubes: iris.cube.CubeList,

869 coords: list[iris.coords.Coord],

870 ensemble_coord: str,

871 filename: str,

872 title: str,

873 **kwargs,

874):

875 """Plot and save a 1D line series.

876

877 Parameters

878 ----------

879 cubes: Cube or CubeList

880 Cube or CubeList containing the cubes to plot on the y-axis.

881 coords: list[Coord]

882 Coordinates to plot on the x-axis, one per cube.

883 ensemble_coord: str

884 Ensemble coordinate in the cube.

885 filename: str

886 Filename of the plot to write.

887 title: str

888 Plot title.

889 """

890 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

891

892 model_colors_map = _get_model_colors_map(cubes)

893

894 # Store min/max ranges.

895 y_levels = []

896

897 # Check match-up across sequence coords gives consistent sizes

898 _validate_cubes_coords(cubes, coords)

899

900 for cube, coord in zip(cubes, coords, strict=True):

901 label = None

902 color = "black"

903 if model_colors_map:

904 label = cube.attributes.get("model_name")

905 color = model_colors_map.get(label)

906 for cube_slice in cube.slices_over(ensemble_coord):

907 # Label with (control) if part of an ensemble or not otherwise.

908 if cube_slice.coord(ensemble_coord).points == [0]:

909 iplt.plot(

910 coord,

911 cube_slice,

912 color=color,

913 marker="o",

914 ls="-",

915 lw=3,

916 label=f"{label} (control)"

917 if len(cube.coord(ensemble_coord).points) > 1

918 else label,

919 )

920 # Label with (perturbed) if part of an ensemble and not the control.

921 else:

922 iplt.plot(

923 coord,

924 cube_slice,

925 color=color,

926 ls="-",

927 lw=1.5,

928 alpha=0.75,

929 label=f"{label} (member)",

930 )

931

932 # Calculate the global min/max if multiple cubes are given.

933 _, levels, _ = _colorbar_map_levels(cube, axis="y")

934 if levels is not None: 934 ↛ 935line 934 didn't jump to line 935 because the condition on line 934 was never true

935 y_levels.append(min(levels))

936 y_levels.append(max(levels))

937

938 # Get the current axes.

939 ax = plt.gca()

940

941 # Add some labels and tweak the style.

942 # check if cubes[0] works for single cube if not CubeList

943 if coords[0].name() == "time":

944 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)

945 else:

946 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)

947 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)

948 ax.set_title(title, fontsize=16)

949

950 ax.ticklabel_format(axis="y", useOffset=False)

951 ax.tick_params(axis="x", labelrotation=15)

952 ax.tick_params(axis="both", labelsize=12)

953

954 # Set y limits to global min and max, autoscale if colorbar doesn't exist.

955 if y_levels: 955 ↛ 956line 955 didn't jump to line 956 because the condition on line 955 was never true

956 ax.set_ylim(min(y_levels), max(y_levels))

957 # Add zero line.

958 if min(y_levels) < 0.0 and max(y_levels) > 0.0:

959 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)

960 logging.debug(

961 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)

962 )

963 else:

964 ax.autoscale()

965

966 # Add gridlines

967 ax.grid(linestyle="--", color="grey", linewidth=1)

968 # Ientify unique labels for legend

969 handles = list(

970 {

971 label: handle

972 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)

973 }.values()

974 )

975 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)

976

977 # Save plot.

978 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

979 logging.info("Saved line plot to %s", filename)

980 plt.close(fig)

981

982

983def _plot_and_save_line_power_spectrum_series(

984 cubes: iris.cube.CubeList,

985 coords: list[iris.coords.Coord],

986 ensemble_coord: str,

987 filename: str,

988 title: str,

989 series_coordinate: str = None,

990 **kwargs,

991):

992 """Plot and save a 1D line series.

993

994 Parameters

995 ----------

996 cubes: Cube or CubeList

997 Cube or CubeList containing the cubes to plot on the y-axis.

998 coords: list[Coord]

999 Coordinates to plot on the x-axis, one per cube.

1000 ensemble_coord: str

1001 Ensemble coordinate in the cube.

1002 filename: str

1003 Filename of the plot to write.

1004 title: str

1005 Plot title.

1006 series_coordinate: str, optional

1007 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.

1008 """

1009 # xn = coords[0].name() # x-axis (e.g. frequency)

1010

1011 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1012 model_colors_map = _get_model_colors_map(cubes)

1013 ax = plt.gca()

1014

1015 # Store min/max ranges.

1016 y_levels = []

1017

1018 line_marker = None

1019 line_width = 1

1020

1021 # for cube, coord in zip(cubes, coords, strict=True):

1022 for cube in iter_maybe(cubes):

1023 # next 2 lines replace chunk of code.

1024 xcoord = select_series_coord(cube, series_coordinate)

1025 xname = xcoord.points

1026

1027 yfield = cube.data # power spectrum

1028 label = None

1029 color = "black"

1030 if model_colors_map: 1030 ↛ 1033line 1030 didn't jump to line 1033 because the condition on line 1030 was always true

1031 label = cube.attributes.get("model_name")

1032 color = model_colors_map.get(label)

1033 for cube_slice in cube.slices_over(ensemble_coord):

1034 # Label with (control) if part of an ensemble or not otherwise.

1035 if cube_slice.coord(ensemble_coord).points == [0]: 1035 ↛ 1049line 1035 didn't jump to line 1049 because the condition on line 1035 was always true

1036 ax.plot(

1037 xname,

1038 yfield,

1039 color=color,

1040 marker=line_marker,

1041 ls="-",

1042 lw=line_width,

1043 label=f"{label} (control)"

1044 if len(cube.coord(ensemble_coord).points) > 1

1045 else label,

1046 )

1047 # Label with (perturbed) if part of an ensemble and not the control.

1048 else:

1049 ax.plot(

1050 xname,

1051 yfield,

1052 color=color,

1053 ls="-",

1054 lw=1.5,

1055 alpha=0.75,

1056 label=f"{label} (member)",

1057 )

1058

1059 # Calculate the global min/max if multiple cubes are given.

1060 _, levels, _ = _colorbar_map_levels(cube, axis="y")

1061 if levels is not None: 1061 ↛ 1062line 1061 didn't jump to line 1062 because the condition on line 1061 was never true

1062 y_levels.append(min(levels))

1063 y_levels.append(max(levels))

1064

1065 # Add some labels and tweak the style.

1066 # check if cubes[0] works for single cube if not CubeList

1067

1068 # title = f"{title}\n [{coords.units.title(coords.points[0])}]"

1069 title = f"{title}"

1070 ax.set_title(title, fontsize=16)

1071

1072 # Set appropriate x-axis label based on coordinate

1073 if series_coordinate == "wavelength" or ( 1073 ↛ 1076line 1073 didn't jump to line 1076 because the condition on line 1073 was never true

1074 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"

1075 ):

1076 ax.set_xlabel("Wavelength (km)", fontsize=14)

1077 elif series_coordinate == "physical_wavenumber" or ( 1077 ↛ 1080line 1077 didn't jump to line 1080 because the condition on line 1077 was never true

1078 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"

1079 ):

1080 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

1081 else: # frequency or check units

1082 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1082 ↛ 1083line 1082 didn't jump to line 1083 because the condition on line 1082 was never true

1083 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

1084 else:

1085 ax.set_xlabel("Wavenumber", fontsize=14)

1086

1087 ax.set_ylabel("Power Spectral Density", fontsize=14)

1088 ax.tick_params(axis="both", labelsize=12)

1089

1090 # Set y limits to global min and max, autoscale if colorbar doesn't exist.

1091

1092 # Set log-log scale

1093 ax.set_xscale("log")

1094 ax.set_yscale("log")

1095

1096 # Add gridlines

1097 ax.grid(linestyle="--", color="grey", linewidth=1)

1098 # Ientify unique labels for legend

1099 handles = list(

1100 {

1101 label: handle

1102 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)

1103 }.values()

1104 )

1105 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)

1106

1107 # Save plot.

1108 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1109 logging.info("Saved line plot to %s", filename)

1110 plt.close(fig)

1111

1112

1113def _plot_and_save_vertical_line_series(

1114 cubes: iris.cube.CubeList,

1115 coords: list[iris.coords.Coord],

1116 ensemble_coord: str,

1117 filename: str,

1118 series_coordinate: str,

1119 title: str,

1120 vmin: float,

1121 vmax: float,

1122 **kwargs,

1123):

1124 """Plot and save a 1D line series in vertical.

1125

1126 Parameters

1127 ----------

1128 cubes: CubeList

1129 1 dimensional Cube or CubeList of the data to plot on x-axis.

1130 coord: list[Coord]

1131 Coordinates to plot on the y-axis, one per cube.

1132 ensemble_coord: str

1133 Ensemble coordinate in the cube.

1134 filename: str

1135 Filename of the plot to write.

1136 series_coordinate: str

1137 Coordinate to use as vertical axis.

1138 title: str

1139 Plot title.

1140 vmin: float

1141 Minimum value for the x-axis.

1142 vmax: float

1143 Maximum value for the x-axis.

1144 """

1145 # plot the vertical pressure axis using log scale

1146 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1147

1148 model_colors_map = _get_model_colors_map(cubes)

1149

1150 # Check match-up across sequence coords gives consistent sizes

1151 _validate_cubes_coords(cubes, coords)

1152

1153 for cube, coord in zip(cubes, coords, strict=True):

1154 label = None

1155 color = "black"

1156 if model_colors_map: 1156 ↛ 1157line 1156 didn't jump to line 1157 because the condition on line 1156 was never true

1157 label = cube.attributes.get("model_name")

1158 color = model_colors_map.get(label)

1159

1160 for cube_slice in cube.slices_over(ensemble_coord):

1161 # If ensemble data given plot control member with (control)

1162 # unless single forecast.

1163 if cube_slice.coord(ensemble_coord).points == [0]:

1164 iplt.plot(

1165 cube_slice,

1166 coord,

1167 color=color,

1168 marker="o",

1169 ls="-",

1170 lw=3,

1171 label=f"{label} (control)"

1172 if len(cube.coord(ensemble_coord).points) > 1

1173 else label,

1174 )

1175 # If ensemble data given plot perturbed members with (perturbed).

1176 else:

1177 iplt.plot(

1178 cube_slice,

1179 coord,

1180 color=color,

1181 ls="-",

1182 lw=1.5,

1183 alpha=0.75,

1184 label=f"{label} (member)",

1185 )

1186

1187 # Get the current axis

1188 ax = plt.gca()

1189

1190 # Special handling for pressure level data.

1191 if series_coordinate == "pressure": 1191 ↛ 1213line 1191 didn't jump to line 1213 because the condition on line 1191 was always true

1192 # Invert y-axis and set to log scale.

1193 ax.invert_yaxis()

1194 ax.set_yscale("log")

1195

1196 # Define y-ticks and labels for pressure log axis.

1197 y_tick_labels = [

1198 "1000",

1199 "850",

1200 "700",

1201 "500",

1202 "300",

1203 "200",

1204 "100",

1205 ]

1206 y_ticks = [1000, 850, 700, 500, 300, 200, 100]

1207

1208 # Set y-axis limits and ticks.

1209 ax.set_ylim(1100, 100)

1210

1211 # Test if series_coordinate is model level data. The UM data uses

1212 # model_level_number and lfric uses full_levels as coordinate.

1213 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):

1214 # Define y-ticks and labels for vertical axis.

1215 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points

1216 y_tick_labels = [str(int(i)) for i in y_ticks]

1217 ax.set_ylim(min(y_ticks), max(y_ticks))

1218

1219 ax.set_yticks(y_ticks)

1220 ax.set_yticklabels(y_tick_labels)

1221

1222 # Set x-axis limits.

1223 ax.set_xlim(vmin, vmax)

1224 # Mark y=0 if present in plot.

1225 if vmin < 0.0 and vmax > 0.0: 1225 ↛ 1226line 1225 didn't jump to line 1226 because the condition on line 1225 was never true

1226 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)

1227

1228 # Add some labels and tweak the style.

1229 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)

1230 ax.set_xlabel(

1231 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14

1232 )

1233 ax.set_title(title, fontsize=16)

1234 ax.ticklabel_format(axis="x")

1235 ax.tick_params(axis="y")

1236 ax.tick_params(axis="both", labelsize=12)

1237

1238 # Add gridlines

1239 ax.grid(linestyle="--", color="grey", linewidth=1)

1240 # Ientify unique labels for legend

1241 handles = list(

1242 {

1243 label: handle

1244 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)

1245 }.values()

1246 )

1247 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)

1248

1249 # Save plot.

1250 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1251 logging.info("Saved line plot to %s", filename)

1252 plt.close(fig)

1253

1254

1255def _plot_and_save_scatter_plot(

1256 cube_x: iris.cube.Cube | iris.cube.CubeList,

1257 cube_y: iris.cube.Cube | iris.cube.CubeList,

1258 filename: str,

1259 title: str,

1260 one_to_one: bool,

1261 model_names: list[str] = None,

1262 **kwargs,

1263):

1264 """Plot and save a 2D scatter plot.

1265

1266 Parameters

1267 ----------

1268 cube_x: Cube | CubeList

1269 1 dimensional Cube or CubeList of the data to plot on x-axis.

1270 cube_y: Cube | CubeList

1271 1 dimensional Cube or CubeList of the data to plot on y-axis.

1272 filename: str

1273 Filename of the plot to write.

1274 title: str

1275 Plot title.

1276 one_to_one: bool

1277 Whether a 1:1 line is plotted.

1278 """

1279 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1280 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures

1281 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate

1282 # over the pairs simultaneously.

1283

1284 # Ensure cube_x and cube_y are iterable

1285 cube_x_iterable = iter_maybe(cube_x)

1286 cube_y_iterable = iter_maybe(cube_y)

1287

1288 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):

1289 iplt.scatter(cube_x_iter, cube_y_iter)

1290 if one_to_one is True:

1291 plt.plot(

1292 [

1293 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),

1294 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),

1295 ],

1296 [

1297 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),

1298 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),

1299 ],

1300 "k",

1301 linestyle="--",

1302 )

1303 ax = plt.gca()

1304

1305 # Add some labels and tweak the style.

1306 if model_names is None:

1307 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)

1308 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)

1309 else:

1310 # Add the model names, these should be order of base (x) and other (y).

1311 ax.set_xlabel(

1312 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14

1313 )

1314 ax.set_ylabel(

1315 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14

1316 )

1317 ax.set_title(title, fontsize=16)

1318 ax.ticklabel_format(axis="y", useOffset=False)

1319 ax.tick_params(axis="x", labelrotation=15)

1320 ax.tick_params(axis="both", labelsize=12)

1321 ax.autoscale()

1322

1323 # Save plot.

1324 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1325 logging.info("Saved scatter plot to %s", filename)

1326 plt.close(fig)

1327

1328

1329def _plot_and_save_vector_plot(

1330 cube_u: iris.cube.Cube,

1331 cube_v: iris.cube.Cube,

1332 filename: str,

1333 title: str,

1334 method: Literal["contourf", "pcolormesh"],

1335 **kwargs,

1336):

1337 """Plot and save a 2D vector plot.

1338

1339 Parameters

1340 ----------

1341 cube_u: Cube

1342 2 dimensional Cube of u component of the data.

1343 cube_v: Cube

1344 2 dimensional Cube of v component of the data.

1345 filename: str

1346 Filename of the plot to write.

1347 title: str

1348 Plot title.

1349 """

1350 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1351

1352 # Create a cube containing the magnitude of the vector field.

1353 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5

1354 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")

1355

1356 # Specify the color bar

1357 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)

1358

1359 # Setup plot map projection, extent and coastlines and borderlines.

1360 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)

1361

1362 if method == "contourf":

1363 # Filled contour plot of the field.

1364 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)

1365 elif method == "pcolormesh":

1366 try:

1367 vmin = min(levels)

1368 vmax = max(levels)

1369 except TypeError:

1370 vmin, vmax = None, None

1371 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax

1372 # if levels are defined.

1373 if norm is not None:

1374 vmin = None

1375 vmax = None

1376 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)

1377 else:

1378 raise ValueError(f"Unknown plotting method: {method}")

1379

1380 # Check to see if transect, and if so, adjust y axis.

1381 if is_transect(cube_vec_mag):

1382 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:

1383 axes.invert_yaxis()

1384 axes.set_yscale("log")

1385 axes.set_ylim(1100, 100)

1386 # If both model_level_number and level_height exists, iplt can construct

1387 # plot as a function of height above orography (NOT sea level).

1388 elif {"model_level_number", "level_height"}.issubset(

1389 {coord.name() for coord in cube_vec_mag.coords()}

1390 ):

1391 axes.set_yscale("log")

1392

1393 axes.set_title(

1394 f"{title}\n"

1395 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"

1396 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"

1397 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"

1398 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",

1399 fontsize=16,

1400 )

1401

1402 else:

1403 # Add title.

1404 axes.set_title(title, fontsize=16)

1405

1406 # Add watermark with min/max/mean. Currently not user togglable.

1407 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.

1408 axes.annotate(

1409 f"Min: {np.min(cube_vec_mag.data):.3g} Max: {np.max(cube_vec_mag.data):.3g} Mean: {np.mean(cube_vec_mag.data):.3g}",

1410 xy=(1, -0.05),

1411 xycoords="axes fraction",

1412 xytext=(-5, 5),

1413 textcoords="offset points",

1414 ha="right",

1415 va="bottom",

1416 size=11,

1417 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),

1418 )

1419

1420 # Add colour bar.

1421 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)

1422 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)

1423 # add ticks and tick_labels for every levels if less than 20 levels exist

1424 if levels is not None and len(levels) < 20:

1425 cbar.set_ticks(levels)

1426 cbar.set_ticklabels([f"{level:.1f}" for level in levels])

1427

1428 # 30 barbs along the longest axis of the plot, or a barb per point for data

1429 # with less than 30 points.

1430 step = max(max(cube_u.shape) // 30, 1)

1431 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")

1432

1433 # Save plot.

1434 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1435 logging.info("Saved vector plot to %s", filename)

1436 plt.close(fig)

1437

1438

1439def _plot_and_save_histogram_series(

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

1441 filename: str,

1442 title: str,

1443 vmin: float,

1444 vmax: float,

1445 **kwargs,

1446):

1447 """Plot and save a histogram series.

1448

1449 Parameters

1450 ----------

1451 cubes: Cube or CubeList

1452 2 dimensional Cube or CubeList of the data to plot as histogram.

1453 filename: str

1454 Filename of the plot to write.

1455 title: str

1456 Plot title.

1457 vmin: float

1458 minimum for colorbar

1459 vmax: float

1460 maximum for colorbar

1461 """

1462 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1463 ax = plt.gca()

1464

1465 model_colors_map = _get_model_colors_map(cubes)

1466

1467 # Set default that histograms will produce probability density function

1468 # at each bin (integral over range sums to 1).

1469 density = True

1470

1471 for cube in iter_maybe(cubes):

1472 # Easier to check title (where var name originates)

1473 # than seeing if long names exist etc.

1474 # Exception case, where distribution better fits log scales/bins.

1475 if "surface_microphysical" in title:

1476 if "amount" in title: 1476 ↛ 1478line 1476 didn't jump to line 1478 because the condition on line 1476 was never true

1477 # Compute histogram following Klingaman et al. (2017): ASoP

1478 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))

1479 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)

1480 density = False

1481 else:

1482 bins = 10.0 ** (

1483 np.arange(-10, 27, 1) / 10.0

1484 ) # Suggestion from RMED toolbox.

1485 bins = np.insert(bins, 0, 0)

1486 ax.set_yscale("log")

1487 vmin = bins[1]

1488 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.

1489 ax.set_xscale("log")

1490 elif "lightning" in title:

1491 bins = [0, 1, 2, 3, 4, 5]

1492 else:

1493 bins = np.linspace(vmin, vmax, 51)

1494 logging.debug(

1495 "Plotting histogram with %s bins %s - %s.",

1496 np.size(bins),

1497 np.min(bins),

1498 np.max(bins),

1499 )

1500

1501 # Reshape cube data into a single array to allow for a single histogram.

1502 # Otherwise we plot xdim histograms stacked.

1503 cube_data_1d = (cube.data).flatten()

1504

1505 label = None

1506 color = "black"

1507 if model_colors_map: 1507 ↛ 1508line 1507 didn't jump to line 1508 because the condition on line 1507 was never true

1508 label = cube.attributes.get("model_name")

1509 color = model_colors_map[label]

1510 x, y = np.histogram(cube_data_1d, bins=bins, density=density)

1511

1512 # Compute area under curve.

1513 if "surface_microphysical" in title and "amount" in title: 1513 ↛ 1514line 1513 didn't jump to line 1514 because the condition on line 1513 was never true

1514 bin_mean = (bins[:-1] + bins[1:]) / 2.0

1515 x = x * bin_mean / x.sum()

1516 x = x[1:]

1517 y = y[1:]

1518

1519 ax.plot(

1520 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label

1521 )

1522

1523 # Add some labels and tweak the style.

1524 ax.set_title(title, fontsize=16)

1525 ax.set_xlabel(

1526 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14

1527 )

1528 ax.set_ylabel("Normalised probability density", fontsize=14)

1529 if "surface_microphysical" in title and "amount" in title: 1529 ↛ 1530line 1529 didn't jump to line 1530 because the condition on line 1529 was never true

1530 ax.set_ylabel(

1531 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14

1532 )

1533 ax.set_xlim(vmin, vmax)

1534 ax.tick_params(axis="both", labelsize=12)

1535

1536 # Overlay grid-lines onto histogram plot.

1537 ax.grid(linestyle="--", color="grey", linewidth=1)

1538 if model_colors_map: 1538 ↛ 1539line 1538 didn't jump to line 1539 because the condition on line 1538 was never true

1539 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)

1540

1541 # Save plot.

1542 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1543 logging.info("Saved histogram plot to %s", filename)

1544 plt.close(fig)

1545

1546

1547def _plot_and_save_postage_stamp_histogram_series(

1548 cube: iris.cube.Cube,

1549 filename: str,

1550 title: str,

1551 stamp_coordinate: str,

1552 vmin: float,

1553 vmax: float,

1554 **kwargs,

1555):

1556 """Plot and save postage (ensemble members) stamps for a histogram series.

1557

1558 Parameters

1559 ----------

1560 cube: Cube

1561 2 dimensional Cube of the data to plot as histogram.

1562 filename: str

1563 Filename of the plot to write.

1564 title: str

1565 Plot title.

1566 stamp_coordinate: str

1567 Coordinate that becomes different plots.

1568 vmin: float

1569 minimum for pdf x-axis

1570 vmax: float

1571 maximum for pdf x-axis

1572 """

1573 # Use the smallest square grid that will fit the members.

1574 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))

1575

1576 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1577 # Make a subplot for each member.

1578 for member, subplot in zip(

1579 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False

1580 ):

1581 # Implicit interface is much easier here, due to needing to have the

1582 # cartopy GeoAxes generated.

1583 plt.subplot(grid_size, grid_size, subplot)

1584 # Reshape cube data into a single array to allow for a single histogram.

1585 # Otherwise we plot xdim histograms stacked.

1586 member_data_1d = (member.data).flatten()

1587 plt.hist(member_data_1d, density=True, stacked=True)

1588 ax = plt.gca()

1589 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")

1590 ax.set_xlim(vmin, vmax)

1591

1592 # Overall figure title.

1593 fig.suptitle(title, fontsize=16)

1594

1595 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1596 logging.info("Saved histogram postage stamp plot to %s", filename)

1597 plt.close(fig)

1598

1599

1600def _plot_and_save_postage_stamps_in_single_plot_histogram_series(

1601 cube: iris.cube.Cube,

1602 filename: str,

1603 title: str,

1604 stamp_coordinate: str,

1605 vmin: float,

1606 vmax: float,

1607 **kwargs,

1608):

1609 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")

1610 ax.set_title(title, fontsize=16)

1611 ax.set_xlim(vmin, vmax)

1612 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)

1613 ax.set_ylabel("normalised probability density", fontsize=14)

1614 # Loop over all slices along the stamp_coordinate

1615 for member in cube.slices_over(stamp_coordinate):

1616 # Flatten the member data to 1D

1617 member_data_1d = member.data.flatten()

1618 # Plot the histogram using plt.hist

1619 plt.hist(

1620 member_data_1d,

1621 density=True,

1622 stacked=True,

1623 label=f"Member #{member.coord(stamp_coordinate).points[0]}",

1624 )

1625

1626 # Add a legend

1627 ax.legend(fontsize=16)

1628

1629 # Save the figure to a file

1630 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1631 logging.info("Saved histogram postage stamp plot to %s", filename)

1632

1633 # Close the figure

1634 plt.close(fig)

1635

1636

1637def _plot_and_save_scattermap_plot(

1638 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs

1639):

1640 """Plot and save a geographical scatter plot.

1641

1642 Parameters

1643 ----------

1644 cube: Cube

1645 1 dimensional Cube of the data points with auxiliary latitude and

1646 longitude coordinates,

1647 filename: str

1648 Filename of the plot to write.

1649 title: str

1650 Plot title.

1651 projection: str

1652 Mapping projection to be used by cartopy.

1653 """

1654 # Setup plot details, size, resolution, etc.

1655 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

1656 if projection is not None:

1657 # Apart from the default, the only projection we currently support is

1658 # a stereographic projection over the North Pole.

1659 if projection == "NP_Stereo":

1660 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))

1661 else:

1662 raise ValueError(f"Unknown projection: {projection}")

1663 else:

1664 axes = plt.axes(projection=ccrs.PlateCarree())

1665

1666 # Scatter plot of the field. The marker size is chosen to give

1667 # symbols that decrease in size as the number of observations

1668 # increases, although the fraction of the figure covered by

1669 # symbols increases roughly as N^(1/2), disregarding overlaps,

1670 # and has been selected for the default figure size of (10, 10).

1671 # Should this be changed, the marker size should be adjusted in

1672 # proportion to the area of the figure.

1673 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))

1674 klon = None

1675 klat = None

1676 for kc in range(len(cube.aux_coords)):

1677 if cube.aux_coords[kc].standard_name == "latitude":

1678 klat = kc

1679 elif cube.aux_coords[kc].standard_name == "longitude":

1680 klon = kc

1681 scatter_map = iplt.scatter(

1682 cube.aux_coords[klon],

1683 cube.aux_coords[klat],

1684 c=cube.data[:],

1685 s=mrk_size,

1686 cmap="jet",

1687 edgecolors="k",

1688 )

1689

1690 # Add coastlines and borderlines.

1691 try:

1692 axes.coastlines(resolution="10m")

1693 axes.add_feature(cfeature.BORDERS)

1694 except AttributeError:

1695 pass

1696

1697 # Add title.

1698 axes.set_title(title, fontsize=16)

1699

1700 # Add colour bar.

1701 cbar = fig.colorbar(scatter_map)

1702 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)

1703

1704 # Save plot.

1705 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

1706 logging.info("Saved geographical scatter plot to %s", filename)

1707 plt.close(fig)

1708

1709

1710def _spatial_plot(

1711 method: Literal["contourf", "pcolormesh"],

1712 cube: iris.cube.Cube,

1713 filename: str | None,

1714 sequence_coordinate: str,

1715 stamp_coordinate: str,

1716 overlay_cube: iris.cube.Cube | None = None,

1717 contour_cube: iris.cube.Cube | None = None,

1718 **kwargs,

1719):

1720 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.

1721

1722 A 2D spatial field can be plotted, but if the sequence_coordinate is present

1723 then a sequence of plots will be produced. Similarly if the stamp_coordinate

1724 is present then postage stamp plots will be produced.

1725

1726 If an overlay_cube and/or contour_cube are specified, multiple variables can

1727 be overplotted on the same figure.

1728

1729 Parameters

1730 ----------

1731 method: "contourf" | "pcolormesh"

1732 The plotting method to use.

1733 cube: Cube

1734 Iris cube of the data to plot. It should have two spatial dimensions,

1735 such as lat and lon, and may also have a another two dimension to be

1736 plotted sequentially and/or as postage stamp plots.

1737 filename: str | None

1738 Name of the plot to write, used as a prefix for plot sequences. If None

1739 uses the recipe name.

1740 sequence_coordinate: str

1741 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

1742 This coordinate must exist in the cube.

1743 stamp_coordinate: str

1744 Coordinate about which to plot postage stamp plots. Defaults to

1745 ``"realization"``.

1746 overlay_cube: Cube | None, optional

1747 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube

1748 contour_cube: Cube | None, optional

1749 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube

1750

1751 Raises

1752 ------

1753 ValueError

1754 If the cube doesn't have the right dimensions.

1755 TypeError

1756 If the cube isn't a single cube.

1757 """

1758 recipe_title = get_recipe_metadata().get("title", "Untitled")

1759

1760 # Ensure we've got a single cube.

1761 cube = _check_single_cube(cube)

1762

1763 # Make postage stamp plots if stamp_coordinate exists and has more than a

1764 # single point.

1765 plotting_func = _plot_and_save_spatial_plot

1766 try:

1767 if cube.coord(stamp_coordinate).shape[0] > 1:

1768 plotting_func = _plot_and_save_postage_stamp_spatial_plot

1769 except iris.exceptions.CoordinateNotFoundError:

1770 pass

1771

1772 # Produce a geographical scatter plot if the data have a

1773 # dimension called observation or model_obs_error

1774 if any( 1774 ↛ 1778line 1774 didn't jump to line 1778 because the condition on line 1774 was never true

1775 crd.var_name == "station" or crd.var_name == "model_obs_error"

1776 for crd in cube.coords()

1777 ):

1778 plotting_func = _plot_and_save_scattermap_plot

1779

1780 # Must have a sequence coordinate.

1781 try:

1782 cube.coord(sequence_coordinate)

1783 except iris.exceptions.CoordinateNotFoundError as err:

1784 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err

1785

1786 # Create a plot for each value of the sequence coordinate.

1787 plot_index = []

1788 nplot = np.size(cube.coord(sequence_coordinate).points)

1789

1790 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):

1791 # Set plot titles and filename

1792 seq_coord = cube_slice.coord(sequence_coordinate)

1793 plot_title, plot_filename = _set_title_and_filename(

1794 seq_coord, nplot, recipe_title, filename

1795 )

1796

1797 # Extract sequence slice for overlay_cube and contour_cube if required.

1798 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)

1799 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)

1800

1801 # Do the actual plotting.

1802 plotting_func(

1803 cube_slice,

1804 filename=plot_filename,

1805 stamp_coordinate=stamp_coordinate,

1806 title=plot_title,

1807 method=method,

1808 overlay_cube=overlay_slice,

1809 contour_cube=contour_slice,

1810 **kwargs,

1811 )

1812 plot_index.append(plot_filename)

1813

1814 # Add list of plots to plot metadata.

1815 complete_plot_index = _append_to_plot_index(plot_index)

1816

1817 # Make a page to display the plots.

1818 _make_plot_html_page(complete_plot_index)

1819

1820

1821def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):

1822 """Get colourmap for mask.

1823

1824 If "mask_for_" appears anywhere in the name of a cube this function will be called

1825 regardless of the name of the variable to ensure a consistent plot.

1826

1827 Parameters

1828 ----------

1829 cube: Cube

1830 Cube of variable for which the colorbar information is desired.

1831 axis: "x", "y", optional

1832 Select the levels for just this axis of a line plot. The min and max

1833 can be set by xmin/xmax or ymin/ymax respectively. For variables where

1834 setting a universal range is not desirable (e.g. temperature), users

1835 can set ymin/ymax values to "auto" in the colorbar definitions file.

1836 Where no additional xmin/xmax or ymin/ymax values are provided, the

1837 axis bounds default to use the vmin/vmax values provided.

1838

1839 Returns

1840 -------

1841 cmap:

1842 Matplotlib colormap.

1843 levels:

1844 List of levels to use for plotting. For continuous plots the min and max

1845 should be taken as the range.

1846 norm:

1847 BoundaryNorm information.

1848 """

1849 if "difference" not in cube.long_name:

1850 if axis:

1851 levels = [0, 1]

1852 # Complete settings based on levels.

1853 return None, levels, None

1854 else:

1855 # Define the levels and colors.

1856 levels = [0, 1, 2]

1857 colors = ["white", "dodgerblue"]

1858 # Create a custom color map.

1859 cmap = mcolors.ListedColormap(colors)

1860 # Normalize the levels.

1861 norm = mcolors.BoundaryNorm(levels, cmap.N)

1862 logging.debug("Colourmap for %s.", cube.long_name)

1863 return cmap, levels, norm

1864 else:

1865 if axis:

1866 levels = [-1, 1]

1867 return None, levels, None

1868 else:

1869 # Search for if mask difference, set to +/- 0.5 as values plotted <

1870 # not <=.

1871 levels = [-2, -0.5, 0.5, 2]

1872 colors = ["goldenrod", "white", "teal"]

1873 cmap = mcolors.ListedColormap(colors)

1874 norm = mcolors.BoundaryNorm(levels, cmap.N)

1875 logging.debug("Colourmap for %s.", cube.long_name)

1876 return cmap, levels, norm

1877

1878

1879def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):

1880 """Get a custom colorbar for a cube in the Beaufort Scale.

1881

1882 Specific variable ranges can be separately set in user-supplied definition

1883 for x- or y-axis limits, or indicate where automated range preferred.

1884

1885 Parameters

1886 ----------

1887 cube: Cube

1888 Cube of variable with Beaufort Scale in name.

1889 axis: "x", "y", optional

1890 Select the levels for just this axis of a line plot. The min and max

1891 can be set by xmin/xmax or ymin/ymax respectively. For variables where

1892 setting a universal range is not desirable (e.g. temperature), users

1893 can set ymin/ymax values to "auto" in the colorbar definitions file.

1894 Where no additional xmin/xmax or ymin/ymax values are provided, the

1895 axis bounds default to use the vmin/vmax values provided.

1896

1897 Returns

1898 -------

1899 cmap:

1900 Matplotlib colormap.

1901 levels:

1902 List of levels to use for plotting. For continuous plots the min and max

1903 should be taken as the range.

1904 norm:

1905 BoundaryNorm information.

1906 """

1907 if "difference" not in cube.long_name:

1908 if axis:

1909 levels = [0, 12]

1910 return None, levels, None

1911 else:

1912 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]

1913 colors = [

1914 "black",

1915 (0, 0, 0.6),

1916 "blue",

1917 "cyan",

1918 "green",

1919 "yellow",

1920 (1, 0.5, 0),

1921 "red",

1922 "pink",

1923 "magenta",

1924 "purple",

1925 "maroon",

1926 "white",

1927 ]

1928 cmap = mcolors.ListedColormap(colors)

1929 norm = mcolors.BoundaryNorm(levels, cmap.N)

1930 logging.info("change colormap for Beaufort Scale colorbar.")

1931 return cmap, levels, norm

1932 else:

1933 if axis:

1934 levels = [-4, 4]

1935 return None, levels, None

1936 else:

1937 levels = [

1938 -3.5,

1939 -2.5,

1940 -1.5,

1941 -0.5,

1942 0.5,

1943 1.5,

1944 2.5,

1945 3.5,

1946 ]

1947 cmap = plt.get_cmap("bwr", 8)

1948 norm = mcolors.BoundaryNorm(levels, cmap.N)

1949 return cmap, levels, norm

1950

1951

1952def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):

1953 """Return altered colourmap for temperature with change in units to Celsius.

1954

1955 If "Celsius" appears anywhere in the name of a cube this function will be called.

1956

1957 Parameters

1958 ----------

1959 cube: Cube

1960 Cube of variable for which the colorbar information is desired.

1961 cmap: Matplotlib colormap.

1962 levels: List

1963 List of levels to use for plotting. For continuous plots the min and max

1964 should be taken as the range.

1965 norm: BoundaryNorm.

1966

1967 Returns

1968 -------

1969 cmap: Matplotlib colormap.

1970 levels: List

1971 List of levels to use for plotting. For continuous plots the min and max

1972 should be taken as the range.

1973 norm: BoundaryNorm.

1974 """

1975 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])

1976 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:

1977 levels = np.array(levels)

1978 levels -= 273

1979 levels = levels.tolist()

1980 else:

1981 # Do nothing keep the existing colourbar attributes

1982 levels = levels

1983 cmap = cmap

1984 norm = norm

1985 return cmap, levels, norm

1986

1987

1988def _custom_colormap_probability(

1989 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None

1990):

1991 """Get a custom colorbar for a probability cube.

1992

1993 Specific variable ranges can be separately set in user-supplied definition

1994 for x- or y-axis limits, or indicate where automated range preferred.

1995

1996 Parameters

1997 ----------

1998 cube: Cube

1999 Cube of variable with probability in name.

2000 axis: "x", "y", optional

2001 Select the levels for just this axis of a line plot. The min and max

2002 can be set by xmin/xmax or ymin/ymax respectively. For variables where

2003 setting a universal range is not desirable (e.g. temperature), users

2004 can set ymin/ymax values to "auto" in the colorbar definitions file.

2005 Where no additional xmin/xmax or ymin/ymax values are provided, the

2006 axis bounds default to use the vmin/vmax values provided.

2007

2008 Returns

2009 -------

2010 cmap:

2011 Matplotlib colormap.

2012 levels:

2013 List of levels to use for plotting. For continuous plots the min and max

2014 should be taken as the range.

2015 norm:

2016 BoundaryNorm information.

2017 """

2018 if axis:

2019 levels = [0, 1]

2020 return None, levels, None

2021 else:

2022 cmap = mcolors.ListedColormap(

2023 [

2024 "#FFFFFF",

2025 "#636363",

2026 "#e1dada",

2027 "#B5CAFF",

2028 "#8FB3FF",

2029 "#7F97FF",

2030 "#ABCF63",

2031 "#E8F59E",

2032 "#FFFA14",

2033 "#FFD121",

2034 "#FFA30A",

2035 ]

2036 )

2037 levels = [0.0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]

2038 norm = mcolors.BoundaryNorm(levels, cmap.N)

2039 return cmap, levels, norm

2040

2041

2042def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):

2043 """Return a custom colourmap for the current recipe."""

2044 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])

2045 if (

2046 any("surface_microphysical" in name for name in varnames)

2047 and "difference" not in cube.long_name

2048 and "mask" not in cube.long_name

2049 ):

2050 # Define the levels and colors

2051 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]

2052 colors = [

2053 "w",

2054 (0, 0, 0.6),

2055 "b",

2056 "c",

2057 "g",

2058 "y",

2059 (1, 0.5, 0),

2060 "r",

2061 "pink",

2062 "m",

2063 "purple",

2064 "maroon",

2065 "gray",

2066 ]

2067 # Create a custom colormap

2068 cmap = mcolors.ListedColormap(colors)

2069 # Normalize the levels

2070 norm = mcolors.BoundaryNorm(levels, cmap.N)

2071 logging.info("change colormap for surface_microphysical variable colorbar.")

2072 else:

2073 # do nothing and keep existing colorbar attributes

2074 cmap = cmap

2075 levels = levels

2076 norm = norm

2077 return cmap, levels, norm

2078

2079

2080def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):

2081 """Return custom colourmap for aviation colour state.

2082

2083 If "aviation_colour_state" appears anywhere in the name of a cube

2084 this function will be called.

2085

2086 Parameters

2087 ----------

2088 cube: Cube

2089 Cube of variable for which the colorbar information is desired.

2090

2091 Returns

2092 -------

2093 cmap: Matplotlib colormap.

2094 levels: List

2095 List of levels to use for plotting. For continuous plots the min and max

2096 should be taken as the range.

2097 norm: BoundaryNorm.

2098 """

2099 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]

2100 colors = [

2101 "#87ceeb",

2102 "#ffffff",

2103 "#8ced69",

2104 "#ffff00",

2105 "#ffd700",

2106 "#ffa500",

2107 "#fe3620",

2108 ]

2109 # Create a custom colormap

2110 cmap = mcolors.ListedColormap(colors)

2111 # Normalise the levels

2112 norm = mcolors.BoundaryNorm(levels, cmap.N)

2113 return cmap, levels, norm

2114

2115

2116def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):

2117 """Return a custom colourmap for the current recipe."""

2118 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])

2119 if (

2120 any("visibility_in_air" in name for name in varnames)

2121 and "difference" not in cube.long_name

2122 and "mask" not in cube.long_name

2123 ):

2124 # Define the levels and colors (in km)

2125 levels = [0, 0.05, 0.1, 0.2, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 50.0, 70.0, 100.0]

2126 norm = mcolors.BoundaryNorm(levels, cmap.N)

2127 colours = [

2128 "#8f00d6",

2129 "#d10000",

2130 "#ff9700",

2131 "#ffff00",

2132 "#00007f",

2133 "#6c9ccd",

2134 "#aae8ff",

2135 "#37a648",

2136 "#8edc64",

2137 "#c5ffc5",

2138 "#dcdcdc",

2139 "#ffffff",

2140 ]

2141 # Create a custom colormap

2142 cmap = mcolors.ListedColormap(colours)

2143 # Normalize the levels

2144 norm = mcolors.BoundaryNorm(levels, cmap.N)

2145 logging.info("change colormap for visibility_in_air variable colorbar.")

2146 else:

2147 # do nothing and keep existing colorbar attributes

2148 cmap = cmap

2149 levels = levels

2150 norm = norm

2151 return cmap, levels, norm

2152

2153

2154def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:

2155 """Return number of models based on cube attributes."""

2156 model_names = list(

2157 filter(

2158 lambda x: x is not None,

2159 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},

2160 )

2161 )

2162 if not model_names:

2163 logging.debug("Missing model names. Will assume single model.")

2164 return 1

2165 else:

2166 return len(model_names)

2167

2168

2169def _validate_cube_shape(

2170 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int

2171) -> None:

2172 """Check all cubes have a model name."""

2173 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2173 ↛ 2174line 2173 didn't jump to line 2174 because the condition on line 2173 was never true

2174 raise ValueError(

2175 f"The number of model names ({num_models}) should equal the number "

2176 f"of cubes ({len(cube)})."

2177 )

2178

2179

2180def _validate_cubes_coords(

2181 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]

2182) -> None:

2183 """Check same number of cubes as sequence coordinate for zip functions."""

2184 if len(cubes) != len(coords): 2184 ↛ 2185line 2184 didn't jump to line 2185 because the condition on line 2184 was never true

2185 raise ValueError(

2186 f"The number of CubeList entries ({len(cubes)}) should equal the number "

2187 f"of sequence coordinates ({len(coords)})."

2188 f"Check that number of time entries in input data are consistent if "

2189 f"performing time-averaging steps prior to plotting outputs."

2190 )

2191

2192

2193####################

2194# Public functions #

2195####################

2196

2197

2198def spatial_contour_plot(

2199 cube: iris.cube.Cube,

2200 filename: str = None,

2201 sequence_coordinate: str = "time",

2202 stamp_coordinate: str = "realization",

2203 **kwargs,

2204) -> iris.cube.Cube:

2205 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.

2206

2207 A 2D spatial field can be plotted, but if the sequence_coordinate is present

2208 then a sequence of plots will be produced. Similarly if the stamp_coordinate

2209 is present then postage stamp plots will be produced.

2210

2211 Parameters

2212 ----------

2213 cube: Cube

2214 Iris cube of the data to plot. It should have two spatial dimensions,

2215 such as lat and lon, and may also have a another two dimension to be

2216 plotted sequentially and/or as postage stamp plots.

2217 filename: str, optional

2218 Name of the plot to write, used as a prefix for plot sequences. Defaults

2219 to the recipe name.

2220 sequence_coordinate: str, optional

2221 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

2222 This coordinate must exist in the cube.

2223 stamp_coordinate: str, optional

2224 Coordinate about which to plot postage stamp plots. Defaults to

2225 ``"realization"``.

2226

2227 Returns

2228 -------

2229 Cube

2230 The original cube (so further operations can be applied).

2231

2232 Raises

2233 ------

2234 ValueError

2235 If the cube doesn't have the right dimensions.

2236 TypeError

2237 If the cube isn't a single cube.

2238 """

2239 _spatial_plot(

2240 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs

2241 )

2242 return cube

2243

2244

2245def spatial_pcolormesh_plot(

2246 cube: iris.cube.Cube,

2247 filename: str = None,

2248 sequence_coordinate: str = "time",

2249 stamp_coordinate: str = "realization",

2250 **kwargs,

2251) -> iris.cube.Cube:

2252 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.

2253

2254 A 2D spatial field can be plotted, but if the sequence_coordinate is present

2255 then a sequence of plots will be produced. Similarly if the stamp_coordinate

2256 is present then postage stamp plots will be produced.

2257

2258 This function is significantly faster than ``spatial_contour_plot``,

2259 especially at high resolutions, and should be preferred unless contiguous

2260 contour areas are important.

2261

2262 Parameters

2263 ----------

2264 cube: Cube

2265 Iris cube of the data to plot. It should have two spatial dimensions,

2266 such as lat and lon, and may also have a another two dimension to be

2267 plotted sequentially and/or as postage stamp plots.

2268 filename: str, optional

2269 Name of the plot to write, used as a prefix for plot sequences. Defaults

2270 to the recipe name.

2271 sequence_coordinate: str, optional

2272 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

2273 This coordinate must exist in the cube.

2274 stamp_coordinate: str, optional

2275 Coordinate about which to plot postage stamp plots. Defaults to

2276 ``"realization"``.

2277

2278 Returns

2279 -------

2280 Cube

2281 The original cube (so further operations can be applied).

2282

2283 Raises

2284 ------

2285 ValueError

2286 If the cube doesn't have the right dimensions.

2287 TypeError

2288 If the cube isn't a single cube.

2289 """

2290 _spatial_plot(

2291 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs

2292 )

2293 return cube

2294

2295

2296def spatial_multi_pcolormesh_plot(

2297 cube: iris.cube.Cube,

2298 overlay_cube: iris.cube.Cube,

2299 contour_cube: iris.cube.Cube,

2300 filename: str = None,

2301 sequence_coordinate: str = "time",

2302 stamp_coordinate: str = "realization",

2303 **kwargs,

2304) -> iris.cube.Cube:

2305 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.

2306

2307 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present

2308 then a sequence of plots will be produced. Similarly if the stamp_coordinate

2309 is present then postage stamp plots will be produced.

2310

2311 If specified, a masked overlay_cube can be overplotted on top of the base cube.

2312

2313 If specified, contours of a contour_cube can be overplotted on top of those.

2314

2315 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.

2316

2317 This function is significantly faster than ``spatial_contour_plot``,

2318 especially at high resolutions, and should be preferred unless contiguous

2319 contour areas are important.

2320

2321 Parameters

2322 ----------

2323 cube: Cube

2324 Iris cube of the data to plot. It should have two spatial dimensions,

2325 such as lat and lon, and may also have a another two dimension to be

2326 plotted sequentially and/or as postage stamp plots.

2327 overlay_cube: Cube

2328 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,

2329 such as lat and lon, and may also have a another two dimension to be

2330 plotted sequentially and/or as postage stamp plots. This is likely to be a masked cube in order not to hide the underlying basis cube.

2331 contour_cube: Cube

2332 Iris cube of the data to plot as a contour overlay on top of basis cube and overlay_cube. It should have two spatial dimensions,

2333 such as lat and lon, and may also have a another two dimension to be

2334 plotted sequentially and/or as postage stamp plots.

2335 filename: str, optional

2336 Name of the plot to write, used as a prefix for plot sequences. Defaults

2337 to the recipe name.

2338 sequence_coordinate: str, optional

2339 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

2340 This coordinate must exist in the cube.

2341 stamp_coordinate: str, optional

2342 Coordinate about which to plot postage stamp plots. Defaults to

2343 ``"realization"``.

2344

2345 Returns

2346 -------

2347 Cube

2348 The original cube (so further operations can be applied).

2349

2350 Raises

2351 ------

2352 ValueError

2353 If the cube doesn't have the right dimensions.

2354 TypeError

2355 If the cube isn't a single cube.

2356 """

2357 _spatial_plot(

2358 "pcolormesh",

2359 cube,

2360 filename,

2361 sequence_coordinate,

2362 stamp_coordinate,

2363 overlay_cube=overlay_cube,

2364 contour_cube=contour_cube,

2365 )

2366 return cube, overlay_cube, contour_cube

2367

2368

2369# TODO: Expand function to handle ensemble data.

2370# line_coordinate: str, optional

2371# Coordinate about which to plot multiple lines. Defaults to

2372# ``"realization"``.

2373def plot_line_series(

2374 cube: iris.cube.Cube | iris.cube.CubeList,

2375 filename: str = None,

2376 series_coordinate: str = "time",

2377 sequence_coordinate: str = "time",

2378 # add the following for ensembles

2379 stamp_coordinate: str = "realization",

2380 single_plot: bool = False,

2381 **kwargs,

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

2383 """Plot a line plot for the specified coordinate.

2384

2385 The Cube or CubeList must be 1D.

2386

2387 Parameters

2388 ----------

2389 iris.cube | iris.cube.CubeList

2390 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.

2391 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.

2392 We do not support different data such as temperature and humidity in the same CubeList for plotting.

2393 filename: str, optional

2394 Name of the plot to write, used as a prefix for plot sequences. Defaults

2395 to the recipe name.

2396 series_coordinate: str, optional

2397 Coordinate about which to make a series. Defaults to ``"time"``. This

2398 coordinate must exist in the cube.

2399

2400 Returns

2401 -------

2402 iris.cube.Cube | iris.cube.CubeList

2403 The original Cube or CubeList (so further operations can be applied).

2404 plotted data.

2405

2406 Raises

2407 ------

2408 ValueError

2409 If the cubes don't have the right dimensions.

2410 TypeError

2411 If the cube isn't a Cube or CubeList.

2412 """

2413 stamp_coordinate = "realization"

2414 # Ensure we have a name for the plot file.

2415 recipe_title = get_recipe_metadata().get("title", "Untitled")

2416

2417 num_models = _get_num_models(cube)

2418

2419 _validate_cube_shape(cube, num_models)

2420

2421 # Iterate over all cubes and extract coordinate to plot.

2422 cubes = iter_maybe(cube)

2423

2424 coords = []

2425 for cube in cubes:

2426 try:

2427 coords.append(cube.coord(series_coordinate))

2428 except iris.exceptions.CoordinateNotFoundError as err:

2429 raise ValueError(

2430 f"Cube must have a {series_coordinate} coordinate."

2431 ) from err

2432 if cube.coords("realization"): 2432 ↛ 2436line 2432 didn't jump to line 2436 because the condition on line 2432 was always true

2433 if cube.ndim > 3: 2433 ↛ 2434line 2433 didn't jump to line 2434 because the condition on line 2433 was never true

2434 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")

2435 else:

2436 raise ValueError("Cube must have a realization coordinate.")

2437

2438 plot_index = []

2439

2440 # Check if this is a spectral plot by looking for spectral coordinates

2441 is_spectral_plot = series_coordinate in [

2442 "frequency",

2443 "physical_wavenumber",

2444 "wavelength",

2445 ]

2446

2447 if is_spectral_plot:

2448 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series

2449 # coordinate frequency/wavenumber.

2450 # If several power spectra are plotted with time as sequence_coordinate for the

2451 # time slider option.

2452

2453 # Internal plotting function.

2454 plotting_func = _plot_and_save_line_power_spectrum_series

2455

2456 for cube in cubes:

2457 try:

2458 cube.coord(sequence_coordinate)

2459 except iris.exceptions.CoordinateNotFoundError as err:

2460 raise ValueError(

2461 f"Cube must have a {sequence_coordinate} coordinate."

2462 ) from err

2463

2464 if num_models == 1: 2464 ↛ 2478line 2464 didn't jump to line 2478 because the condition on line 2464 was always true

2465 # check for ensembles

2466 if ( 2466 ↛ 2470line 2466 didn't jump to line 2470 because the condition on line 2466 was never true

2467 stamp_coordinate in [c.name() for c in cubes[0].coords()]

2468 and cubes[0].coord(stamp_coordinate).shape[0] > 1

2469 ):

2470 if single_plot:

2471 # Plot spectra, mean and ensemble spread on 1 plot

2472 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series

2473 else:

2474 # Plot postage stamps

2475 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series

2476 cube_iterables = cubes[0].slices_over(sequence_coordinate)

2477 else:

2478 all_points = sorted(

2479 set(

2480 itertools.chain.from_iterable(

2481 cb.coord(sequence_coordinate).points for cb in cubes

2482 )

2483 )

2484 )

2485 all_slices = list(

2486 itertools.chain.from_iterable(

2487 cb.slices_over(sequence_coordinate) for cb in cubes

2488 )

2489 )

2490 # Matched slices (matched by seq coord point; it may happen that

2491 # evaluated models do not cover the same seq coord range, hence matching

2492 # necessary)

2493 cube_iterables = [

2494 iris.cube.CubeList(

2495 s

2496 for s in all_slices

2497 if s.coord(sequence_coordinate).points[0] == point

2498 )

2499 for point in all_points

2500 ]

2501

2502 nplot = np.size(cube.coord(sequence_coordinate).points)

2503

2504 # Create a plot for each value of the sequence coordinate. Allowing for

2505 # multiple cubes in a CubeList to be plotted in the same plot for similar

2506 # sequence values. Passing a CubeList into the internal plotting function

2507 # for similar values of the sequence coordinate. cube_slice can be an

2508 # iris.cube.Cube or an iris.cube.CubeList.

2509

2510 for cube_slice in cube_iterables:

2511 # Normalize cube_slice to a list of cubes

2512 if isinstance(cube_slice, iris.cube.CubeList): 2512 ↛ 2513line 2512 didn't jump to line 2513 because the condition on line 2512 was never true

2513 cubes = list(cube_slice)

2514 elif isinstance(cube_slice, iris.cube.Cube): 2514 ↛ 2517line 2514 didn't jump to line 2517 because the condition on line 2514 was always true

2515 cubes = [cube_slice]

2516 else:

2517 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")

2518

2519 # Use sequence value so multiple sequences can merge.

2520 seq_coord = cube_slice[0].coord(sequence_coordinate)

2521 plot_title, plot_filename = _set_title_and_filename(

2522 seq_coord, nplot, recipe_title, filename

2523 )

2524

2525 # Format the coordinate value in a unit appropriate way.

2526 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"

2527

2528 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs

2529 if nplot == 1 and seq_coord.has_bounds: 2529 ↛ 2534line 2529 didn't jump to line 2534 because the condition on line 2529 was always true

2530 if np.size(seq_coord.bounds) > 1: 2530 ↛ 2531line 2530 didn't jump to line 2531 because the condition on line 2530 was never true

2531 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"

2532

2533 # Do the actual plotting.

2534 plotting_func(

2535 cube_slice,

2536 coords,

2537 stamp_coordinate,

2538 plot_filename,

2539 title,

2540 series_coordinate,

2541 )

2542

2543 plot_index.append(plot_filename)

2544 else:

2545 # Format the title and filename using plotted series coordinate

2546 nplot = 1

2547 seq_coord = coords[0]

2548 plot_title, plot_filename = _set_title_and_filename(

2549 seq_coord, nplot, recipe_title, filename

2550 )

2551 # Do the actual plotting for all other series coordinate options.

2552 _plot_and_save_line_series(

2553 cubes, coords, stamp_coordinate, plot_filename, plot_title

2554 )

2555

2556 plot_index.append(plot_filename)

2557

2558 # append plot to list of plots

2559 complete_plot_index = _append_to_plot_index(plot_index)

2560

2561 # Make a page to display the plots.

2562 _make_plot_html_page(complete_plot_index)

2563

2564 return cube

2565

2566

2567def plot_vertical_line_series(

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

2569 filename: str = None,

2570 series_coordinate: str = "model_level_number",

2571 sequence_coordinate: str = "time",

2572 # line_coordinate: str = "realization",

2573 **kwargs,

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

2575 """Plot a line plot against a type of vertical coordinate.

2576

2577 The Cube or CubeList must be 1D.

2578

2579 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present

2580 then a sequence of plots will be produced.

2581

2582 Parameters

2583 ----------

2584 iris.cube | iris.cube.CubeList

2585 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.

2586 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.

2587 We do not support different data such as temperature and humidity in the same CubeList for plotting.

2588 filename: str, optional

2589 Name of the plot to write, used as a prefix for plot sequences. Defaults

2590 to the recipe name.

2591 series_coordinate: str, optional

2592 Coordinate to plot on the y-axis. Can be ``pressure`` or

2593 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``

2594 for LFRic. Defaults to ``model_level_number``.

2595 This coordinate must exist in the cube.

2596 sequence_coordinate: str, optional

2597 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

2598 This coordinate must exist in the cube.

2599

2600 Returns

2601 -------

2602 iris.cube.Cube | iris.cube.CubeList

2603 The original Cube or CubeList (so further operations can be applied).

2604 Plotted data.

2605

2606 Raises

2607 ------

2608 ValueError

2609 If the cubes doesn't have the right dimensions.

2610 TypeError

2611 If the cube isn't a Cube or CubeList.

2612 """

2613 # Ensure we have a name for the plot file.

2614 recipe_title = get_recipe_metadata().get("title", "Untitled")

2615

2616 cubes = iter_maybe(cubes)

2617 # Initialise empty list to hold all data from all cubes in a CubeList

2618 all_data = []

2619

2620 # Store min/max ranges for x range.

2621 x_levels = []

2622

2623 num_models = _get_num_models(cubes)

2624

2625 _validate_cube_shape(cubes, num_models)

2626

2627 # Iterate over all cubes in cube or CubeList and plot.

2628 coords = []

2629 for cube in cubes:

2630 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.

2631 try:

2632 coords.append(cube.coord(series_coordinate))

2633 except iris.exceptions.CoordinateNotFoundError as err:

2634 raise ValueError(

2635 f"Cube must have a {series_coordinate} coordinate."

2636 ) from err

2637

2638 try:

2639 if cube.ndim > 1 or not cube.coords("realization"): 2639 ↛ 2647line 2639 didn't jump to line 2647 because the condition on line 2639 was always true

2640 cube.coord(sequence_coordinate)

2641 except iris.exceptions.CoordinateNotFoundError as err:

2642 raise ValueError(

2643 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."

2644 ) from err

2645

2646 # Get minimum and maximum from levels information.

2647 _, levels, _ = _colorbar_map_levels(cube, axis="x")

2648 if levels is not None: 2648 ↛ 2652line 2648 didn't jump to line 2652 because the condition on line 2648 was always true

2649 x_levels.append(min(levels))

2650 x_levels.append(max(levels))

2651 else:

2652 all_data.append(cube.data)

2653

2654 if len(x_levels) == 0: 2654 ↛ 2656line 2654 didn't jump to line 2656 because the condition on line 2654 was never true

2655 # Combine all data into a single NumPy array

2656 combined_data = np.concatenate(all_data)

2657

2658 # Set the lower and upper limit for the x-axis to ensure all plots have

2659 # same range. This needs to read the whole cube over the range of the

2660 # sequence and if applicable postage stamp coordinate.

2661 vmin = np.floor(combined_data.min())

2662 vmax = np.ceil(combined_data.max())

2663 else:

2664 vmin = min(x_levels)

2665 vmax = max(x_levels)

2666

2667 # Matching the slices (matching by seq coord point; it may happen that

2668 # evaluated models do not cover the same seq coord range, hence matching

2669 # necessary)

2670 def filter_cube_iterables(cube_iterables) -> bool:

2671 return len(cube_iterables) == len(coords)

2672

2673 cube_iterables = filter(

2674 filter_cube_iterables,

2675 (

2676 iris.cube.CubeList(

2677 s

2678 for s in itertools.chain.from_iterable(

2679 cb.slices_over(sequence_coordinate) for cb in cubes

2680 )

2681 if s.coord(sequence_coordinate).points[0] == point

2682 )

2683 for point in sorted(

2684 set(

2685 itertools.chain.from_iterable(

2686 cb.coord(sequence_coordinate).points for cb in cubes

2687 )

2688 )

2689 )

2690 ),

2691 )

2692

2693 # Create a plot for each value of the sequence coordinate.

2694 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for

2695 # similar sequence values. Passing a CubeList into the internal plotting function

2696 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube

2697 # or an iris.cube.CubeList.

2698 plot_index = []

2699 nplot = np.size(cubes[0].coord(sequence_coordinate).points)

2700 for cubes_slice in cube_iterables:

2701 # Format the coordinate value in a unit appropriate way.

2702 seq_coord = cubes_slice[0].coord(sequence_coordinate)

2703 plot_title, plot_filename = _set_title_and_filename(

2704 seq_coord, nplot, recipe_title, filename

2705 )

2706

2707 # Do the actual plotting.

2708 _plot_and_save_vertical_line_series(

2709 cubes_slice,

2710 coords,

2711 "realization",

2712 plot_filename,

2713 series_coordinate,

2714 title=plot_title,

2715 vmin=vmin,

2716 vmax=vmax,

2717 )

2718 plot_index.append(plot_filename)

2719

2720 # Add list of plots to plot metadata.

2721 complete_plot_index = _append_to_plot_index(plot_index)

2722

2723 # Make a page to display the plots.

2724 _make_plot_html_page(complete_plot_index)

2725

2726 return cubes

2727

2728

2729def qq_plot(

2730 cubes: iris.cube.CubeList,

2731 coordinates: list[str],

2732 percentiles: list[float],

2733 model_names: list[str],

2734 filename: str = None,

2735 one_to_one: bool = True,

2736 **kwargs,

2737) -> iris.cube.CubeList:

2738 """Plot a Quantile-Quantile plot between two models for common time points.

2739

2740 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are

2741 collapsed within the operator over all specified coordinates such as

2742 grid_latitude, grid_longitude, vertical levels, but also realisation representing

2743 ensemble members to ensure a 1D cube (array).

2744

2745 Parameters

2746 ----------

2747 cubes: iris.cube.CubeList

2748 Two cubes of the same variable with different models.

2749 coordinate: list[str]

2750 The list of coordinates to collapse over. This list should be

2751 every coordinate within the cube to result in a 1D cube around

2752 the percentile coordinate.

2753 percent: list[float]

2754 A list of percentiles to appear in the plot.

2755 model_names: list[str]

2756 A list of model names to appear on the axis of the plot.

2757 filename: str, optional

2758 Filename of the plot to write.

2759 one_to_one: bool, optional

2760 If True a 1:1 line is plotted; if False it is not. Default is True.

2761

2762 Raises

2763 ------

2764 ValueError

2765 When the cubes are not compatible.

2766

2767 Notes

2768 -----

2769 The quantile-quantile plot is a variant on the scatter plot representing

2770 two datasets by their quantiles (percentiles) for common time points.

2771 This plot does not use a theoretical distribution to compare against, but

2772 compares percentiles of two datasets. This plot does

2773 not use all raw data points, but plots the selected percentiles (quantiles) of

2774 each variable instead for the two datasets, thereby normalising the data for a

2775 direct comparison between the selected percentiles of the two dataset distributions.

2776

2777 Quantile-quantile plots are valuable for comparing against

2778 observations and other models. Identical percentiles between the variables

2779 will lie on the one-to-one line implying the values correspond well to each

2780 other. Where there is a deviation from the one-to-one line a range of

2781 possibilities exist depending on how and where the data is shifted (e.g.,

2782 Wilks 2011 [Wilks2011]_).

2783

2784 For distributions above the one-to-one line the distribution is left-skewed;

2785 below is right-skewed. A distinct break implies a bimodal distribution, and

2786 closer values/values further apart at the tails imply poor representation of

2787 the extremes.

2788

2789 References

2790 ----------

2791 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric

2792 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.

2793 """

2794 # Check cubes using same functionality as the difference operator.

2795 if len(cubes) != 2:

2796 raise ValueError("cubes should contain exactly 2 cubes.")

2797 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))

2798 other: Cube = cubes.extract_cube(

2799 iris.Constraint(

2800 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes

2801 )

2802 )

2803

2804 # Get spatial coord names.

2805 base_lat_name, base_lon_name = get_cube_yxcoordname(base)

2806 other_lat_name, other_lon_name = get_cube_yxcoordname(other)

2807

2808 # Ensure cubes to compare are on common differencing grid.

2809 # This is triggered if either

2810 # i) latitude and longitude shapes are not the same. Note grid points

2811 # are not compared directly as these can differ through rounding

2812 # errors.

2813 # ii) or variables are known to often sit on different grid staggering

2814 # in different models (e.g. cell center vs cell edge), as is the case

2815 # for UM and LFRic comparisons.

2816 # In future greater choice of regridding method might be applied depending

2817 # on variable type. Linear regridding can in general be appropriate for smooth

2818 # variables. Care should be taken with interpretation of differences

2819 # given this dependency on regridding.

2820 if (

2821 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape

2822 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape

2823 ) or (

2824 base.long_name

2825 in [

2826 "eastward_wind_at_10m",

2827 "northward_wind_at_10m",

2828 "northward_wind_at_cell_centres",

2829 "eastward_wind_at_cell_centres",

2830 "zonal_wind_at_pressure_levels",

2831 "meridional_wind_at_pressure_levels",

2832 "potential_vorticity_at_pressure_levels",

2833 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",

2834 ]

2835 ):

2836 logging.debug(

2837 "Linear regridding base cube to other grid to compute differences"

2838 )

2839 base = regrid_onto_cube(base, other, method="Linear")

2840

2841 # Extract just common time points.

2842 base, other = _extract_common_time_points(base, other)

2843

2844 # Equalise attributes so we can merge.

2845 fully_equalise_attributes([base, other])

2846 logging.debug("Base: %s\nOther: %s", base, other)

2847

2848 # Collapse cubes.

2849 base = collapse(

2850 base,

2851 coordinate=coordinates,

2852 method="PERCENTILE",

2853 additional_percent=percentiles,

2854 )

2855 other = collapse(

2856 other,

2857 coordinate=coordinates,

2858 method="PERCENTILE",

2859 additional_percent=percentiles,

2860 )

2861

2862 # Ensure we have a name for the plot file.

2863 recipe_title = get_recipe_metadata().get("title", "Untitled")

2864 title = f"{recipe_title}"

2865

2866 if filename is None:

2867 filename = slugify(recipe_title)

2868

2869 # Add file extension.

2870 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"

2871

2872 # Do the actual plotting on a scatter plot

2873 _plot_and_save_scatter_plot(

2874 base, other, plot_filename, title, one_to_one, model_names

2875 )

2876

2877 # Add list of plots to plot metadata.

2878 plot_index = _append_to_plot_index([plot_filename])

2879

2880 # Make a page to display the plots.

2881 _make_plot_html_page(plot_index)

2882

2883 return iris.cube.CubeList([base, other])

2884

2885

2886def scatter_plot(

2887 cube_x: iris.cube.Cube | iris.cube.CubeList,

2888 cube_y: iris.cube.Cube | iris.cube.CubeList,

2889 filename: str = None,

2890 one_to_one: bool = True,

2891 **kwargs,

2892) -> iris.cube.CubeList:

2893 """Plot a scatter plot between two variables.

2894

2895 Both cubes must be 1D.

2896

2897 Parameters

2898 ----------

2899 cube_x: Cube | CubeList

2900 1 dimensional Cube of the data to plot on y-axis.

2901 cube_y: Cube | CubeList

2902 1 dimensional Cube of the data to plot on x-axis.

2903 filename: str, optional

2904 Filename of the plot to write.

2905 one_to_one: bool, optional

2906 If True a 1:1 line is plotted; if False it is not. Default is True.

2907

2908 Returns

2909 -------

2910 cubes: CubeList

2911 CubeList of the original x and y cubes for further processing.

2912

2913 Raises

2914 ------

2915 ValueError

2916 If the cube doesn't have the right dimensions and cubes not the same

2917 size.

2918 TypeError

2919 If the cube isn't a single cube.

2920

2921 Notes

2922 -----

2923 Scatter plots are used for determining if there is a relationship between

2924 two variables. Positive relations have a slope going from bottom left to top

2925 right; Negative relations have a slope going from top left to bottom right.

2926 """

2927 # Iterate over all cubes in cube or CubeList and plot.

2928 for cube_iter in iter_maybe(cube_x):

2929 # Check cubes are correct shape.

2930 cube_iter = _check_single_cube(cube_iter)

2931 if cube_iter.ndim > 1:

2932 raise ValueError("cube_x must be 1D.")

2933

2934 # Iterate over all cubes in cube or CubeList and plot.

2935 for cube_iter in iter_maybe(cube_y):

2936 # Check cubes are correct shape.

2937 cube_iter = _check_single_cube(cube_iter)

2938 if cube_iter.ndim > 1:

2939 raise ValueError("cube_y must be 1D.")

2940

2941 # Ensure we have a name for the plot file.

2942 recipe_title = get_recipe_metadata().get("title", "Untitled")

2943 title = f"{recipe_title}"

2944

2945 if filename is None:

2946 filename = slugify(recipe_title)

2947

2948 # Add file extension.

2949 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"

2950

2951 # Do the actual plotting.

2952 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)

2953

2954 # Add list of plots to plot metadata.

2955 plot_index = _append_to_plot_index([plot_filename])

2956

2957 # Make a page to display the plots.

2958 _make_plot_html_page(plot_index)

2959

2960 return iris.cube.CubeList([cube_x, cube_y])

2961

2962

2963def vector_plot(

2964 cube_u: iris.cube.Cube,

2965 cube_v: iris.cube.Cube,

2966 filename: str = None,

2967 sequence_coordinate: str = "time",

2968 **kwargs,

2969) -> iris.cube.CubeList:

2970 """Plot a vector plot based on the input u and v components."""

2971 recipe_title = get_recipe_metadata().get("title", "Untitled")

2972

2973 # Cubes must have a matching sequence coordinate.

2974 try:

2975 # Check that the u and v cubes have the same sequence coordinate.

2976 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2976 ↛ anywhereline 2976 didn't jump anywhere: it always raised an exception.

2977 raise ValueError("Coordinates do not match.")

2978 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:

2979 raise ValueError(

2980 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"

2981 ) from err

2982

2983 # Create a plot for each value of the sequence coordinate.

2984 plot_index = []

2985 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)

2986 for cube_u_slice, cube_v_slice in zip(

2987 cube_u.slices_over(sequence_coordinate),

2988 cube_v.slices_over(sequence_coordinate),

2989 strict=True,

2990 ):

2991 # Format the coordinate value in a unit appropriate way.

2992 seq_coord = cube_u_slice.coord(sequence_coordinate)

2993 plot_title, plot_filename = _set_title_and_filename(

2994 seq_coord, nplot, recipe_title, filename

2995 )

2996

2997 # Do the actual plotting.

2998 _plot_and_save_vector_plot(

2999 cube_u_slice,

3000 cube_v_slice,

3001 filename=plot_filename,

3002 title=plot_title,

3003 method="contourf",

3004 )

3005 plot_index.append(plot_filename)

3006

3007 # Add list of plots to plot metadata.

3008 complete_plot_index = _append_to_plot_index(plot_index)

3009

3010 # Make a page to display the plots.

3011 _make_plot_html_page(complete_plot_index)

3012

3013 return iris.cube.CubeList([cube_u, cube_v])

3014

3015

3016def plot_histogram_series(

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

3018 filename: str = None,

3019 sequence_coordinate: str = "time",

3020 stamp_coordinate: str = "realization",

3021 single_plot: bool = False,

3022 **kwargs,

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

3024 """Plot a histogram plot for each vertical level provided.

3025

3026 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)

3027 is present then a sequence of plots will be produced using the time slider

3028 functionality to scroll through histograms against time. If a

3029 stamp_coordinate is present then postage stamp plots will be produced. If

3030 stamp_coordinate and single_plot is True, all postage stamp plots will be

3031 plotted in a single plot instead of separate postage stamp plots.

3032

3033 Parameters

3034 ----------

3035 cubes: Cube | iris.cube.CubeList

3036 Iris cube or CubeList of the data to plot. It should have a single dimension other

3037 than the stamp coordinate.

3038 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.

3039 We do not support different data such as temperature and humidity in the same CubeList for plotting.

3040 filename: str, optional

3041 Name of the plot to write, used as a prefix for plot sequences. Defaults

3042 to the recipe name.

3043 sequence_coordinate: str, optional

3044 Coordinate about which to make a plot sequence. Defaults to ``"time"``.

3045 This coordinate must exist in the cube and will be used for the time

3046 slider.

3047 stamp_coordinate: str, optional

3048 Coordinate about which to plot postage stamp plots. Defaults to

3049 ``"realization"``.

3050 single_plot: bool, optional

3051 If True, all postage stamp plots will be plotted in a single plot. If

3052 False, each postage stamp plot will be plotted separately. Is only valid

3053 if stamp_coordinate exists and has more than a single point.

3054

3055 Returns

3056 -------

3057 iris.cube.Cube | iris.cube.CubeList

3058 The original Cube or CubeList (so further operations can be applied).

3059 Plotted data.

3060

3061 Raises

3062 ------

3063 ValueError

3064 If the cube doesn't have the right dimensions.

3065 TypeError

3066 If the cube isn't a Cube or CubeList.

3067 """

3068 recipe_title = get_recipe_metadata().get("title", "Untitled")

3069

3070 cubes = iter_maybe(cubes)

3071 # Ensure we have a name for the plot file.

3072 if filename is None:

3073 filename = slugify(recipe_title)

3074

3075 # Internal plotting function.

3076 plotting_func = _plot_and_save_histogram_series

3077

3078 num_models = _get_num_models(cubes)

3079

3080 _validate_cube_shape(cubes, num_models)

3081

3082 # If several histograms are plotted with time as sequence_coordinate for the

3083 # time slider option.

3084 for cube in cubes:

3085 try:

3086 cube.coord(sequence_coordinate)

3087 except iris.exceptions.CoordinateNotFoundError as err:

3088 raise ValueError(

3089 f"Cube must have a {sequence_coordinate} coordinate."

3090 ) from err

3091

3092 # Get minimum and maximum from levels information.

3093 levels = None

3094 for cube in cubes: 3094 ↛ 3110line 3094 didn't jump to line 3110 because the loop on line 3094 didn't complete

3095 # First check if user-specified "auto" range variable.

3096 # This maintains the value of levels as None, so proceed.

3097 _, levels, _ = _colorbar_map_levels(cube, axis="y")

3098 if levels is None:

3099 break

3100 # If levels is changed, recheck to use the vmin,vmax or

3101 # levels-based ranges for histogram plots.

3102 _, levels, _ = _colorbar_map_levels(cube)

3103 logging.debug("levels: %s", levels)

3104 if levels is not None: 3104 ↛ 3094line 3104 didn't jump to line 3094 because the condition on line 3104 was always true

3105 vmin = min(levels)

3106 vmax = max(levels)

3107 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)

3108 break

3109

3110 if levels is None:

3111 vmin = min(cb.data.min() for cb in cubes)

3112 vmax = max(cb.data.max() for cb in cubes)

3113

3114 # Make postage stamp plots if stamp_coordinate exists and has more than a

3115 # single point. If single_plot is True:

3116 # -- all postage stamp plots will be plotted in a single plot instead of

3117 # separate postage stamp plots.

3118 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are

3119 # produced per single model only

3120

3121 if num_models == 1: 3121 ↛ 3134line 3121 didn't jump to line 3134 because the condition on line 3121 was always true

3122 if ( 3122 ↛ 3126line 3122 didn't jump to line 3126 because the condition on line 3122 was never true

3123 stamp_coordinate in [c.name() for c in cubes[0].coords()]

3124 and cubes[0].coord(stamp_coordinate).shape[0] > 1

3125 ):

3126 if single_plot:

3127 plotting_func = (

3128 _plot_and_save_postage_stamps_in_single_plot_histogram_series

3129 )

3130 else:

3131 plotting_func = _plot_and_save_postage_stamp_histogram_series

3132 cube_iterables = cubes[0].slices_over(sequence_coordinate)

3133 else:

3134 all_points = sorted(

3135 set(

3136 itertools.chain.from_iterable(

3137 cb.coord(sequence_coordinate).points for cb in cubes

3138 )

3139 )

3140 )

3141 all_slices = list(

3142 itertools.chain.from_iterable(

3143 cb.slices_over(sequence_coordinate) for cb in cubes

3144 )

3145 )

3146 # Matched slices (matched by seq coord point; it may happen that

3147 # evaluated models do not cover the same seq coord range, hence matching

3148 # necessary)

3149 cube_iterables = [

3150 iris.cube.CubeList(

3151 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point

3152 )

3153 for point in all_points

3154 ]

3155

3156 plot_index = []

3157 nplot = np.size(cube.coord(sequence_coordinate).points)

3158 # Create a plot for each value of the sequence coordinate. Allowing for

3159 # multiple cubes in a CubeList to be plotted in the same plot for similar

3160 # sequence values. Passing a CubeList into the internal plotting function

3161 # for similar values of the sequence coordinate. cube_slice can be an

3162 # iris.cube.Cube or an iris.cube.CubeList.

3163 for cube_slice in cube_iterables:

3164 single_cube = cube_slice

3165 if isinstance(cube_slice, iris.cube.CubeList): 3165 ↛ 3166line 3165 didn't jump to line 3166 because the condition on line 3165 was never true

3166 single_cube = cube_slice[0]

3167

3168 # Set plot titles and filename, based on sequence coordinate

3169 seq_coord = single_cube.coord(sequence_coordinate)

3170 # Use time coordinate in title and filename if single histogram output.

3171 if sequence_coordinate == "realization" and nplot == 1: 3171 ↛ 3172line 3171 didn't jump to line 3172 because the condition on line 3171 was never true

3172 seq_coord = single_cube.coord("time")

3173 plot_title, plot_filename = _set_title_and_filename(

3174 seq_coord, nplot, recipe_title, filename

3175 )

3176

3177 # Do the actual plotting.

3178

3179 plotting_func(

3180 cube_slice,

3181 filename=plot_filename,

3182 stamp_coordinate=stamp_coordinate,

3183 title=plot_title,

3184 vmin=vmin,

3185 vmax=vmax,

3186 )

3187 plot_index.append(plot_filename)

3188

3189 # Add list of plots to plot metadata.

3190 complete_plot_index = _append_to_plot_index(plot_index)

3191

3192 # Make a page to display the plots.

3193 _make_plot_html_page(complete_plot_index)

3194

3195 return cubes

3196

3197

3198def _plot_and_save_postage_stamp_power_spectrum_series(

3199 cubes: iris.cube.Cube,

3200 coords: list[iris.coords.Coord],

3201 stamp_coordinate: str,

3202 filename: str,

3203 title: str,

3204 series_coordinate: str = None,

3205 **kwargs,

3206):

3207 """Plot and save postage (ensemble members) stamps for a power spectrum series.

3208

3209 Parameters

3210 ----------

3211 cubes: Cube or CubeList

3212 Cube or Cubelist of the power spectrum data.

3213 coords: list[Coord]

3214 Coordinates to plot on the x-axis, one per cube.

3215 stamp_coordinate: str

3216 Coordinate that becomes different plots.

3217 filename: str

3218 Filename of the plot to write.

3219 title: str

3220 Plot title.

3221 series_coordinate: str, optional

3222 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.

3223

3224 """

3225 # Use the smallest square grid that will fit the members.

3226 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))

3227

3228 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")

3229 model_colors_map = _get_model_colors_map(cubes)

3230 # ax = plt.gca()

3231 # Make a subplot for each member.

3232 for member, subplot in zip(

3233 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False

3234 ):

3235 ax = plt.subplot(grid_size, grid_size, subplot)

3236

3237 # Store min/max ranges.

3238 y_levels = []

3239

3240 line_marker = None

3241 line_width = 1

3242

3243 for cube in iter_maybe(member):

3244 xcoord = select_series_coord(cube, series_coordinate)

3245 xname = xcoord.points

3246

3247 yfield = cube.data # power spectrum

3248 label = None

3249 color = "black"

3250 if model_colors_map: 3250 ↛ 3251line 3250 didn't jump to line 3251 because the condition on line 3250 was never true

3251 label = cube.attributes.get("model_name")

3252 color = model_colors_map.get(label)

3253

3254 if member.coord(stamp_coordinate).points == [0]:

3255 ax.plot(

3256 xname,

3257 yfield,

3258 color=color,

3259 marker=line_marker,

3260 ls="-",

3261 lw=line_width,

3262 label=f"{label} (control)"

3263 if len(cube.coord(stamp_coordinate).points) > 1

3264 else label,

3265 )

3266 # Label with member if part of an ensemble and not the control.

3267 else:

3268 ax.plot(

3269 xname,

3270 yfield,

3271 color=color,

3272 ls="-",

3273 lw=1.5,

3274 alpha=0.75,

3275 label=f"{label} (member)",

3276 )

3277

3278 # Calculate the global min/max if multiple cubes are given.

3279 _, levels, _ = _colorbar_map_levels(cube, axis="y")

3280 if levels is not None: 3280 ↛ 3281line 3280 didn't jump to line 3281 because the condition on line 3280 was never true

3281 y_levels.append(min(levels))

3282 y_levels.append(max(levels))

3283

3284 # Add some labels and tweak the style.

3285 title = f"{title}"

3286 ax.set_title(title, fontsize=16)

3287

3288 # Set appropriate x-axis label based on coordinate

3289 if series_coordinate == "wavelength" or ( 3289 ↛ 3292line 3289 didn't jump to line 3292 because the condition on line 3289 was never true

3290 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"

3291 ):

3292 ax.set_xlabel("Wavelength (km)", fontsize=14)

3293 elif series_coordinate == "physical_wavenumber" or ( 3293 ↛ 3298line 3293 didn't jump to line 3298 because the condition on line 3293 was always true

3294 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"

3295 ):

3296 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

3297 else: # frequency or check units

3298 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":

3299 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

3300 else:

3301 ax.set_xlabel("Wavenumber", fontsize=14)

3302

3303 ax.set_ylabel("Power Spectral Density", fontsize=14)

3304 ax.tick_params(axis="both", labelsize=12)

3305

3306 # Set log-log scale

3307 ax.set_xscale("log")

3308 ax.set_yscale("log")

3309

3310 # Add gridlines

3311 ax.grid(linestyle="--", color="grey", linewidth=1)

3312 # Ientify unique labels for legend

3313 handles = list(

3314 {

3315 label: handle

3316 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)

3317 }.values()

3318 )

3319 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)

3320

3321 ax = plt.gca()

3322 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")

3323

3324 # Overall figure title.

3325 fig.suptitle(title, fontsize=16)

3326

3327 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

3328 logging.info("Saved histogram postage stamp plot to %s", filename)

3329 plt.close(fig)

3330

3331

3332def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(

3333 cubes: iris.cube.Cube,

3334 coords: list[iris.coords.Coord],

3335 stamp_coordinate: str,

3336 filename: str,

3337 title: str,

3338 series_coordinate: str = None,

3339 **kwargs,

3340):

3341 """Plot and save power spectra for ensemble members in single plot.

3342

3343 Parameters

3344 ----------

3345 cubes: Cube or CubeList

3346 Cube or Cubelist of the power spectrum data.

3347 coords: list[Coord]

3348 Coordinates to plot on the x-axis, one per cube.

3349 stamp_coordinate: str

3350 Coordinate that becomes different plots.

3351 filename: str

3352 Filename of the plot to write.

3353 title: str

3354 Plot title.

3355 series_coordinate: str, optional

3356 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.

3357

3358 """

3359 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")

3360 model_colors_map = _get_model_colors_map(cubes)

3361

3362 line_marker = None

3363 line_width = 1

3364

3365 # Compute ensemble statistics to show spread

3366 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)

3367 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)

3368 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)

3369

3370 xcoord_global = mean_cube.coord(series_coordinate)

3371 x_global = xcoord_global.points

3372

3373 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):

3374 xcoord = select_series_coord(member, series_coordinate)

3375 xname = xcoord.points

3376

3377 yfield = member.data # power spectrum

3378 color = "black"

3379 if model_colors_map: 3379 ↛ 3383line 3379 didn't jump to line 3383 because the condition on line 3379 was always true

3380 label = member.attributes.get("model_name") if i == 0 else None

3381 color = model_colors_map.get(label)

3382

3383 if member.coord(stamp_coordinate).points == [0]:

3384 ax.plot(

3385 xname,

3386 yfield,

3387 color=color,

3388 marker=line_marker,

3389 ls="-",

3390 lw=line_width,

3391 label=f"{label} (control)"

3392 if len(member.coord(stamp_coordinate).points) > 1

3393 else label,

3394 )

3395 # Label with member number if part of an ensemble and not the control.

3396 else:

3397 ax.plot(

3398 xname,

3399 yfield,

3400 color=color,

3401 ls="-",

3402 lw=1.5,

3403 alpha=0.75,

3404 label=label,

3405 )

3406

3407 # Set appropriate x-axis label based on coordinate

3408 if series_coordinate == "wavelength" or ( 3408 ↛ 3411line 3408 didn't jump to line 3411 because the condition on line 3408 was never true

3409 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"

3410 ):

3411 ax.set_xlabel("Wavelength (km)", fontsize=14)

3412 elif series_coordinate == "physical_wavenumber" or ( 3412 ↛ 3417line 3412 didn't jump to line 3417 because the condition on line 3412 was always true

3413 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"

3414 ):

3415 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

3416 else: # frequency or check units

3417 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":

3418 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)

3419 else:

3420 ax.set_xlabel("Wavenumber", fontsize=14)

3421

3422 # Add ensemble spread shading

3423 ax.fill_between(

3424 x_global,

3425 min_cube.data,

3426 max_cube.data,

3427 color="grey",

3428 alpha=0.3,

3429 label="Ensemble spread",

3430 )

3431

3432 # Add ensemble mean line

3433 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")

3434

3435 ax.set_ylabel("Power Spectral Density", fontsize=14)

3436 ax.tick_params(axis="both", labelsize=12)

3437

3438 # Set y limits to global min and max, autoscale if colorbar doesn't exist.

3439 # Set log-log scale

3440 ax.set_xscale("log")

3441 ax.set_yscale("log")

3442

3443 # Add gridlines

3444 ax.grid(linestyle="--", color="grey", linewidth=1)

3445 # Identify unique labels for legend

3446 handles = list(

3447 {

3448 label: handle

3449 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)

3450 }.values()

3451 )

3452 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)

3453

3454 # Add a legend

3455 ax.legend(fontsize=16)

3456

3457 # Figure title.

3458 ax.set_title(title, fontsize=16)

3459

3460 # Save the figure to a file

3461 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())

3462

3463 # Close the figure

3464 plt.close(fig)

Coverage for src / CSET / operators / plot.py: 81%

1072 statements