Coverage for src/CSET/operators/plot.py: 80%
1025 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-20 14:29 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24from typing import Literal
26import cartopy.crs as ccrs
27import cartopy.feature as cfeature
28import iris
29import iris.coords
30import iris.cube
31import iris.exceptions
32import iris.plot as iplt
33import matplotlib as mpl
34import matplotlib.pyplot as plt
35import numpy as np
36from cartopy.mpl.geoaxes import GeoAxes
37from iris.cube import Cube
38from markdown_it import MarkdownIt
39from mpl_toolkits.axes_grid1.inset_locator import inset_axes
41from CSET._common import (
42 filename_slugify,
43 get_recipe_metadata,
44 iter_maybe,
45 render_file,
46 slugify,
47)
48from CSET.operators._colormaps import (
49 colorbar_map_levels,
50 get_model_colors_map,
51)
52from CSET.operators._utils import (
53 check_sequence_coordinate,
54 check_single_cube,
55 check_stamp_coordinate,
56 fully_equalise_attributes,
57 get_cube_yxcoordname,
58 get_num_models,
59 is_transect,
60 slice_over_maybe,
61 validate_cube_shape,
62 validate_cubes_coords,
63)
64from CSET.operators.collapse import collapse
65from CSET.operators.misc import _extract_common_time_points
66from CSET.operators.regrid import regrid_onto_cube
68# Use a non-interactive plotting backend.
69mpl.use("agg")
72############################
73# Private helper functions #
74############################
77def _append_to_plot_index(plot_index: list) -> list:
78 """Add plots into the plot index, returning the complete plot index."""
79 with open("meta.json", "r+t", encoding="UTF-8") as fp:
80 fcntl.flock(fp, fcntl.LOCK_EX)
81 fp.seek(0)
82 meta = json.load(fp)
83 complete_plot_index = meta.get("plots", [])
84 complete_plot_index = complete_plot_index + plot_index
85 meta["plots"] = complete_plot_index
86 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
87 os.getenv("DO_CASE_AGGREGATION")
88 ):
89 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
90 fp.seek(0)
91 fp.truncate()
92 json.dump(meta, fp, indent=2)
93 return complete_plot_index
96def _make_plot_html_page(plots: list):
97 """Create a HTML page to display a plot image."""
98 # Debug check that plots actually contains some strings.
99 assert isinstance(plots[0], str)
101 # Load HTML template file.
102 operator_files = importlib.resources.files()
103 template_file = operator_files.joinpath("_plot_page_template.html")
105 # Get some metadata.
106 meta = get_recipe_metadata()
107 title = meta.get("title", "Untitled")
108 description = MarkdownIt().render(meta.get("description", "*No description.*"))
110 # Prepare template variables.
111 variables = {
112 "title": title,
113 "description": description,
114 "initial_plot": plots[0],
115 "plots": plots,
116 "title_slug": slugify(title),
117 }
119 # Render template.
120 html = render_file(template_file, **variables)
122 # Save completed HTML.
123 with open("index.html", "wt", encoding="UTF-8") as fp:
124 fp.write(html)
127def _setup_spatial_map(
128 cube: iris.cube.Cube,
129 figure,
130 cmap,
131 grid_size: tuple[int, int] | None = None,
132 subplot: int | None = None,
133):
134 """Define map projections, extent and add coastlines and borderlines for spatial plots.
136 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
137 is specified, and map extent defined based on the input data.
139 Parameters
140 ----------
141 cube: Cube
142 2 dimensional (lat and lon) Cube of the data to plot.
143 figure:
144 Matplotlib Figure object holding all plot elements.
145 cmap:
146 Matplotlib colormap.
147 grid_size: (int, int), optional
148 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
149 subplot: int, optional
150 Subplot index if multiple spatial subplots in figure.
152 Returns
153 -------
154 axes:
155 Matplotlib GeoAxes definition.
156 """
157 # Identify min/max plot bounds.
158 try:
159 lat_axis, lon_axis = get_cube_yxcoordname(cube)
160 xmin = np.nanmin(cube.coord(lon_axis).points)
161 xmax = np.nanmax(cube.coord(lon_axis).points)
162 ymin = np.nanmin(cube.coord(lat_axis).points)
163 ymax = np.nanmax(cube.coord(lat_axis).points)
165 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
166 if np.abs(xmax - xmin) > 180.0:
167 xmin = xmin - 180.0
168 xmax = xmax - 180.0
169 logging.debug("Adjusting plot bounds to fit global extent.")
171 # Consider map projection orientation.
172 # Adapting orientation enables plotting across international dateline.
173 # Users can adapt the default central_longitude if alternative projections views.
174 if xmax > 180.0 or xmin < -180.0:
175 central_longitude = 180.0
176 else:
177 central_longitude = 0.0
179 # Define spatial map projection.
180 coord_system = cube.coord(lat_axis).coord_system
181 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
182 # Define rotated pole map projection for rotated pole inputs.
183 projection = ccrs.RotatedPole(
184 pole_longitude=coord_system.grid_north_pole_longitude,
185 pole_latitude=coord_system.grid_north_pole_latitude,
186 central_rotated_longitude=central_longitude,
187 )
188 crs = projection
189 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 189 ↛ 191line 189 didn't jump to line 191 because the condition on line 189 was never true
190 # Define Transverse Mercator projection for TM inputs.
191 projection = ccrs.TransverseMercator(
192 central_longitude=coord_system.longitude_of_central_meridian,
193 central_latitude=coord_system.latitude_of_projection_origin,
194 false_easting=coord_system.false_easting,
195 false_northing=coord_system.false_northing,
196 scale_factor=coord_system.scale_factor_at_central_meridian,
197 )
198 crs = projection
199 else:
200 # Assume polar projection for regional grids encompassing N. Pole
201 if ymin > 20.0 and ymax > 80.0:
202 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
203 elif ymin < -80.0 and ymax < -20.0:
204 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
205 # Define regular map projection for non-rotated pole inputs.
206 # Alternatives might include e.g. for global model outputs:
207 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
208 # projection = ccrs.NearsidePerspective(
209 # central_longitude=180.0,
210 # central_latitude=0,
211 # satellite_height=35785831,
212 # )
213 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
214 else:
215 projection = ccrs.PlateCarree(central_longitude=central_longitude)
216 crs = ccrs.PlateCarree()
218 # Define axes for plot (or subplot) with required map projection.
219 if subplot is not None:
220 axes = figure.add_subplot(
221 grid_size[0], grid_size[1], subplot, projection=projection
222 )
223 else:
224 axes = figure.add_subplot(projection=projection)
226 # Add coastlines and borderlines if cube contains x and y map coordinates.
227 # Avoid adding lines for 2D masked data or specific fixed ancillary spatial plots.
228 if (cube.ndim > 1 and iris.util.is_masked(cube.data)) or any(
229 name in cube.name() for name in ["land_", "orography", "altitude"]
230 ):
231 pass
232 else:
233 if cmap.name in ["viridis", "Greys"]:
234 coastcol = "magenta"
235 else:
236 coastcol = "black"
237 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
238 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
239 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
241 # Add gridlines.
242 gl = axes.gridlines(
243 alpha=0.3,
244 draw_labels=True,
245 dms=False,
246 x_inline=False,
247 y_inline=False,
248 )
249 gl.top_labels = False
250 gl.right_labels = False
251 if subplot:
252 gl.bottom_labels = False
253 gl.left_labels = False
254 if subplot % grid_size[1] == 1:
255 gl.left_labels = True
256 if subplot > ((grid_size[0] - 1) * grid_size[1]): 256 ↛ 261line 256 didn't jump to line 261 because the condition on line 256 was always true
257 gl.bottom_labels = True
259 # If is lat/lon spatial map, fix extent to keep plot tight.
260 # Specifying crs within set_extent helps ensure only data region is shown.
261 if isinstance(coord_system, iris.coord_systems.GeogCS):
262 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
264 except ValueError:
265 # Skip if not both x and y map coordinates.
266 axes = figure.gca()
267 pass
269 return axes
272def _get_plot_resolution() -> int:
273 """Get resolution of rasterised plots in pixels per inch."""
274 return get_recipe_metadata().get("plot_resolution", 100)
277def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
278 """Return title and filename based on start and end points or bounds."""
279 if use_bounds and seq_coord.has_bounds():
280 vals = seq_coord.bounds.flatten()
281 else:
282 vals = seq_coord.points
283 start = seq_coord.units.title(vals[0])
284 end = seq_coord.units.title(vals[-1])
286 if start == end:
287 sequence_title = f"\n [{start}]"
288 sequence_fname = f"_{filename_slugify(start)}"
289 else:
290 sequence_title = f"\n [{start} to {end}]"
291 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
293 # Do not include time if coord set to zero.
294 if (
295 seq_coord.units == "hours since 0001-01-01 00:00:00"
296 and vals[0] == 0
297 and vals[-1] == 0
298 ):
299 sequence_title = ""
300 sequence_fname = ""
302 return sequence_title, sequence_fname
305def _set_title_and_filename(
306 seq_coord: iris.coords.Coord,
307 nplot: int,
308 recipe_title: str,
309 filename: str,
310):
311 """Set plot title and filename based on cube coordinate.
313 Parameters
314 ----------
315 sequence_coordinate: iris.coords.Coord
316 Coordinate about which to make a plot sequence.
317 nplot: int
318 Number of output plots to generate - controls title/naming.
319 recipe_title: str
320 Default plot title, potentially to update.
321 filename: str
322 Input plot filename, potentially to update.
324 Returns
325 -------
326 plot_title: str
327 Output formatted plot title string, based on plotted data.
328 plot_filename: str
329 Output formatted plot filename string.
330 """
331 ndim = seq_coord.ndim
332 npoints = np.size(seq_coord.points)
333 sequence_title = ""
334 sequence_fname = ""
336 # Case 1: Multiple dimension sequence input - list number of aggregated cases
337 # (e.g. aggregation histogram plots)
338 if ndim > 1:
339 ncase = np.shape(seq_coord)[0]
340 sequence_title = f"\n [{ncase} cases]"
341 sequence_fname = f"_{ncase}cases"
343 # Case 2: Single dimension input
344 else:
345 # Single sequence point
346 if npoints == 1:
347 if nplot > 1:
348 # Default labels for sequence inputs
349 sequence_value = seq_coord.units.title(seq_coord.points[0])
350 sequence_title = f"\n [{sequence_value}]"
351 sequence_fname = f"_{filename_slugify(sequence_value)}"
352 else:
353 # Aggregated attribute available where input collapsed over aggregation
354 try:
355 ncase = seq_coord.attributes["number_reference_times"]
356 sequence_title = f"\n [{ncase} cases]"
357 sequence_fname = f"_{ncase}cases"
358 except KeyError:
359 sequence_title, sequence_fname = _get_start_end_strings(
360 seq_coord, use_bounds=seq_coord.has_bounds()
361 )
362 # Multiple sequence (e.g. time) points
363 else:
364 sequence_title, sequence_fname = _get_start_end_strings(
365 seq_coord, use_bounds=False
366 )
368 # Set plot title and filename
369 plot_title = f"{recipe_title}{sequence_title}"
371 # Set plot filename, defaulting to user input if provided.
372 if filename is None:
373 filename = slugify(recipe_title)
374 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
375 else:
376 if nplot > 1:
377 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
378 else:
379 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
381 return plot_title, plot_filename
384def _select_series_coord(cube, series_coordinate):
385 """Determine the grid coordinates to use to calculate grid spacing."""
386 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
387 if series_coordinate in spacing_coordinates: 387 ↛ 393line 387 didn't jump to line 393 because the condition on line 387 was always true
388 # Try the requested coordinate first then the fallbacks in order.
389 fallbacks = [series_coordinate] + [
390 c for c in spacing_coordinates if c != series_coordinate
391 ]
392 else:
393 fallbacks = {series_coordinate}
395 # Try each possible coordinate.
396 for coord in fallbacks:
397 try:
398 return cube.coord(coord)
399 except iris.exceptions.CoordinateNotFoundError:
400 logging.debug("Coordinate %s not found.", coord)
402 # If we get here, none of the fallback options were found.
403 raise iris.exceptions.CoordinateNotFoundError(
404 f"No valid coordinate found for '{series_coordinate}' "
405 f"or fallback options {fallbacks}"
406 )
409def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
410 """Control postage stamp plot output titles based on stamp coordinate."""
411 if stamp_coord.name() == "realization":
412 mtitle = "Member"
413 else:
414 mtitle = stamp_coord.name().capitalize()
416 if stamp_coord.name() == "time":
417 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
418 else:
419 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
421 return mtitle
424def _set_axis_range(cubes):
425 """Get minimum and maximum from levels information."""
426 levels = None
427 for cube in cubes: 427 ↛ 443line 427 didn't jump to line 443 because the loop on line 427 didn't complete
428 # First check if user-specified "auto" range variable.
429 # This maintains the value of levels as None, so proceed.
430 _, levels, _ = colorbar_map_levels(cube, axis="y")
431 if levels is None:
432 break
433 # If levels is changed, recheck to use the vmin,vmax or
434 # levels-based ranges for histogram plots.
435 _, levels, _ = colorbar_map_levels(cube)
436 logging.debug("levels: %s", levels)
437 if levels is not None: 437 ↛ 427line 437 didn't jump to line 427 because the condition on line 437 was always true
438 vmin = min(levels)
439 vmax = max(levels)
440 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
441 break
443 if levels is None:
444 vmin = min(cb.data.min() for cb in cubes)
445 vmax = max(cb.data.max() for cb in cubes)
447 return vmin, vmax
450def _find_matched_slices(cubes, sequence_coordinate):
451 """Identify matched cubes in CubeList by sequence_coordinate values.
453 Ensures common points are compared for multiple cube inputs.
454 """
455 all_points = sorted(
456 set(
457 itertools.chain.from_iterable(
458 cb.coord(sequence_coordinate).points for cb in cubes
459 )
460 )
461 )
462 all_slices = list(
463 itertools.chain.from_iterable(
464 cb.slices_over(sequence_coordinate) for cb in cubes
465 )
466 )
467 # Matched slices (matched by seq coord point; it may happen that
468 # evaluated models do not cover the same seq coord range, hence matching
469 # necessary)
470 cube_iterables = [
471 iris.cube.CubeList(
472 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
473 )
474 for point in all_points
475 ]
477 return cube_iterables
480def _plot_and_save_spatial_plot(
481 cube: iris.cube.Cube,
482 filename: str,
483 title: str,
484 method: Literal["contourf", "pcolormesh", "scatter"],
485 overlay_cube: iris.cube.Cube | None = None,
486 contour_cube: iris.cube.Cube | None = None,
487 point_cube: iris.cube.Cube | None = None,
488 **kwargs,
489):
490 """Plot and save a spatial plot.
492 Parameters
493 ----------
494 cube: Cube
495 2 dimensional (lat and lon) Cube of the data to plot.
496 filename: str
497 Filename of the plot to write.
498 title: str
499 Plot title.
500 method: "contourf" | "pcolormesh" | "scatter"
501 The plotting method to use
502 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
503 overlay_cube: Cube, optional
504 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
505 contour_cube: Cube, optional
506 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
507 point_cube: Cube, optional
508 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
509 """
510 # Setup plot details, size, resolution, etc.
511 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
513 # Specify the color bar
514 cmap, levels, norm = colorbar_map_levels(cube)
516 # If overplotting, set required colorbars
517 if overlay_cube:
518 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
519 if contour_cube:
520 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
522 # Setup plot map projection, extent and coastlines and borderlines.
523 axes = _setup_spatial_map(cube, fig, cmap)
525 # Set colorscale bounds
526 try:
527 vmin = min(levels)
528 vmax = max(levels)
529 except TypeError:
530 vmin, vmax = None, None
531 # Ensure to use norm and not vmin/vmax if levels are defined.
532 if norm is not None:
533 vmin = None
534 vmax = None
535 logging.debug("Plotting using defined levels.")
537 # Plot the field.
538 if method == "contourf":
539 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
540 elif method == "pcolormesh":
541 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
542 elif method == "scatter": 542 ↛ 550line 542 didn't jump to line 550 because the condition on line 542 was never true
543 # Scatter plot of the field. The marker size is chosen to give
544 # symbols that decrease in size as the number of data points
545 # increases, although the fraction of the figure covered by
546 # symbols increases roughly as N^(1/2), disregarding overlaps,
547 # and has been selected for the default figure size of (10, 10).
548 # Should this be changed, the marker size should be adjusted in
549 # proportion to the area of the figure.
550 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
551 lat_axis, lon_axis = get_cube_yxcoordname(cube)
552 plot = iplt.scatter(
553 cube.coord(lon_axis),
554 cube.coord(lat_axis),
555 c=cube.data[:],
556 s=mrk_size,
557 cmap=cmap,
558 edgecolors="k",
559 norm=norm,
560 vmin=vmin,
561 vmax=vmax,
562 )
563 else:
564 raise ValueError(f"Unknown plotting method: {method}")
566 # Overplot overlay field, if required
567 if overlay_cube:
568 try:
569 over_vmin = min(over_levels)
570 over_vmax = max(over_levels)
571 except TypeError:
572 over_vmin, over_vmax = None, None
573 if over_norm is not None: 573 ↛ 574line 573 didn't jump to line 574 because the condition on line 573 was never true
574 over_vmin = None
575 over_vmax = None
576 overlay = iplt.pcolormesh(
577 overlay_cube,
578 cmap=over_cmap,
579 norm=over_norm,
580 alpha=0.8,
581 vmin=over_vmin,
582 vmax=over_vmax,
583 )
584 # Overplot contour field, if required, with contour labelling.
585 if contour_cube:
586 contour = iplt.contour(
587 contour_cube,
588 colors="darkgray",
589 levels=cntr_levels,
590 norm=cntr_norm,
591 alpha=0.5,
592 linestyles="--",
593 linewidths=1,
594 )
595 plt.clabel(contour)
596 # Overplot valid elements of point-based field, if required.
597 # Check for non-masked points only to avoid plotting missing data.
598 if point_cube: 598 ↛ 599line 598 didn't jump to line 599 because the condition on line 598 was never true
599 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
600 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
601 lon_coord = point_cube.coord(lon_axis)
602 lat_coord = point_cube.coord(lat_axis)
603 valid = ~point_cube.data.mask
604 valid_lon = iris.coords.AuxCoord(
605 lon_coord.points[valid],
606 standard_name=lon_coord.standard_name,
607 units=lon_coord.units,
608 coord_system=lon_coord.coord_system,
609 )
610 valid_lat = iris.coords.AuxCoord(
611 lat_coord.points[valid],
612 standard_name=lat_coord.standard_name,
613 units=lat_coord.units,
614 coord_system=lat_coord.coord_system,
615 )
616 iplt.scatter(
617 valid_lon,
618 valid_lat,
619 c=point_cube.data[valid],
620 s=mrk_size,
621 cmap=cmap,
622 edgecolors="k",
623 norm=norm,
624 vmin=vmin,
625 vmax=vmax,
626 )
628 # Check to see if transect, and if so, adjust y axis.
629 if is_transect(cube):
630 if "pressure" in [coord.name() for coord in cube.coords()]:
631 axes.invert_yaxis()
632 axes.set_yscale("log")
633 axes.set_ylim(1100, 100)
634 # If both model_level_number and level_height exists, iplt can construct
635 # plot as a function of height above orography (NOT sea level).
636 elif {"model_level_number", "level_height"}.issubset( 636 ↛ 641line 636 didn't jump to line 641 because the condition on line 636 was always true
637 {coord.name() for coord in cube.coords()}
638 ):
639 axes.set_yscale("log")
641 axes.set_title(
642 f"{title}\n"
643 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
644 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
645 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
646 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
647 fontsize=16,
648 )
650 # Inset code
651 axins = inset_axes(
652 axes,
653 width="20%",
654 height="20%",
655 loc="upper right",
656 axes_class=GeoAxes,
657 axes_kwargs={"map_projection": ccrs.PlateCarree()},
658 )
660 # Slightly transparent to reduce plot blocking.
661 axins.patch.set_alpha(0.4)
663 axins.coastlines(resolution="50m")
664 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
666 SLat, SLon, ELat, ELon = (
667 float(coord) for coord in cube.attributes["transect_coords"].split("_")
668 )
670 # Draw line between them
671 axins.plot(
672 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
673 )
675 # Plot points (note: lon, lat order for Cartopy)
676 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
677 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
679 lon_min, lon_max = sorted([SLon, ELon])
680 lat_min, lat_max = sorted([SLat, ELat])
682 # Midpoints
683 lon_mid = (lon_min + lon_max) / 2
684 lat_mid = (lat_min + lat_max) / 2
686 # Maximum half-range
687 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
688 if half_range == 0: # points identical → provide small default 688 ↛ 692line 688 didn't jump to line 692 because the condition on line 688 was always true
689 half_range = 1
691 # Set square extent
692 axins.set_extent(
693 [
694 lon_mid - half_range,
695 lon_mid + half_range,
696 lat_mid - half_range,
697 lat_mid + half_range,
698 ],
699 crs=ccrs.PlateCarree(),
700 )
702 # Ensure square aspect
703 axins.set_aspect("equal")
705 else:
706 # Add title.
707 axes.set_title(title, fontsize=16)
709 # Adjust padding if spatial plot or transect
710 if is_transect(cube):
711 yinfopad = -0.1
712 ycbarpad = 0.1
713 else:
714 yinfopad = 0.01
715 ycbarpad = 0.042
717 # Add watermark with min/max/mean. Currently not user togglable.
718 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
719 axes.annotate(
720 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
721 xy=(0.025, yinfopad),
722 xycoords="axes fraction",
723 xytext=(-5, 5),
724 textcoords="offset points",
725 ha="left",
726 va="bottom",
727 size=11,
728 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
729 )
731 # Add secondary colour bar for overlay_cube field if required.
732 if overlay_cube:
733 cbarB = fig.colorbar(
734 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
735 )
736 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
737 # add ticks and tick_labels for every levels if less than 20 levels exist
738 if over_levels is not None and len(over_levels) < 20: 738 ↛ 739line 738 didn't jump to line 739 because the condition on line 738 was never true
739 cbarB.set_ticks(over_levels)
740 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
741 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
742 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
743 logging.debug("Set secondary colorbar ticks and labels.")
745 # Add main colour bar.
746 cbar = fig.colorbar(
747 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
748 )
750 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
751 # add ticks and tick_labels for every levels if less than 20 levels exist
752 if levels is not None and len(levels) < 20:
753 cbar.set_ticks(levels)
754 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
755 if "rainfall" or "snowfall" or "visibility" in cube.name(): 755 ↛ 758line 755 didn't jump to line 758 because the condition on line 755 was always true
756 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
757 # Tick labels for rainfall rates from Nimrod radar data.
758 if "rainfall rate composite" in cube.name(): 758 ↛ 759line 758 didn't jump to line 759 because the condition on line 758 was never true
759 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
760 # Tick labels for rain accumulations from Nimrod radar data.
761 if "rain accumulation" in cube.name(): 761 ↛ 762line 761 didn't jump to line 762 because the condition on line 761 was never true
762 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
763 if "wts accumulation" in cube.name(): 763 ↛ 764line 763 didn't jump to line 764 because the condition on line 763 was never true
764 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
765 cbar.minorticks_off()
766 cbar.set_ticks(tick_levels)
767 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
768 cbar.set_label(label=f"{cube.name()}", size=14)
769 # Tick labels for model rainfall data.
770 if "surface_microphysical" in cube.name(): 770 ↛ 773line 770 didn't jump to line 773 because the condition on line 770 was always true
771 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
772 # Tick labels for Nimrod weights data.
773 logging.debug("Set colorbar ticks and labels.")
775 # Save plot.
776 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
777 logging.info("Saved spatial plot to %s", filename)
778 plt.close(fig)
781def _plot_and_save_postage_stamp_spatial_plot(
782 cube: iris.cube.Cube,
783 filename: str,
784 stamp_coordinate: str,
785 title: str,
786 method: Literal["contourf", "pcolormesh"],
787 overlay_cube: iris.cube.Cube | None = None,
788 contour_cube: iris.cube.Cube | None = None,
789 **kwargs,
790):
791 """Plot postage stamp spatial plots from an ensemble.
793 Parameters
794 ----------
795 cube: Cube
796 Iris cube of data to be plotted. It must have the stamp coordinate.
797 filename: str
798 Filename of the plot to write.
799 stamp_coordinate: str
800 Coordinate that becomes different plots.
801 method: "contourf" | "pcolormesh"
802 The plotting method to use.
803 overlay_cube: Cube, optional
804 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
805 contour_cube: Cube, optional
806 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
808 Raises
809 ------
810 ValueError
811 If the cube doesn't have the right dimensions.
812 """
813 # Use the smallest square grid that will fit the members.
814 nmember = len(cube.coord(stamp_coordinate).points)
815 grid_rows = int(math.sqrt(nmember))
816 grid_size = math.ceil(nmember / grid_rows)
818 fig = plt.figure(
819 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
820 )
822 # Specify the color bar
823 cmap, levels, norm = colorbar_map_levels(cube)
824 # If overplotting, set required colorbars
825 if overlay_cube: 825 ↛ 826line 825 didn't jump to line 826 because the condition on line 825 was never true
826 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
827 if contour_cube: 827 ↛ 828line 827 didn't jump to line 828 because the condition on line 827 was never true
828 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
830 # Make a subplot for each member.
831 for member, subplot in zip(
832 cube.slices_over(stamp_coordinate),
833 range(1, grid_size * grid_rows + 1),
834 strict=False,
835 ):
836 # Setup subplot map projection, extent and coastlines and borderlines.
837 axes = _setup_spatial_map(
838 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
839 )
840 if method == "contourf":
841 # Filled contour plot of the field.
842 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
843 elif method == "pcolormesh":
844 if levels is not None:
845 vmin = min(levels)
846 vmax = max(levels)
847 else:
848 raise TypeError("Unknown vmin and vmax range.")
849 vmin, vmax = None, None
850 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
851 # if levels are defined.
852 if norm is not None: 852 ↛ 853line 852 didn't jump to line 853 because the condition on line 852 was never true
853 vmin = None
854 vmax = None
855 # pcolormesh plot of the field.
856 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
857 else:
858 raise ValueError(f"Unknown plotting method: {method}")
860 # Overplot overlay field, if required
861 if overlay_cube: 861 ↛ 862line 861 didn't jump to line 862 because the condition on line 861 was never true
862 try:
863 over_vmin = min(over_levels)
864 over_vmax = max(over_levels)
865 except TypeError:
866 over_vmin, over_vmax = None, None
867 if over_norm is not None:
868 over_vmin = None
869 over_vmax = None
870 iplt.pcolormesh(
871 overlay_cube[member.coord(stamp_coordinate).points[0]],
872 cmap=over_cmap,
873 norm=over_norm,
874 alpha=0.6,
875 vmin=over_vmin,
876 vmax=over_vmax,
877 )
878 # Overplot contour field, if required
879 if contour_cube: 879 ↛ 880line 879 didn't jump to line 880 because the condition on line 879 was never true
880 iplt.contour(
881 contour_cube[member.coord(stamp_coordinate).points[0]],
882 colors="darkgray",
883 levels=cntr_levels,
884 norm=cntr_norm,
885 alpha=0.6,
886 linestyles="--",
887 linewidths=1,
888 )
889 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
890 axes.set_title(f"{mtitle}")
892 # Put the shared colorbar in its own axes.
893 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
894 colorbar = fig.colorbar(
895 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
896 )
897 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
899 # Overall figure title.
900 fig.suptitle(title, fontsize=16)
902 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
903 logging.info("Saved contour postage stamp plot to %s", filename)
904 plt.close(fig)
907def _plot_and_save_line_series(
908 cubes: iris.cube.CubeList,
909 coords: list[iris.coords.Coord],
910 ensemble_coord: str,
911 filename: str,
912 title: str,
913 **kwargs,
914):
915 """Plot and save a 1D line series.
917 Parameters
918 ----------
919 cubes: Cube or CubeList
920 Cube or CubeList containing the cubes to plot on the y-axis.
921 coords: list[Coord]
922 Coordinates to plot on the x-axis, one per cube.
923 ensemble_coord: str
924 Ensemble coordinate in the cube.
925 filename: str
926 Filename of the plot to write.
927 title: str
928 Plot title.
929 """
930 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
932 model_colors_map = get_model_colors_map(cubes)
934 # Store min/max ranges.
935 y_levels = []
937 # Check match-up across sequence coords gives consistent sizes
938 validate_cubes_coords(cubes, coords)
940 for cube, coord in zip(cubes, coords, strict=True):
941 label = None
942 color = "black"
943 if model_colors_map:
944 label = cube.attributes.get("model_name")
945 color = model_colors_map.get(label)
946 for cube_slice in cube.slices_over(ensemble_coord):
947 # Label with (control) if part of an ensemble or not otherwise.
948 if cube_slice.coord(ensemble_coord).points == [0]:
949 iplt.plot(
950 coord,
951 cube_slice,
952 color=color,
953 marker="o",
954 ls="-",
955 lw=3,
956 label=f"{label} (control)"
957 if len(cube.coord(ensemble_coord).points) > 1
958 else label,
959 )
960 # Label with (perturbed) if part of an ensemble and not the control.
961 else:
962 iplt.plot(
963 coord,
964 cube_slice,
965 color=color,
966 ls="-",
967 lw=1.5,
968 alpha=0.75,
969 label=f"{label} (member)",
970 )
972 # Calculate the global min/max if multiple cubes are given.
973 _, levels, _ = colorbar_map_levels(cube, axis="y")
974 if levels is not None: 974 ↛ 975line 974 didn't jump to line 975 because the condition on line 974 was never true
975 y_levels.append(min(levels))
976 y_levels.append(max(levels))
978 # Get the current axes.
979 ax = plt.gca()
981 # Add some labels and tweak the style.
982 # check if cubes[0] works for single cube if not CubeList
983 if coords[0].name() == "time":
984 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
985 else:
986 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
987 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
988 ax.set_title(title, fontsize=16)
990 ax.ticklabel_format(axis="y", useOffset=False)
991 ax.tick_params(axis="x", labelrotation=15)
992 ax.tick_params(axis="both", labelsize=12)
994 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
995 if y_levels: 995 ↛ 996line 995 didn't jump to line 996 because the condition on line 995 was never true
996 ax.set_ylim(min(y_levels), max(y_levels))
997 # Add zero line.
998 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
999 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1000 logging.debug(
1001 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1002 )
1003 else:
1004 ax.autoscale()
1006 # Add gridlines
1007 ax.grid(linestyle="--", color="grey", linewidth=1)
1008 # Ientify unique labels for legend
1009 handles = list(
1010 {
1011 label: handle
1012 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1013 }.values()
1014 )
1015 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1017 # Save plot.
1018 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1019 logging.info("Saved line plot to %s", filename)
1020 plt.close(fig)
1023def _plot_and_save_line_power_spectrum_series(
1024 cubes: iris.cube.Cube | iris.cube.CubeList,
1025 coords: list[iris.coords.Coord],
1026 ensemble_coord: str,
1027 filename: str,
1028 title: str,
1029 series_coordinate: str,
1030 **kwargs,
1031):
1032 """Plot and save a 1D line series.
1034 Parameters
1035 ----------
1036 cubes: Cube or CubeList
1037 Cube or CubeList containing the cubes to plot on the y-axis.
1038 coords: list[Coord]
1039 Coordinates to plot on the x-axis, one per cube.
1040 ensemble_coord: str
1041 Ensemble coordinate in the cube.
1042 filename: str
1043 Filename of the plot to write.
1044 title: str
1045 Plot title.
1046 series_coordinate: str
1047 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1048 """
1049 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1050 model_colors_map = get_model_colors_map(cubes)
1051 ax = plt.gca()
1053 # Store min/max ranges.
1054 y_levels = []
1056 line_marker = None
1057 line_width = 1
1059 for cube in iter_maybe(cubes):
1060 # next 2 lines replace chunk of code.
1061 xcoord = _select_series_coord(cube, series_coordinate)
1062 xname = xcoord.points
1064 yfield = cube.data # power spectrum
1065 label = None
1066 color = "black"
1067 if model_colors_map: 1067 ↛ 1070line 1067 didn't jump to line 1070 because the condition on line 1067 was always true
1068 label = cube.attributes.get("model_name")
1069 color = model_colors_map.get(label)
1070 for cube_slice in cube.slices_over(ensemble_coord):
1071 # Label with (control) if part of an ensemble or not otherwise.
1072 if cube_slice.coord(ensemble_coord).points == [0]: 1072 ↛ 1086line 1072 didn't jump to line 1086 because the condition on line 1072 was always true
1073 ax.plot(
1074 xname,
1075 yfield,
1076 color=color,
1077 marker=line_marker,
1078 ls="-",
1079 lw=line_width,
1080 label=f"{label} (control)"
1081 if len(cube.coord(ensemble_coord).points) > 1
1082 else label,
1083 )
1084 # Label with (perturbed) if part of an ensemble and not the control.
1085 else:
1086 ax.plot(
1087 xname,
1088 yfield,
1089 color=color,
1090 ls="-",
1091 lw=1.5,
1092 alpha=0.75,
1093 label=f"{label} (member)",
1094 )
1096 # Calculate the global min/max if multiple cubes are given.
1097 _, levels, _ = colorbar_map_levels(cube, axis="y")
1098 if levels is not None: 1098 ↛ 1099line 1098 didn't jump to line 1099 because the condition on line 1098 was never true
1099 y_levels.append(min(levels))
1100 y_levels.append(max(levels))
1102 # Add some labels and tweak the style.
1104 title = f"{title}"
1105 ax.set_title(title, fontsize=16)
1107 # Set appropriate x-axis label based on coordinate
1108 if series_coordinate == "wavelength" or ( 1108 ↛ 1111line 1108 didn't jump to line 1111 because the condition on line 1108 was never true
1109 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1110 ):
1111 ax.set_xlabel("Wavelength (km)", fontsize=14)
1112 elif series_coordinate == "physical_wavenumber" or ( 1112 ↛ 1115line 1112 didn't jump to line 1115 because the condition on line 1112 was never true
1113 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1114 ):
1115 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1116 else: # frequency or check units
1117 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1117 ↛ 1118line 1117 didn't jump to line 1118 because the condition on line 1117 was never true
1118 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1119 else:
1120 ax.set_xlabel("Wavenumber", fontsize=14)
1122 ax.set_ylabel("Power Spectral Density", fontsize=14)
1123 ax.tick_params(axis="both", labelsize=12)
1125 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1127 # Set log-log scale
1128 ax.set_xscale("log")
1129 ax.set_yscale("log")
1131 # Add gridlines
1132 ax.grid(linestyle="--", color="grey", linewidth=1)
1133 # Ientify unique labels for legend
1134 handles = list(
1135 {
1136 label: handle
1137 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1138 }.values()
1139 )
1140 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1142 # Save plot.
1143 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1144 logging.info("Saved line plot to %s", filename)
1145 plt.close(fig)
1148def _plot_and_save_vertical_line_series(
1149 cubes: iris.cube.CubeList,
1150 coords: list[iris.coords.Coord],
1151 ensemble_coord: str,
1152 filename: str,
1153 series_coordinate: str,
1154 title: str,
1155 vmin: float,
1156 vmax: float,
1157 **kwargs,
1158):
1159 """Plot and save a 1D line series in vertical.
1161 Parameters
1162 ----------
1163 cubes: CubeList
1164 1 dimensional Cube or CubeList of the data to plot on x-axis.
1165 coord: list[Coord]
1166 Coordinates to plot on the y-axis, one per cube.
1167 ensemble_coord: str
1168 Ensemble coordinate in the cube.
1169 filename: str
1170 Filename of the plot to write.
1171 series_coordinate: str
1172 Coordinate to use as vertical axis.
1173 title: str
1174 Plot title.
1175 vmin: float
1176 Minimum value for the x-axis.
1177 vmax: float
1178 Maximum value for the x-axis.
1179 """
1180 # plot the vertical pressure axis using log scale
1181 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1183 model_colors_map = get_model_colors_map(cubes)
1185 # Check match-up across sequence coords gives consistent sizes
1186 validate_cubes_coords(cubes, coords)
1188 for cube, coord in zip(cubes, coords, strict=True):
1189 label = None
1190 color = "black"
1191 if model_colors_map: 1191 ↛ 1192line 1191 didn't jump to line 1192 because the condition on line 1191 was never true
1192 label = cube.attributes.get("model_name")
1193 color = model_colors_map.get(label)
1195 for cube_slice in cube.slices_over(ensemble_coord):
1196 # If ensemble data given plot control member with (control)
1197 # unless single forecast.
1198 if cube_slice.coord(ensemble_coord).points == [0]:
1199 iplt.plot(
1200 cube_slice,
1201 coord,
1202 color=color,
1203 marker="o",
1204 ls="-",
1205 lw=3,
1206 label=f"{label} (control)"
1207 if len(cube.coord(ensemble_coord).points) > 1
1208 else label,
1209 )
1210 # If ensemble data given plot perturbed members with (perturbed).
1211 else:
1212 iplt.plot(
1213 cube_slice,
1214 coord,
1215 color=color,
1216 ls="-",
1217 lw=1.5,
1218 alpha=0.75,
1219 label=f"{label} (member)",
1220 )
1222 # Get the current axis
1223 ax = plt.gca()
1225 # Special handling for pressure level data.
1226 if series_coordinate == "pressure": 1226 ↛ 1248line 1226 didn't jump to line 1248 because the condition on line 1226 was always true
1227 # Invert y-axis and set to log scale.
1228 ax.invert_yaxis()
1229 ax.set_yscale("log")
1231 # Define y-ticks and labels for pressure log axis.
1232 y_tick_labels = [
1233 "1000",
1234 "850",
1235 "700",
1236 "500",
1237 "300",
1238 "200",
1239 "100",
1240 ]
1241 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1243 # Set y-axis limits and ticks.
1244 ax.set_ylim(1100, 100)
1246 # Test if series_coordinate is model level data. The UM data uses
1247 # model_level_number and lfric uses full_levels as coordinate.
1248 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1249 # Define y-ticks and labels for vertical axis.
1250 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1251 y_tick_labels = [str(int(i)) for i in y_ticks]
1252 ax.set_ylim(min(y_ticks), max(y_ticks))
1254 ax.set_yticks(y_ticks)
1255 ax.set_yticklabels(y_tick_labels)
1257 # Set x-axis limits.
1258 ax.set_xlim(vmin, vmax)
1259 # Mark y=0 if present in plot.
1260 if vmin < 0.0 and vmax > 0.0: 1260 ↛ 1261line 1260 didn't jump to line 1261 because the condition on line 1260 was never true
1261 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1263 # Add some labels and tweak the style.
1264 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1265 ax.set_xlabel(
1266 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1267 )
1268 ax.set_title(title, fontsize=16)
1269 ax.ticklabel_format(axis="x")
1270 ax.tick_params(axis="y")
1271 ax.tick_params(axis="both", labelsize=12)
1273 # Add gridlines
1274 ax.grid(linestyle="--", color="grey", linewidth=1)
1275 # Ientify unique labels for legend
1276 handles = list(
1277 {
1278 label: handle
1279 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1280 }.values()
1281 )
1282 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1284 # Save plot.
1285 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1286 logging.info("Saved line plot to %s", filename)
1287 plt.close(fig)
1290def _plot_and_save_scatter_plot(
1291 cube_x: iris.cube.Cube | iris.cube.CubeList,
1292 cube_y: iris.cube.Cube | iris.cube.CubeList,
1293 filename: str,
1294 title: str,
1295 one_to_one: bool,
1296 model_names: list[str] = None,
1297 **kwargs,
1298):
1299 """Plot and save a 2D scatter plot.
1301 Parameters
1302 ----------
1303 cube_x: Cube | CubeList
1304 1 dimensional Cube or CubeList of the data to plot on x-axis.
1305 cube_y: Cube | CubeList
1306 1 dimensional Cube or CubeList of the data to plot on y-axis.
1307 filename: str
1308 Filename of the plot to write.
1309 title: str
1310 Plot title.
1311 one_to_one: bool
1312 Whether a 1:1 line is plotted.
1313 """
1314 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1315 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1316 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1317 # over the pairs simultaneously.
1319 # Ensure cube_x and cube_y are iterable
1320 cube_x_iterable = iter_maybe(cube_x)
1321 cube_y_iterable = iter_maybe(cube_y)
1323 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1324 iplt.scatter(cube_x_iter, cube_y_iter)
1325 if one_to_one is True:
1326 plt.plot(
1327 [
1328 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1329 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1330 ],
1331 [
1332 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1333 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1334 ],
1335 "k",
1336 linestyle="--",
1337 )
1338 ax = plt.gca()
1340 # Add some labels and tweak the style.
1341 if model_names is None:
1342 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1343 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1344 else:
1345 # Add the model names, these should be order of base (x) and other (y).
1346 ax.set_xlabel(
1347 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1348 )
1349 ax.set_ylabel(
1350 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1351 )
1352 ax.set_title(title, fontsize=16)
1353 ax.ticklabel_format(axis="y", useOffset=False)
1354 ax.tick_params(axis="x", labelrotation=15)
1355 ax.tick_params(axis="both", labelsize=12)
1356 ax.autoscale()
1358 # Save plot.
1359 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1360 logging.info("Saved scatter plot to %s", filename)
1361 plt.close(fig)
1364def _plot_and_save_vector_plot(
1365 cube_u: iris.cube.Cube,
1366 cube_v: iris.cube.Cube,
1367 filename: str,
1368 title: str,
1369 method: Literal["contourf", "pcolormesh"],
1370 **kwargs,
1371):
1372 """Plot and save a 2D vector plot.
1374 Parameters
1375 ----------
1376 cube_u: Cube
1377 2 dimensional Cube of u component of the data.
1378 cube_v: Cube
1379 2 dimensional Cube of v component of the data.
1380 filename: str
1381 Filename of the plot to write.
1382 title: str
1383 Plot title.
1384 """
1385 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1386 # Create a cube containing the magnitude of the vector field.
1387 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1388 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1389 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1390 cube_vec_mag.rename(
1391 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1392 )
1394 # Specify the color bar
1395 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1397 # Setup plot map projection, extent and coastlines and borderlines.
1398 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1400 if method == "contourf":
1401 # Filled contour plot of the field.
1402 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1403 elif method == "pcolormesh":
1404 try:
1405 vmin = min(levels)
1406 vmax = max(levels)
1407 except TypeError:
1408 vmin, vmax = None, None
1409 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1410 # if levels are defined.
1411 if norm is not None:
1412 vmin = None
1413 vmax = None
1414 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1415 else:
1416 raise ValueError(f"Unknown plotting method: {method}")
1418 # Check to see if transect, and if so, adjust y axis.
1419 if is_transect(cube_vec_mag):
1420 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1421 axes.invert_yaxis()
1422 axes.set_yscale("log")
1423 axes.set_ylim(1100, 100)
1424 # If both model_level_number and level_height exists, iplt can construct
1425 # plot as a function of height above orography (NOT sea level).
1426 elif {"model_level_number", "level_height"}.issubset(
1427 {coord.name() for coord in cube_vec_mag.coords()}
1428 ):
1429 axes.set_yscale("log")
1431 axes.set_title(
1432 f"{title}\n"
1433 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1434 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1435 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1436 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1437 fontsize=16,
1438 )
1440 else:
1441 # Add title.
1442 axes.set_title(title, fontsize=16)
1444 # Add watermark with min/max/mean. Currently not user togglable.
1445 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1446 axes.annotate(
1447 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}",
1448 xy=(0.05, -0.05),
1449 xycoords="axes fraction",
1450 xytext=(-5, 5),
1451 textcoords="offset points",
1452 ha="right",
1453 va="bottom",
1454 size=11,
1455 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1456 )
1458 # Add colour bar.
1459 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1460 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1461 # add ticks and tick_labels for every levels if less than 20 levels exist
1462 if levels is not None and len(levels) < 20:
1463 cbar.set_ticks(levels)
1464 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1466 # 30 barbs along the longest axis of the plot, or a barb per point for data
1467 # with less than 30 points.
1468 step = max(max(cube_u.shape) // 30, 1)
1469 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1471 # Save plot.
1472 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1473 logging.info("Saved vector plot to %s", filename)
1474 plt.close(fig)
1477def _plot_and_save_histogram_series(
1478 cubes: iris.cube.Cube | iris.cube.CubeList,
1479 filename: str,
1480 title: str,
1481 vmin: float,
1482 vmax: float,
1483 **kwargs,
1484):
1485 """Plot and save a histogram series.
1487 Parameters
1488 ----------
1489 cubes: Cube or CubeList
1490 2 dimensional Cube or CubeList of the data to plot as histogram.
1491 filename: str
1492 Filename of the plot to write.
1493 title: str
1494 Plot title.
1495 vmin: float
1496 minimum for colorbar
1497 vmax: float
1498 maximum for colorbar
1499 """
1500 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1501 ax = plt.gca()
1503 model_colors_map = get_model_colors_map(cubes)
1505 # Set default that histograms will produce probability density function
1506 # at each bin (integral over range sums to 1).
1507 density = True
1509 for cube in iter_maybe(cubes):
1510 # Easier to check title (where var name originates)
1511 # than seeing if long names exist etc.
1512 # Exception case, where distribution better fits log scales/bins.
1513 if (
1514 ("surface_microphysical" in title)
1515 or ("rain accumulation" in title)
1516 or ("Rainfall rate Composite" in title)
1517 or ("Nimrod_5min" in title)
1518 ):
1519 if "amount" in title: 1519 ↛ 1521line 1519 didn't jump to line 1521 because the condition on line 1519 was never true
1520 # Compute histogram following Klingaman et al. (2017): ASoP
1521 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1522 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1523 density = False
1524 else:
1525 bins = 10.0 ** (
1526 np.arange(-10, 27, 1) / 10.0
1527 ) # Suggestion from RMED toolbox.
1528 bins = np.insert(bins, 0, 0)
1529 ax.set_yscale("log")
1530 vmin = bins[1]
1531 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1532 ax.set_xscale("log")
1533 elif "lightning" in title:
1534 bins = [0, 1, 2, 3, 4, 5]
1535 else:
1536 bins = np.linspace(vmin, vmax, 51)
1537 logging.debug(
1538 "Plotting histogram with %s bins %s - %s.",
1539 np.size(bins),
1540 np.min(bins),
1541 np.max(bins),
1542 )
1544 # Reshape cube data into a single array to allow for a single histogram.
1545 # Otherwise we plot xdim histograms stacked.
1546 cube_data_1d = (cube.data).flatten()
1548 label = None
1549 color = "black"
1550 if model_colors_map:
1551 label = cube.attributes.get("model_name")
1552 color = model_colors_map[label]
1553 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1555 # Compute area under curve.
1556 if ( 1556 ↛ 1562line 1556 didn't jump to line 1562 because the condition on line 1556 was never true
1557 ("surface_microphysical" in title and "amount" in title)
1558 or ("rain_accumulation" in title)
1559 or ("Rainfall rate Composite" in title)
1560 or ("Nimrod_5min" in title)
1561 ):
1562 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1563 x = x * bin_mean / x.sum()
1564 x = x[1:]
1565 y = y[1:]
1567 ax.plot(
1568 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1569 )
1571 # Add some labels and tweak the style.
1572 ax.set_title(title, fontsize=16)
1573 ax.set_xlabel(
1574 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1575 )
1576 ax.set_ylabel("Normalised probability density", fontsize=14)
1577 if ( 1577 ↛ 1582line 1577 didn't jump to line 1582 because the condition on line 1577 was never true
1578 ("surface_microphysical" in title and "amount" in title)
1579 or ("rain accumulation" in title)
1580 or ("Nimrod_5min" in title)
1581 ):
1582 ax.set_ylabel(
1583 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1584 )
1585 ax.set_xlim(vmin, vmax)
1586 ax.tick_params(axis="both", labelsize=12)
1588 # Overlay grid-lines onto histogram plot.
1589 ax.grid(linestyle="--", color="grey", linewidth=1)
1590 if model_colors_map:
1591 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1593 # Save plot.
1594 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1595 logging.info("Saved histogram plot to %s", filename)
1596 plt.close(fig)
1599def _plot_and_save_postage_stamp_histogram_series(
1600 cube: iris.cube.Cube,
1601 filename: str,
1602 title: str,
1603 stamp_coordinate: str,
1604 vmin: float,
1605 vmax: float,
1606 **kwargs,
1607):
1608 """Plot and save postage (ensemble members) stamps for a histogram series.
1610 Parameters
1611 ----------
1612 cube: Cube
1613 2 dimensional Cube of the data to plot as histogram.
1614 filename: str
1615 Filename of the plot to write.
1616 title: str
1617 Plot title.
1618 stamp_coordinate: str
1619 Coordinate that becomes different plots.
1620 vmin: float
1621 minimum for pdf x-axis
1622 vmax: float
1623 maximum for pdf x-axis
1624 """
1625 # Use the smallest square grid that will fit the members.
1626 nmember = len(cube.coord(stamp_coordinate).points)
1627 grid_rows = int(math.sqrt(nmember))
1628 grid_size = math.ceil(nmember / grid_rows)
1630 fig = plt.figure(
1631 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1632 )
1633 # Make a subplot for each member.
1634 for member, subplot in zip(
1635 cube.slices_over(stamp_coordinate),
1636 range(1, grid_size * grid_rows + 1),
1637 strict=False,
1638 ):
1639 # Implicit interface is much easier here, due to needing to have the
1640 # cartopy GeoAxes generated.
1641 plt.subplot(grid_rows, grid_size, subplot)
1642 # Reshape cube data into a single array to allow for a single histogram.
1643 # Otherwise we plot xdim histograms stacked.
1644 member_data_1d = (member.data).flatten()
1645 plt.hist(member_data_1d, density=True, stacked=True)
1646 axes = plt.gca()
1647 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1648 axes.set_title(f"{mtitle}")
1649 axes.set_xlim(vmin, vmax)
1651 # Overall figure title.
1652 fig.suptitle(title, fontsize=16)
1654 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1655 logging.info("Saved histogram postage stamp plot to %s", filename)
1656 plt.close(fig)
1659def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1660 cube: iris.cube.Cube,
1661 filename: str,
1662 title: str,
1663 stamp_coordinate: str,
1664 vmin: float,
1665 vmax: float,
1666 **kwargs,
1667):
1668 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1669 ax.set_title(title, fontsize=16)
1670 ax.set_xlim(vmin, vmax)
1671 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1672 ax.set_ylabel("normalised probability density", fontsize=14)
1673 # Loop over all slices along the stamp_coordinate
1674 for member in cube.slices_over(stamp_coordinate):
1675 # Flatten the member data to 1D
1676 member_data_1d = member.data.flatten()
1677 # Plot the histogram using plt.hist
1678 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1679 plt.hist(
1680 member_data_1d,
1681 density=True,
1682 stacked=True,
1683 label=f"{mtitle}",
1684 )
1686 # Add a legend
1687 ax.legend(fontsize=16)
1689 # Save the figure to a file
1690 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1691 logging.info("Saved histogram postage stamp plot to %s", filename)
1693 # Close the figure
1694 plt.close(fig)
1697def _spatial_plot(
1698 method: Literal["contourf", "pcolormesh", "scatter"],
1699 cube: iris.cube.Cube,
1700 filename: str | None,
1701 sequence_coordinate: str,
1702 stamp_coordinate: str,
1703 overlay_cube: iris.cube.Cube | None = None,
1704 contour_cube: iris.cube.Cube | None = None,
1705 point_cube: iris.cube.Cube | None = None,
1706 **kwargs,
1707):
1708 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1710 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1711 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1712 is present then postage stamp plots will be produced.
1714 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1715 be overplotted on the same figure.
1717 Parameters
1718 ----------
1719 method: "contourf" | "pcolormesh" | "scatter"
1720 The plotting method to use.
1721 Select choice of "contourf" or "pcolormesh" for gridded data.
1722 Use "scatter" for point-based data.
1723 cube: Cube
1724 Iris cube of the data to plot. It should have two spatial dimensions,
1725 such as lat and lon, and may also have a another two dimension to be
1726 plotted sequentially and/or as postage stamp plots.
1727 filename: str | None
1728 Name of the plot to write, used as a prefix for plot sequences. If None
1729 uses the recipe name.
1730 sequence_coordinate: str
1731 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1732 This coordinate must exist in the cube.
1733 stamp_coordinate: str
1734 Coordinate about which to plot postage stamp plots. Defaults to
1735 ``"realization"``.
1736 overlay_cube: Cube | None, optional
1737 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1738 contour_cube: Cube | None, optional
1739 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1740 point_cube: Cube | None, optional
1741 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
1743 Raises
1744 ------
1745 ValueError
1746 If the cube doesn't have the right dimensions.
1747 TypeError
1748 If the cube isn't a single cube.
1749 """
1750 # Ensure we've got a single cube.
1751 cube = check_single_cube(cube)
1753 # Set title based on recipe metadata or use cube name
1754 recipe_title = get_recipe_metadata().get("title", cube.name())
1756 # Check if there is a valid stamp coordinate in cube dimensions.
1757 if stamp_coordinate == "realization": 1757 ↛ 1762line 1757 didn't jump to line 1762 because the condition on line 1757 was always true
1758 stamp_coordinate = check_stamp_coordinate(cube)
1760 # Make postage stamp plots if stamp_coordinate exists and has more than a
1761 # single point.
1762 plotting_func = _plot_and_save_spatial_plot
1763 try:
1764 if cube.coord(stamp_coordinate).shape[0] > 1:
1765 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1766 except iris.exceptions.CoordinateNotFoundError:
1767 pass
1769 # Produce a geographical scatter plot if the data have a
1770 # dimension called observation or model_obs_error
1771 if any( 1771 ↛ 1775line 1771 didn't jump to line 1775 because the condition on line 1771 was never true
1772 crd.var_name == "station" or crd.var_name == "model_obs_error"
1773 for crd in cube.coords()
1774 ):
1775 plotting_func = _plot_and_save_spatial_plot
1776 method = "scatter"
1778 # Must have a sequence coordinate.
1779 try:
1780 cube.coord(sequence_coordinate)
1781 except iris.exceptions.CoordinateNotFoundError as err:
1782 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1784 # Create a plot for each value of the sequence coordinate.
1785 plot_index = []
1786 nplot = np.size(cube.coord(sequence_coordinate).points)
1788 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1789 # Set plot titles and filename
1790 seq_coord = cube_slice.coord(sequence_coordinate)
1791 plot_title, plot_filename = _set_title_and_filename(
1792 seq_coord, nplot, recipe_title, filename
1793 )
1795 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1796 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1797 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1798 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1800 # Do the actual plotting.
1801 plotting_func(
1802 cube_slice,
1803 filename=plot_filename,
1804 stamp_coordinate=stamp_coordinate,
1805 title=plot_title,
1806 method=method,
1807 overlay_cube=overlay_slice,
1808 contour_cube=contour_slice,
1809 point_cube=point_slice,
1810 **kwargs,
1811 )
1812 plot_index.append(plot_filename)
1814 # Add list of plots to plot metadata.
1815 complete_plot_index = _append_to_plot_index(plot_index)
1817 # Make a page to display the plots.
1818 _make_plot_html_page(complete_plot_index)
1821####################
1822# Public functions #
1823####################
1826def spatial_contour_plot(
1827 cube: iris.cube.Cube,
1828 filename: str = None,
1829 sequence_coordinate: str = "time",
1830 stamp_coordinate: str = "realization",
1831 **kwargs,
1832) -> iris.cube.Cube:
1833 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1835 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1836 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1837 is present then postage stamp plots will be produced.
1839 Parameters
1840 ----------
1841 cube: Cube
1842 Iris cube of the data to plot. It should have two spatial dimensions,
1843 such as lat and lon, and may also have a another two dimension to be
1844 plotted sequentially and/or as postage stamp plots.
1845 filename: str, optional
1846 Name of the plot to write, used as a prefix for plot sequences. Defaults
1847 to the recipe name.
1848 sequence_coordinate: str, optional
1849 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1850 This coordinate must exist in the cube.
1851 stamp_coordinate: str, optional
1852 Coordinate about which to plot postage stamp plots. Defaults to
1853 ``"realization"``.
1855 Returns
1856 -------
1857 Cube
1858 The original cube (so further operations can be applied).
1860 Raises
1861 ------
1862 ValueError
1863 If the cube doesn't have the right dimensions.
1864 TypeError
1865 If the cube isn't a single cube.
1866 """
1867 _spatial_plot(
1868 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1869 )
1870 return cube
1873def spatial_pcolormesh_plot(
1874 cube: iris.cube.Cube,
1875 filename: str = None,
1876 sequence_coordinate: str = "time",
1877 stamp_coordinate: str = "realization",
1878 **kwargs,
1879) -> iris.cube.Cube:
1880 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1882 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1883 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1884 is present then postage stamp plots will be produced.
1886 This function is significantly faster than ``spatial_contour_plot``,
1887 especially at high resolutions, and should be preferred unless contiguous
1888 contour areas are important.
1890 Parameters
1891 ----------
1892 cube: Cube
1893 Iris cube of the data to plot. It should have two spatial dimensions,
1894 such as lat and lon, and may also have a another two dimension to be
1895 plotted sequentially and/or as postage stamp plots.
1896 filename: str, optional
1897 Name of the plot to write, used as a prefix for plot sequences. Defaults
1898 to the recipe name.
1899 sequence_coordinate: str, optional
1900 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1901 This coordinate must exist in the cube.
1902 stamp_coordinate: str, optional
1903 Coordinate about which to plot postage stamp plots. Defaults to
1904 ``"realization"``.
1906 Returns
1907 -------
1908 Cube
1909 The original cube (so further operations can be applied).
1911 Raises
1912 ------
1913 ValueError
1914 If the cube doesn't have the right dimensions.
1915 TypeError
1916 If the cube isn't a single cube.
1917 """
1918 _spatial_plot(
1919 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1920 )
1921 return cube
1924def spatial_multi_pcolormesh_plot(
1925 cube: iris.cube.Cube,
1926 overlay_cube: iris.cube.Cube | None = None,
1927 contour_cube: iris.cube.Cube | None = None,
1928 point_cube: iris.cube.Cube | None = None,
1929 filename: str = None,
1930 sequence_coordinate: str = "time",
1931 stamp_coordinate: str = "realization",
1932 **kwargs,
1933) -> iris.cube.Cube:
1934 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1936 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1937 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1938 is present then postage stamp plots will be produced.
1940 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1942 If specified, contours of a contour_cube can be overplotted on top of those.
1944 If specified, a spatial scatter map of point_cube can be overplotted.
1946 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1948 This function is significantly faster than ``spatial_contour_plot``,
1949 especially at high resolutions, and should be preferred unless contiguous
1950 contour areas are important.
1952 Parameters
1953 ----------
1954 cube: Cube
1955 Iris cube of the data to plot. It should have two spatial dimensions,
1956 such as lat and lon, and may also have two additional dimensions to be
1957 plotted sequentially and/or as postage stamp plots.
1958 overlay_cube: Cube, optional
1959 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1960 such as lat and lon, and may also have two additional dimensions to be
1961 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.
1962 If not provided, output plot generated without overlay cube.
1963 contour_cube: Cube, optional
1964 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,
1965 such as lat and lon, and may also have two additional dimensions to be
1966 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
1967 point_cube: Cube, optional
1968 Iris cube of the data to plot as a scatter map overlay on top of basis cube (overlay_cube and/or contour_cube). It should have two
1969 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
1970 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
1971 postage stamp plots. If not provided, output plot generated without point-based layer.
1972 filename: str, optional
1973 Name of the plot to write, used as a prefix for plot sequences. Defaults
1974 to the recipe name.
1975 sequence_coordinate: str, optional
1976 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1977 This coordinate must exist in the cube.
1978 stamp_coordinate: str, optional
1979 Coordinate about which to plot postage stamp plots. Defaults to
1980 ``"realization"``.
1982 Returns
1983 -------
1984 Cube
1985 The original cube (so further operations can be applied).
1987 Raises
1988 ------
1989 ValueError
1990 If the cube doesn't have the right dimensions.
1991 TypeError
1992 If the cube isn't a single cube.
1993 """
1994 _spatial_plot(
1995 "pcolormesh",
1996 cube,
1997 filename,
1998 sequence_coordinate,
1999 stamp_coordinate,
2000 overlay_cube=overlay_cube,
2001 contour_cube=contour_cube,
2002 point_cube=point_cube,
2003 )
2004 return cube, overlay_cube, contour_cube, point_cube
2007# TODO: Expand function to handle ensemble data.
2008# line_coordinate: str, optional
2009# Coordinate about which to plot multiple lines. Defaults to
2010# ``"realization"``.
2011def plot_line_series(
2012 cube: iris.cube.Cube | iris.cube.CubeList,
2013 filename: str = None,
2014 series_coordinate: str = "time",
2015 sequence_coordinate: str = "time",
2016 # add the following for ensembles
2017 stamp_coordinate: str = "realization",
2018 single_plot: bool = False,
2019 **kwargs,
2020) -> iris.cube.Cube | iris.cube.CubeList:
2021 """Plot a line plot for the specified coordinate.
2023 The Cube or CubeList must be 1D.
2025 Parameters
2026 ----------
2027 iris.cube | iris.cube.CubeList
2028 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2029 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2030 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2031 filename: str, optional
2032 Name of the plot to write, used as a prefix for plot sequences. Defaults
2033 to the recipe name.
2034 series_coordinate: str, optional
2035 Coordinate about which to make a series. Defaults to ``"time"``. This
2036 coordinate must exist in the cube.
2038 Returns
2039 -------
2040 iris.cube.Cube | iris.cube.CubeList
2041 The original Cube or CubeList (so further operations can be applied).
2043 Raises
2044 ------
2045 ValueError
2046 If the cubes don't have the right dimensions.
2047 TypeError
2048 If the cube isn't a Cube or CubeList.
2049 """
2050 # Ensure we have a name for the plot file.
2051 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2053 num_models = get_num_models(cube)
2055 validate_cube_shape(cube, num_models)
2057 # Iterate over all cubes and extract coordinate to plot.
2058 cubes = iter_maybe(cube)
2059 coords = []
2060 for cube in cubes:
2061 try:
2062 coords.append(cube.coord(series_coordinate))
2063 except iris.exceptions.CoordinateNotFoundError as err:
2064 raise ValueError(
2065 f"Cube must have a {series_coordinate} coordinate."
2066 ) from err
2067 if cube.coords("realization"): 2067 ↛ 2071line 2067 didn't jump to line 2071 because the condition on line 2067 was always true
2068 if cube.ndim > 3: 2068 ↛ 2069line 2068 didn't jump to line 2069 because the condition on line 2068 was never true
2069 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2070 else:
2071 raise ValueError("Cube must have a realization coordinate.")
2073 plot_index = []
2075 # Check if this is a spectral plot by looking for spectral coordinates
2076 is_spectral_plot = series_coordinate in [
2077 "frequency",
2078 "physical_wavenumber",
2079 "wavelength",
2080 ]
2082 if is_spectral_plot:
2083 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2084 # coordinate frequency/wavenumber.
2085 # If several power spectra are plotted with time as sequence_coordinate for the
2086 # time slider option.
2088 # Internal plotting function.
2089 plotting_func = _plot_and_save_line_power_spectrum_series
2091 for cube in cubes:
2092 try:
2093 cube.coord(sequence_coordinate)
2094 except iris.exceptions.CoordinateNotFoundError as err:
2095 raise ValueError(
2096 f"Cube must have a {sequence_coordinate} coordinate."
2097 ) from err
2099 if num_models == 1: 2099 ↛ 2113line 2099 didn't jump to line 2113 because the condition on line 2099 was always true
2100 # check for ensembles
2101 if ( 2101 ↛ 2105line 2101 didn't jump to line 2105 because the condition on line 2101 was never true
2102 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2103 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2104 ):
2105 if single_plot:
2106 # Plot spectra, mean and ensemble spread on 1 plot
2107 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2108 else:
2109 # Plot postage stamps
2110 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2111 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2112 else:
2113 all_points = sorted(
2114 set(
2115 itertools.chain.from_iterable(
2116 cb.coord(sequence_coordinate).points for cb in cubes
2117 )
2118 )
2119 )
2120 all_slices = list(
2121 itertools.chain.from_iterable(
2122 cb.slices_over(sequence_coordinate) for cb in cubes
2123 )
2124 )
2125 # Matched slices (matched by seq coord point; it may happen that
2126 # evaluated models do not cover the same seq coord range, hence matching
2127 # necessary)
2128 cube_iterables = [
2129 iris.cube.CubeList(
2130 s
2131 for s in all_slices
2132 if s.coord(sequence_coordinate).points[0] == point
2133 )
2134 for point in all_points
2135 ]
2137 nplot = np.size(cube.coord(sequence_coordinate).points)
2139 # Create a plot for each value of the sequence coordinate. Allowing for
2140 # multiple cubes in a CubeList to be plotted in the same plot for similar
2141 # sequence values. Passing a CubeList into the internal plotting function
2142 # for similar values of the sequence coordinate. cube_slice can be an
2143 # iris.cube.Cube or an iris.cube.CubeList.
2145 for cube_slice in cube_iterables:
2146 # Normalize cube_slice to a list of cubes
2147 if isinstance(cube_slice, iris.cube.CubeList): 2147 ↛ 2148line 2147 didn't jump to line 2148 because the condition on line 2147 was never true
2148 cubes = list(cube_slice)
2149 elif isinstance(cube_slice, iris.cube.Cube): 2149 ↛ 2152line 2149 didn't jump to line 2152 because the condition on line 2149 was always true
2150 cubes = [cube_slice]
2151 else:
2152 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2154 # Use sequence value so multiple sequences can merge.
2155 seq_coord = cube_slice[0].coord(sequence_coordinate)
2156 plot_title, plot_filename = _set_title_and_filename(
2157 seq_coord, nplot, recipe_title, filename
2158 )
2160 # Format the coordinate value in a unit appropriate way.
2161 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2163 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2164 if nplot == 1 and seq_coord.has_bounds: 2164 ↛ 2169line 2164 didn't jump to line 2169 because the condition on line 2164 was always true
2165 if np.size(seq_coord.bounds) > 1: 2165 ↛ 2166line 2165 didn't jump to line 2166 because the condition on line 2165 was never true
2166 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2168 # Do the actual plotting.
2169 plotting_func(
2170 cube_slice,
2171 coords,
2172 stamp_coordinate,
2173 plot_filename,
2174 title,
2175 series_coordinate,
2176 )
2178 plot_index.append(plot_filename)
2179 else:
2180 # Format the title and filename using plotted series coordinate
2181 nplot = 1
2182 seq_coord = coords[0]
2183 plot_title, plot_filename = _set_title_and_filename(
2184 seq_coord, nplot, recipe_title, filename
2185 )
2186 # Do the actual plotting for all other series coordinate options.
2187 _plot_and_save_line_series(
2188 cubes, coords, stamp_coordinate, plot_filename, plot_title
2189 )
2191 plot_index.append(plot_filename)
2193 # append plot to list of plots
2194 complete_plot_index = _append_to_plot_index(plot_index)
2196 # Make a page to display the plots.
2197 _make_plot_html_page(complete_plot_index)
2199 return cube
2202def plot_vertical_line_series(
2203 cubes: iris.cube.Cube | iris.cube.CubeList,
2204 filename: str = None,
2205 series_coordinate: str = "model_level_number",
2206 sequence_coordinate: str = "time",
2207 # line_coordinate: str = "realization",
2208 **kwargs,
2209) -> iris.cube.Cube | iris.cube.CubeList:
2210 """Plot a line plot against a type of vertical coordinate.
2212 The Cube or CubeList must be 1D.
2214 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2215 then a sequence of plots will be produced.
2217 Parameters
2218 ----------
2219 iris.cube | iris.cube.CubeList
2220 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2221 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2222 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2223 filename: str, optional
2224 Name of the plot to write, used as a prefix for plot sequences. Defaults
2225 to the recipe name.
2226 series_coordinate: str, optional
2227 Coordinate to plot on the y-axis. Can be ``pressure`` or
2228 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2229 for LFRic. Defaults to ``model_level_number``.
2230 This coordinate must exist in the cube.
2231 sequence_coordinate: str, optional
2232 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2233 This coordinate must exist in the cube.
2235 Returns
2236 -------
2237 iris.cube.Cube | iris.cube.CubeList
2238 The original Cube or CubeList (so further operations can be applied).
2239 Plotted data.
2241 Raises
2242 ------
2243 ValueError
2244 If the cubes doesn't have the right dimensions.
2245 TypeError
2246 If the cube isn't a Cube or CubeList.
2247 """
2248 # Ensure we have a name for the plot file.
2249 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2251 cubes = iter_maybe(cubes)
2252 # Initialise empty list to hold all data from all cubes in a CubeList
2253 all_data = []
2255 # Store min/max ranges for x range.
2256 x_levels = []
2258 num_models = get_num_models(cubes)
2260 validate_cube_shape(cubes, num_models)
2262 # Iterate over all cubes in cube or CubeList and plot.
2263 coords = []
2264 for cube in cubes:
2265 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2266 try:
2267 coords.append(cube.coord(series_coordinate))
2268 except iris.exceptions.CoordinateNotFoundError as err:
2269 raise ValueError(
2270 f"Cube must have a {series_coordinate} coordinate."
2271 ) from err
2273 try:
2274 if cube.ndim > 1 or not cube.coords("realization"): 2274 ↛ 2282line 2274 didn't jump to line 2282 because the condition on line 2274 was always true
2275 cube.coord(sequence_coordinate)
2276 except iris.exceptions.CoordinateNotFoundError as err:
2277 raise ValueError(
2278 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2279 ) from err
2281 # Get minimum and maximum from levels information.
2282 _, levels, _ = colorbar_map_levels(cube, axis="x")
2283 if levels is not None: 2283 ↛ 2287line 2283 didn't jump to line 2287 because the condition on line 2283 was always true
2284 x_levels.append(min(levels))
2285 x_levels.append(max(levels))
2286 else:
2287 all_data.append(cube.data)
2289 if len(x_levels) == 0: 2289 ↛ 2291line 2289 didn't jump to line 2291 because the condition on line 2289 was never true
2290 # Combine all data into a single NumPy array
2291 combined_data = np.concatenate(all_data)
2293 # Set the lower and upper limit for the x-axis to ensure all plots have
2294 # same range. This needs to read the whole cube over the range of the
2295 # sequence and if applicable postage stamp coordinate.
2296 vmin = np.floor(combined_data.min())
2297 vmax = np.ceil(combined_data.max())
2298 else:
2299 vmin = min(x_levels)
2300 vmax = max(x_levels)
2302 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2303 # a single profile directly without iterating over a sequence.
2304 sequence_coords = [
2305 cube.coord(sequence_coordinate)
2306 for cube in cubes
2307 if cube.coords(sequence_coordinate)
2308 ]
2309 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2310 np.size(coord.points) > 1 for coord in sequence_coords
2311 )
2312 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2313 np.size(coord.points) == 1 for coord in sequence_coords
2314 )
2316 plot_index = []
2317 if has_sequence_coord: 2317 ↛ 2342line 2317 didn't jump to line 2342 because the condition on line 2317 was always true
2318 # Matching the slices (matching by seq coord point; it may happen that
2319 # evaluated models do not cover the same seq coord range, hence matching
2320 # necessary)
2321 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2322 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2323 for cubes_slice in cube_iterables:
2324 # Format the coordinate value in a unit appropriate way.
2325 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2326 plot_title, plot_filename = _set_title_and_filename(
2327 seq_coord, nplot, recipe_title, filename
2328 )
2330 # Do the actual plotting.
2331 _plot_and_save_vertical_line_series(
2332 cubes_slice,
2333 coords,
2334 "realization",
2335 plot_filename,
2336 series_coordinate,
2337 title=plot_title,
2338 vmin=vmin,
2339 vmax=vmax,
2340 )
2341 plot_index.append(plot_filename)
2342 elif has_scalar_sequence_coord:
2343 # Scalar sequence coordinate (typically aggregated time bounds):
2344 # make one plot and include sequence period in title/filename.
2345 plot_title, plot_filename = _set_title_and_filename(
2346 sequence_coords[0], 1, recipe_title, filename
2347 )
2349 _plot_and_save_vertical_line_series(
2350 cubes,
2351 coords,
2352 "realization",
2353 plot_filename,
2354 series_coordinate,
2355 title=plot_title,
2356 vmin=vmin,
2357 vmax=vmax,
2358 )
2359 plot_index.append(plot_filename)
2360 else:
2361 # 1D case: no sequence coordinate, plot a single profile.
2362 plot_title = recipe_title
2363 if filename:
2364 plot_filename = filename
2365 else:
2366 plot_filename = f"{slugify(plot_title)}.png"
2368 _plot_and_save_vertical_line_series(
2369 cubes,
2370 coords,
2371 "realization",
2372 plot_filename,
2373 series_coordinate,
2374 title=plot_title,
2375 vmin=vmin,
2376 vmax=vmax,
2377 )
2378 plot_index.append(plot_filename)
2380 # Add list of plots to plot metadata.
2381 complete_plot_index = _append_to_plot_index(plot_index)
2383 # Make a page to display the plots.
2384 _make_plot_html_page(complete_plot_index)
2386 return cubes
2389def qq_plot(
2390 cubes: iris.cube.CubeList,
2391 coordinates: list[str],
2392 percentiles: list[float],
2393 model_names: list[str],
2394 filename: str = None,
2395 one_to_one: bool = True,
2396 **kwargs,
2397) -> iris.cube.CubeList:
2398 """Plot a Quantile-Quantile plot between two models for common time points.
2400 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2401 collapsed within the operator over all specified coordinates such as
2402 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2403 ensemble members to ensure a 1D cube (array).
2405 Parameters
2406 ----------
2407 cubes: iris.cube.CubeList
2408 Two cubes of the same variable with different models.
2409 coordinate: list[str]
2410 The list of coordinates to collapse over. This list should be
2411 every coordinate within the cube to result in a 1D cube around
2412 the percentile coordinate.
2413 percent: list[float]
2414 A list of percentiles to appear in the plot.
2415 model_names: list[str]
2416 A list of model names to appear on the axis of the plot.
2417 filename: str, optional
2418 Filename of the plot to write.
2419 one_to_one: bool, optional
2420 If True a 1:1 line is plotted; if False it is not. Default is True.
2422 Raises
2423 ------
2424 ValueError
2425 When the cubes are not compatible.
2427 Notes
2428 -----
2429 The quantile-quantile plot is a variant on the scatter plot representing
2430 two datasets by their quantiles (percentiles) for common time points.
2431 This plot does not use a theoretical distribution to compare against, but
2432 compares percentiles of two datasets. This plot does
2433 not use all raw data points, but plots the selected percentiles (quantiles) of
2434 each variable instead for the two datasets, thereby normalising the data for a
2435 direct comparison between the selected percentiles of the two dataset distributions.
2437 Quantile-quantile plots are valuable for comparing against
2438 observations and other models. Identical percentiles between the variables
2439 will lie on the one-to-one line implying the values correspond well to each
2440 other. Where there is a deviation from the one-to-one line a range of
2441 possibilities exist depending on how and where the data is shifted (e.g.,
2442 Wilks 2011 [Wilks2011]_).
2444 For distributions above the one-to-one line the distribution is left-skewed;
2445 below is right-skewed. A distinct break implies a bimodal distribution, and
2446 closer values/values further apart at the tails imply poor representation of
2447 the extremes.
2449 References
2450 ----------
2451 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2452 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2453 """
2454 # Check cubes using same functionality as the difference operator.
2455 if len(cubes) != 2:
2456 raise ValueError("cubes should contain exactly 2 cubes.")
2457 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2458 other: Cube = cubes.extract_cube(
2459 iris.Constraint(
2460 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2461 )
2462 )
2464 # Get spatial coord names.
2465 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2466 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2468 # Ensure cubes to compare are on common differencing grid.
2469 # This is triggered if either
2470 # i) latitude and longitude shapes are not the same. Note grid points
2471 # are not compared directly as these can differ through rounding
2472 # errors.
2473 # ii) or variables are known to often sit on different grid staggering
2474 # in different models (e.g. cell center vs cell edge), as is the case
2475 # for UM and LFRic comparisons.
2476 # In future greater choice of regridding method might be applied depending
2477 # on variable type. Linear regridding can in general be appropriate for smooth
2478 # variables. Care should be taken with interpretation of differences
2479 # given this dependency on regridding.
2480 if (
2481 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2482 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2483 ) or (
2484 base.long_name
2485 in [
2486 "eastward_wind_at_10m",
2487 "northward_wind_at_10m",
2488 "northward_wind_at_cell_centres",
2489 "eastward_wind_at_cell_centres",
2490 "zonal_wind_at_pressure_levels",
2491 "meridional_wind_at_pressure_levels",
2492 "potential_vorticity_at_pressure_levels",
2493 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2494 ]
2495 ):
2496 logging.debug(
2497 "Linear regridding base cube to other grid to compute differences"
2498 )
2499 base = regrid_onto_cube(base, other, method="Linear")
2501 # Extract just common time points.
2502 base, other = _extract_common_time_points(base, other)
2504 # Equalise attributes so we can merge.
2505 fully_equalise_attributes([base, other])
2506 logging.debug("Base: %s\nOther: %s", base, other)
2508 # Collapse cubes.
2509 base = collapse(
2510 base,
2511 coordinate=coordinates,
2512 method="PERCENTILE",
2513 additional_percent=percentiles,
2514 )
2515 other = collapse(
2516 other,
2517 coordinate=coordinates,
2518 method="PERCENTILE",
2519 additional_percent=percentiles,
2520 )
2522 # Ensure we have a name for the plot file.
2523 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2524 title = f"{recipe_title}"
2526 if filename is None:
2527 filename = slugify(recipe_title)
2529 # Add file extension.
2530 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2532 # Do the actual plotting on a scatter plot
2533 _plot_and_save_scatter_plot(
2534 base, other, plot_filename, title, one_to_one, model_names
2535 )
2537 # Add list of plots to plot metadata.
2538 plot_index = _append_to_plot_index([plot_filename])
2540 # Make a page to display the plots.
2541 _make_plot_html_page(plot_index)
2543 return iris.cube.CubeList([base, other])
2546def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2547 """
2548 Plot a Hinton style triangle/scorecard plot.
2550 This plot type can be useful for summarising high level information, such as comparing
2551 how 'skillful' two models are when verified against observations for a variety of metrics,
2552 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2553 than customisable by the user as input arguments; many have been designed to automatically
2554 scale the plot depending on the number of x and y components.
2556 Parameters
2557 ----------
2558 change: np.ndarray
2559 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2560 size/direction.
2561 signif: np.ndarray
2562 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2563 xaxis_labels: list
2564 List of labels for the xaxis (must match the second dimension length of signif and change,
2565 along with magnitude if not None).
2566 yaxis_labels: list
2567 List of labels for the yaxis (must match the first dimension length of signif and change,
2568 along with magnitude if not None).
2569 magnitude: np.ndarray | None
2570 Optional 2D array, matching the shape of change, signif, which contains numerical values
2571 the user wishes to display under each respective triangle.
2573 Returns
2574 -------
2575 matplotlib axes object to either display or do further modifications to.
2576 """
2577 # Setup colors of triangles
2578 color_pos = "#7CAE00"
2579 color_neg = "#7B68EE"
2581 # Setup cell/text size ratios
2582 figsize = None
2583 cell_size_in = 0.35
2584 text_row_ratio = 0.25
2586 # Ensure arrays, and change to bool for sig.
2587 change = np.asarray(change)
2588 signif = np.asarray(signif).astype(bool)
2589 if magnitude is not None: 2589 ↛ 2590line 2589 didn't jump to line 2590 because the condition on line 2589 was never true
2590 magnitude = np.asarray(magnitude)
2592 # Get the number of x and y elements
2593 ny, nx = change.shape
2595 # Build non-uniform y coordinates
2596 tri_height = 1.0
2597 txt_height = text_row_ratio
2599 tri_y = []
2600 txt_y = []
2601 y_edges = [0.0]
2603 y = 0.0
2604 for _j in range(ny):
2605 tri_y.append(y + tri_height / 2)
2606 y += tri_height
2607 y_edges.append(y)
2609 if magnitude is not None: 2609 ↛ 2610line 2609 didn't jump to line 2610 because the condition on line 2609 was never true
2610 txt_y.append(y + txt_height / 2)
2611 y += txt_height
2612 y_edges.append(y)
2614 total_height = y
2616 # Dynamic figure size
2617 if figsize is None: 2617 ↛ 2622line 2617 didn't jump to line 2622 because the condition on line 2617 was always true
2618 width = nx * cell_size_in
2619 height = total_height * cell_size_in + 2
2620 figsize = (width, height)
2622 fig, ax = plt.subplots(figsize=figsize)
2624 # Setup axes and grid.
2625 ax.set_aspect("equal", adjustable="box")
2626 ax.set_xlim(-0.5, nx - 0.5)
2627 ax.set_ylim(0, total_height)
2629 ax.set_xticks(np.arange(nx))
2630 ax.set_xticklabels(xaxis_labels, rotation=90)
2632 ax.set_yticks(tri_y)
2633 ax.set_yticklabels(yaxis_labels)
2635 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2636 ax.set_yticks(y_edges, minor=True)
2638 ax.set_axisbelow(True)
2639 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2640 ax.grid(False, which="major")
2641 ax.tick_params(which="minor", length=0)
2643 ax.invert_yaxis()
2645 # Compute marker scaling (fixed overlap)
2646 fig.canvas.draw()
2648 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2649 width_in, height_in = bbox.width, bbox.height
2651 cell_w = (width_in * fig.dpi) / nx
2652 cell_h = (height_in * fig.dpi) / total_height
2653 cell_pixels = min(cell_w, cell_h)
2655 max_marker_size = (0.6 * cell_pixels) ** 2
2657 text_fontsize = cell_pixels * 0.15
2659 # Plot triangles + text
2660 for j in range(ny):
2661 for i in range(nx):
2662 val = change[j, i]
2663 if np.isnan(val): 2663 ↛ 2664line 2663 didn't jump to line 2664 because the condition on line 2663 was never true
2664 continue
2666 if abs(val) < 0.01: 2666 ↛ 2667line 2666 didn't jump to line 2667 because the condition on line 2666 was never true
2667 continue
2669 sig = signif[j, i]
2670 size = max_marker_size * abs(val)
2672 # Triangle style
2673 if val >= 0:
2674 marker = "^"
2675 color = color_pos
2676 else:
2677 marker = "v"
2678 color = color_neg
2680 if sig:
2681 edgecolor = "black"
2682 linewidth = 0.6
2683 else:
2684 edgecolor = "none"
2685 linewidth = 0.0
2687 # Triangle
2688 ax.scatter(
2689 i,
2690 tri_y[j],
2691 s=size,
2692 marker=marker,
2693 c=color,
2694 edgecolors=edgecolor,
2695 linewidths=linewidth,
2696 zorder=3,
2697 clip_on=True, # ensures no rendering bleed
2698 )
2700 # Text row
2701 if magnitude is not None: 2701 ↛ 2702line 2701 didn't jump to line 2702 because the condition on line 2701 was never true
2702 mag_val = magnitude[j, i]
2704 if not np.isnan(mag_val):
2705 ax.text(
2706 i,
2707 txt_y[j],
2708 f"{mag_val:.1f}",
2709 ha="center",
2710 va="center",
2711 fontsize=text_fontsize,
2712 color="black",
2713 zorder=4,
2714 )
2716 plt.tight_layout()
2717 return fig, ax
2720def scatter_plot(
2721 cube_x: iris.cube.Cube | iris.cube.CubeList,
2722 cube_y: iris.cube.Cube | iris.cube.CubeList,
2723 filename: str = None,
2724 one_to_one: bool = True,
2725 **kwargs,
2726) -> iris.cube.CubeList:
2727 """Plot a scatter plot between two variables.
2729 Both cubes must be 1D.
2731 Parameters
2732 ----------
2733 cube_x: Cube | CubeList
2734 1 dimensional Cube of the data to plot on y-axis.
2735 cube_y: Cube | CubeList
2736 1 dimensional Cube of the data to plot on x-axis.
2737 filename: str, optional
2738 Filename of the plot to write.
2739 one_to_one: bool, optional
2740 If True a 1:1 line is plotted; if False it is not. Default is True.
2742 Returns
2743 -------
2744 cubes: CubeList
2745 CubeList of the original x and y cubes for further processing.
2747 Raises
2748 ------
2749 ValueError
2750 If the cube doesn't have the right dimensions and cubes not the same
2751 size.
2752 TypeError
2753 If the cube isn't a single cube.
2755 Notes
2756 -----
2757 Scatter plots are used for determining if there is a relationship between
2758 two variables. Positive relations have a slope going from bottom left to top
2759 right; Negative relations have a slope going from top left to bottom right.
2760 """
2761 # Iterate over all cubes in cube or CubeList and plot.
2762 for cube_iter in iter_maybe(cube_x):
2763 # Check cubes are correct shape.
2764 cube_iter = check_single_cube(cube_iter)
2765 if cube_iter.ndim > 1:
2766 raise ValueError("cube_x must be 1D.")
2768 # Iterate over all cubes in cube or CubeList and plot.
2769 for cube_iter in iter_maybe(cube_y):
2770 # Check cubes are correct shape.
2771 cube_iter = check_single_cube(cube_iter)
2772 if cube_iter.ndim > 1:
2773 raise ValueError("cube_y must be 1D.")
2775 # Ensure we have a name for the plot file.
2776 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2777 title = f"{recipe_title}"
2779 if filename is None:
2780 filename = slugify(recipe_title)
2782 # Add file extension.
2783 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2785 # Do the actual plotting.
2786 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2788 # Add list of plots to plot metadata.
2789 plot_index = _append_to_plot_index([plot_filename])
2791 # Make a page to display the plots.
2792 _make_plot_html_page(plot_index)
2794 return iris.cube.CubeList([cube_x, cube_y])
2797def vector_plot(
2798 cube_u: iris.cube.Cube,
2799 cube_v: iris.cube.Cube,
2800 filename: str = None,
2801 sequence_coordinate: str = "time",
2802 **kwargs,
2803) -> iris.cube.CubeList:
2804 """Plot a vector plot based on the input u and v components."""
2805 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2807 # Cubes must have a matching sequence coordinate.
2808 try:
2809 # Check that the u and v cubes have the same sequence coordinate.
2810 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2810 ↛ anywhereline 2810 didn't jump anywhere: it always raised an exception.
2811 raise ValueError("Coordinates do not match.")
2812 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2813 raise ValueError(
2814 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2815 ) from err
2817 # Create a plot for each value of the sequence coordinate.
2818 plot_index = []
2819 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2820 for cube_u_slice, cube_v_slice in zip(
2821 cube_u.slices_over(sequence_coordinate),
2822 cube_v.slices_over(sequence_coordinate),
2823 strict=True,
2824 ):
2825 # Format the coordinate value in a unit appropriate way.
2826 seq_coord = cube_u_slice.coord(sequence_coordinate)
2827 plot_title, plot_filename = _set_title_and_filename(
2828 seq_coord, nplot, recipe_title, filename
2829 )
2831 # Do the actual plotting.
2832 _plot_and_save_vector_plot(
2833 cube_u_slice,
2834 cube_v_slice,
2835 filename=plot_filename,
2836 title=plot_title,
2837 method="pcolormesh",
2838 )
2839 plot_index.append(plot_filename)
2841 # Add list of plots to plot metadata.
2842 complete_plot_index = _append_to_plot_index(plot_index)
2844 # Make a page to display the plots.
2845 _make_plot_html_page(complete_plot_index)
2847 return iris.cube.CubeList([cube_u, cube_v])
2850def plot_histogram_series(
2851 cubes: iris.cube.Cube | iris.cube.CubeList,
2852 filename: str = None,
2853 sequence_coordinate: str = "time",
2854 stamp_coordinate: str = "realization",
2855 single_plot: bool = False,
2856 **kwargs,
2857) -> iris.cube.Cube | iris.cube.CubeList:
2858 """Plot a histogram plot for each vertical level provided.
2860 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2861 is present then a sequence of plots will be produced using the time slider
2862 functionality to scroll through histograms against time. If a
2863 stamp_coordinate is present then postage stamp plots will be produced. If
2864 stamp_coordinate and single_plot is True, all postage stamp plots will be
2865 plotted in a single plot instead of separate postage stamp plots.
2867 Parameters
2868 ----------
2869 cubes: Cube | iris.cube.CubeList
2870 Iris cube or CubeList of the data to plot. It should have a single dimension other
2871 than the stamp coordinate.
2872 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2873 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2874 filename: str, optional
2875 Name of the plot to write, used as a prefix for plot sequences. Defaults
2876 to the recipe name.
2877 sequence_coordinate: str, optional
2878 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2879 This coordinate must exist in the cube and will be used for the time
2880 slider.
2881 stamp_coordinate: str, optional
2882 Coordinate about which to plot postage stamp plots. Defaults to
2883 ``"realization"``.
2884 single_plot: bool, optional
2885 If True, all postage stamp plots will be plotted in a single plot. If
2886 False, each postage stamp plot will be plotted separately. Is only valid
2887 if stamp_coordinate exists and has more than a single point.
2889 Returns
2890 -------
2891 iris.cube.Cube | iris.cube.CubeList
2892 The original Cube or CubeList (so further operations can be applied).
2893 Plotted data.
2895 Raises
2896 ------
2897 ValueError
2898 If the cube doesn't have the right dimensions.
2899 TypeError
2900 If the cube isn't a Cube or CubeList.
2901 """
2902 recipe_title = get_recipe_metadata().get("title", "Histogram")
2904 cubes = iter_maybe(cubes)
2905 # Ensure we have a name for the plot file.
2906 if filename is None:
2907 filename = slugify(recipe_title)
2909 # Internal plotting function.
2910 plotting_func = _plot_and_save_histogram_series
2912 num_models = get_num_models(cubes)
2914 validate_cube_shape(cubes, num_models)
2916 # If several histograms are plotted, check sequence_coordinate
2917 check_sequence_coordinate(cubes, sequence_coordinate)
2919 # Get axis minimum and maximum from levels information.
2920 # If no levels set, derive minima and maxima from data in CubeList.
2921 vmin, vmax = _set_axis_range(cubes)
2923 # Make postage stamp plots if stamp_coordinate exists and has more than a
2924 # single point. If single_plot is True:
2925 # -- all postage stamp plots will be plotted in a single plot instead of
2926 # separate postage stamp plots.
2927 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2928 # produced per single model only
2929 if num_models == 1:
2930 if ( 2930 ↛ 2934line 2930 didn't jump to line 2934 because the condition on line 2930 was never true
2931 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2932 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2933 ):
2934 if single_plot:
2935 plotting_func = (
2936 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2937 )
2938 else:
2939 plotting_func = _plot_and_save_postage_stamp_histogram_series
2940 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2941 else:
2942 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2944 plot_index = []
2945 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2946 # Create a plot for each value of the sequence coordinate. Allowing for
2947 # multiple cubes in a CubeList to be plotted in the same plot for similar
2948 # sequence values. Passing a CubeList into the internal plotting function
2949 # for similar values of the sequence coordinate. cube_slice can be an
2950 # iris.cube.Cube or an iris.cube.CubeList.
2951 for cube_slice in cube_iterables:
2952 single_cube = cube_slice
2953 if isinstance(cube_slice, iris.cube.CubeList):
2954 single_cube = cube_slice[0]
2956 # Ensure valid stamp coordinate in cube dimensions
2957 if stamp_coordinate == "realization": 2957 ↛ 2960line 2957 didn't jump to line 2960 because the condition on line 2957 was always true
2958 stamp_coordinate = check_stamp_coordinate(single_cube)
2959 # Set plot titles and filename, based on sequence coordinate
2960 seq_coord = single_cube.coord(sequence_coordinate)
2961 # Use time coordinate in title and filename if single histogram output.
2962 if sequence_coordinate == "realization" and nplot == 1: 2962 ↛ 2963line 2962 didn't jump to line 2963 because the condition on line 2962 was never true
2963 seq_coord = single_cube.coord("time")
2964 plot_title, plot_filename = _set_title_and_filename(
2965 seq_coord, nplot, recipe_title, filename
2966 )
2968 # Do the actual plotting.
2969 plotting_func(
2970 cube_slice,
2971 filename=plot_filename,
2972 stamp_coordinate=stamp_coordinate,
2973 title=plot_title,
2974 vmin=vmin,
2975 vmax=vmax,
2976 )
2977 plot_index.append(plot_filename)
2979 # Add list of plots to plot metadata.
2980 complete_plot_index = _append_to_plot_index(plot_index)
2982 # Make a page to display the plots.
2983 _make_plot_html_page(complete_plot_index)
2985 return cubes
2988def _plot_and_save_postage_stamp_power_spectrum_series(
2989 cubes: iris.cube.Cube,
2990 coords: list[iris.coords.Coord],
2991 stamp_coordinate: str,
2992 filename: str,
2993 title: str,
2994 series_coordinate: str = None,
2995 **kwargs,
2996):
2997 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2999 Parameters
3000 ----------
3001 cubes: Cube or CubeList
3002 Cube or Cubelist of the power spectrum data.
3003 coords: list[Coord]
3004 Coordinates to plot on the x-axis, one per cube.
3005 stamp_coordinate: str
3006 Coordinate that becomes different plots.
3007 filename: str
3008 Filename of the plot to write.
3009 title: str
3010 Plot title.
3011 series_coordinate: str, optional
3012 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3014 """
3015 # Use the smallest square grid that will fit the members.
3016 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
3018 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3019 model_colors_map = get_model_colors_map(cubes)
3020 # ax = plt.gca()
3021 # Make a subplot for each member.
3022 for member, subplot in zip(
3023 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3024 ):
3025 ax = plt.subplot(grid_size, grid_size, subplot)
3027 # Store min/max ranges.
3028 y_levels = []
3030 line_marker = None
3031 line_width = 1
3033 for cube in iter_maybe(member):
3034 xcoord = _select_series_coord(cube, series_coordinate)
3035 xname = xcoord.points
3037 yfield = cube.data # power spectrum
3038 label = None
3039 color = "black"
3040 if model_colors_map: 3040 ↛ 3041line 3040 didn't jump to line 3041 because the condition on line 3040 was never true
3041 label = cube.attributes.get("model_name")
3042 color = model_colors_map.get(label)
3044 if member.coord(stamp_coordinate).points == [0]:
3045 ax.plot(
3046 xname,
3047 yfield,
3048 color=color,
3049 marker=line_marker,
3050 ls="-",
3051 lw=line_width,
3052 label=f"{label} (control)"
3053 if len(cube.coord(stamp_coordinate).points) > 1
3054 else label,
3055 )
3056 # Label with member if part of an ensemble and not the control.
3057 else:
3058 ax.plot(
3059 xname,
3060 yfield,
3061 color=color,
3062 ls="-",
3063 lw=1.5,
3064 alpha=0.75,
3065 label=f"{label} (member)",
3066 )
3068 # Calculate the global min/max if multiple cubes are given.
3069 _, levels, _ = colorbar_map_levels(cube, axis="y")
3070 if levels is not None: 3070 ↛ 3071line 3070 didn't jump to line 3071 because the condition on line 3070 was never true
3071 y_levels.append(min(levels))
3072 y_levels.append(max(levels))
3074 # Add some labels and tweak the style.
3075 title = f"{title}"
3076 ax.set_title(title, fontsize=16)
3078 # Set appropriate x-axis label based on coordinate
3079 if series_coordinate == "wavelength" or ( 3079 ↛ 3082line 3079 didn't jump to line 3082 because the condition on line 3079 was never true
3080 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3081 ):
3082 ax.set_xlabel("Wavelength (km)", fontsize=14)
3083 elif series_coordinate == "physical_wavenumber" or ( 3083 ↛ 3088line 3083 didn't jump to line 3088 because the condition on line 3083 was always true
3084 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3085 ):
3086 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3087 else: # frequency or check units
3088 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3089 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3090 else:
3091 ax.set_xlabel("Wavenumber", fontsize=14)
3093 ax.set_ylabel("Power Spectral Density", fontsize=14)
3094 ax.tick_params(axis="both", labelsize=12)
3096 # Set log-log scale
3097 ax.set_xscale("log")
3098 ax.set_yscale("log")
3100 # Add gridlines
3101 ax.grid(linestyle="--", color="grey", linewidth=1)
3102 # Ientify unique labels for legend
3103 handles = list(
3104 {
3105 label: handle
3106 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3107 }.values()
3108 )
3109 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3111 ax = plt.gca()
3112 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3114 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3115 logging.info("Saved histogram postage stamp plot to %s", filename)
3116 plt.close(fig)
3119def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3120 cubes: iris.cube.Cube,
3121 coords: list[iris.coords.Coord],
3122 stamp_coordinate: str,
3123 filename: str,
3124 title: str,
3125 series_coordinate: str = None,
3126 **kwargs,
3127):
3128 """Plot and save power spectra for ensemble members in single plot.
3130 Parameters
3131 ----------
3132 cubes: Cube or CubeList
3133 Cube or Cubelist of the power spectrum data.
3134 coords: list[Coord]
3135 Coordinates to plot on the x-axis, one per cube.
3136 stamp_coordinate: str
3137 Coordinate that becomes different plots.
3138 filename: str
3139 Filename of the plot to write.
3140 title: str
3141 Plot title.
3142 series_coordinate: str, optional
3143 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3145 """
3146 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3147 model_colors_map = get_model_colors_map(cubes)
3149 line_marker = None
3150 line_width = 1
3152 # Compute ensemble statistics to show spread
3153 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3154 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3155 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3157 xcoord_global = mean_cube.coord(series_coordinate)
3158 x_global = xcoord_global.points
3160 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3161 xcoord = _select_series_coord(member, series_coordinate)
3162 xname = xcoord.points
3164 yfield = member.data # power spectrum
3165 color = "black"
3166 if model_colors_map: 3166 ↛ 3170line 3166 didn't jump to line 3170 because the condition on line 3166 was always true
3167 label = member.attributes.get("model_name") if i == 0 else None
3168 color = model_colors_map.get(label)
3170 if member.coord(stamp_coordinate).points == [0]:
3171 ax.plot(
3172 xname,
3173 yfield,
3174 color=color,
3175 marker=line_marker,
3176 ls="-",
3177 lw=line_width,
3178 label=f"{label} (control)"
3179 if len(member.coord(stamp_coordinate).points) > 1
3180 else label,
3181 )
3182 # Label with member number if part of an ensemble and not the control.
3183 else:
3184 ax.plot(
3185 xname,
3186 yfield,
3187 color=color,
3188 ls="-",
3189 lw=1.5,
3190 alpha=0.75,
3191 label=label,
3192 )
3194 # Set appropriate x-axis label based on coordinate
3195 if series_coordinate == "wavelength" or ( 3195 ↛ 3198line 3195 didn't jump to line 3198 because the condition on line 3195 was never true
3196 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3197 ):
3198 ax.set_xlabel("Wavelength (km)", fontsize=14)
3199 elif series_coordinate == "physical_wavenumber" or ( 3199 ↛ 3204line 3199 didn't jump to line 3204 because the condition on line 3199 was always true
3200 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3201 ):
3202 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3203 else: # frequency or check units
3204 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3205 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3206 else:
3207 ax.set_xlabel("Wavenumber", fontsize=14)
3209 # Add ensemble spread shading
3210 ax.fill_between(
3211 x_global,
3212 min_cube.data,
3213 max_cube.data,
3214 color="grey",
3215 alpha=0.3,
3216 label="Ensemble spread",
3217 )
3219 # Add ensemble mean line
3220 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3222 ax.set_ylabel("Power Spectral Density", fontsize=14)
3223 ax.tick_params(axis="both", labelsize=12)
3225 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3226 # Set log-log scale
3227 ax.set_xscale("log")
3228 ax.set_yscale("log")
3230 # Add gridlines
3231 ax.grid(linestyle="--", color="grey", linewidth=1)
3232 # Identify unique labels for legend
3233 handles = list(
3234 {
3235 label: handle
3236 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3237 }.values()
3238 )
3239 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3241 # Figure title.
3242 ax.set_title(title, fontsize=16)
3244 # Save the figure to a file
3245 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3247 # Close the figure
3248 plt.close(fig)