Coverage for src / CSET / operators / plot.py: 81%
1110 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-20 13:52 +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 functools
19import importlib.resources
20import itertools
21import json
22import logging
23import math
24import os
25from typing import Literal
27import cartopy.crs as ccrs
28import cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.colors as mcolors
36import matplotlib.pyplot as plt
37import numpy as np
38from cartopy.mpl.geoaxes import GeoAxes
39from iris.cube import Cube
40from markdown_it import MarkdownIt
41from mpl_toolkits.axes_grid1.inset_locator import inset_axes
43from CSET._common import (
44 combine_dicts,
45 filename_slugify,
46 get_recipe_metadata,
47 iter_maybe,
48 render_file,
49 slugify,
50)
51from CSET.operators._utils import (
52 fully_equalise_attributes,
53 get_cube_yxcoordname,
54 is_transect,
55 slice_over_maybe,
56)
57from CSET.operators.collapse import collapse
58from CSET.operators.misc import _extract_common_time_points
59from CSET.operators.regrid import regrid_onto_cube
61# Suppress matplotlib debug output
62# logging.getLogger('matplotlib.ticker').setLevel(logging.WARNING)
64# Use a non-interactive plotting backend.
65mpl.use("agg")
67DEFAULT_DISCRETE_COLORS = mpl.colormaps["tab10"].colors + mpl.colormaps["Accent"].colors
69############################
70# Private helper functions #
71############################
74def _append_to_plot_index(plot_index: list) -> list:
75 """Add plots into the plot index, returning the complete plot index."""
76 with open("meta.json", "r+t", encoding="UTF-8") as fp:
77 fcntl.flock(fp, fcntl.LOCK_EX)
78 fp.seek(0)
79 meta = json.load(fp)
80 complete_plot_index = meta.get("plots", [])
81 complete_plot_index = complete_plot_index + plot_index
82 meta["plots"] = complete_plot_index
83 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
84 os.getenv("DO_CASE_AGGREGATION")
85 ):
86 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
87 fp.seek(0)
88 fp.truncate()
89 json.dump(meta, fp, indent=2)
90 return complete_plot_index
93def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:
94 """Ensure a single cube is given.
96 If a CubeList of length one is given that the contained cube is returned,
97 otherwise an error is raised.
99 Parameters
100 ----------
101 cube: Cube | CubeList
102 The cube to check.
104 Returns
105 -------
106 cube: Cube
107 The checked cube.
109 Raises
110 ------
111 TypeError
112 If the input cube is not a Cube or CubeList of a single Cube.
113 """
114 if isinstance(cube, iris.cube.Cube):
115 return cube
116 if isinstance(cube, iris.cube.CubeList):
117 if len(cube) == 1:
118 return cube[0]
119 raise TypeError("Must have a single cube", cube)
122def _make_plot_html_page(plots: list):
123 """Create a HTML page to display a plot image."""
124 # Debug check that plots actually contains some strings.
125 assert isinstance(plots[0], str)
127 # Load HTML template file.
128 operator_files = importlib.resources.files()
129 template_file = operator_files.joinpath("_plot_page_template.html")
131 # Get some metadata.
132 meta = get_recipe_metadata()
133 title = meta.get("title", "Untitled")
134 description = MarkdownIt().render(meta.get("description", "*No description.*"))
136 # Prepare template variables.
137 variables = {
138 "title": title,
139 "description": description,
140 "initial_plot": plots[0],
141 "plots": plots,
142 "title_slug": slugify(title),
143 }
145 # Render template.
146 html = render_file(template_file, **variables)
148 # Save completed HTML.
149 with open("index.html", "wt", encoding="UTF-8") as fp:
150 fp.write(html)
153@functools.cache
154def _load_colorbar_map(user_colorbar_file: str = None) -> dict:
155 """Load the colorbar definitions from a file.
157 This is a separate function to make it cacheable.
158 """
159 colorbar_file = importlib.resources.files().joinpath("_colorbar_definition.json")
160 with open(colorbar_file, "rt", encoding="UTF-8") as fp:
161 colorbar = json.load(fp)
163 logging.debug("User colour bar file: %s", user_colorbar_file)
164 override_colorbar = {}
165 if user_colorbar_file:
166 try:
167 with open(user_colorbar_file, "rt", encoding="UTF-8") as fp:
168 override_colorbar = json.load(fp)
169 except FileNotFoundError:
170 logging.warning("Colorbar file does not exist. Using default values.")
172 # Overwrite values with the user supplied colorbar definition.
173 colorbar = combine_dicts(colorbar, override_colorbar)
174 return colorbar
177def _get_model_colors_map(cubes: iris.cube.CubeList | iris.cube.Cube) -> dict:
178 """Get an appropriate colors for model lines in line plots.
180 For each model in the list of cubes colors either from user provided
181 color definition file (so-called style file) or from default colors are mapped
182 to model_name attribute.
184 Parameters
185 ----------
186 cubes: CubeList or Cube
187 Cubes with model_name attribute
189 Returns
190 -------
191 model_colors_map:
192 Dictionary mapping model_name attribute to colors
193 """
194 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
195 colorbar = _load_colorbar_map(user_colorbar_file)
196 model_names = sorted(
197 filter(
198 lambda x: x is not None,
199 (cube.attributes.get("model_name", None) for cube in iter_maybe(cubes)),
200 )
201 )
202 if not model_names:
203 return {}
204 use_user_colors = all(mname in colorbar.keys() for mname in model_names)
205 if use_user_colors: 205 ↛ 206line 205 didn't jump to line 206 because the condition on line 205 was never true
206 return {mname: colorbar[mname] for mname in model_names}
208 color_list = itertools.cycle(DEFAULT_DISCRETE_COLORS)
209 return {mname: color for mname, color in zip(model_names, color_list, strict=False)}
212def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
213 """Get an appropriate colorbar for the given cube.
215 For the given variable the appropriate colorbar is looked up from a
216 combination of the built-in CSET colorbar definitions, and any user supplied
217 definitions. As well as varying on variables, these definitions may also
218 exist for specific pressure levels to account for variables with
219 significantly different ranges at different heights. The colorbars also exist
220 for masks and mask differences for considering variable presence diagnostics.
221 Specific variable ranges can be separately set in user-supplied definition
222 for x- or y-axis limits, or indicate where automated range preferred.
224 Parameters
225 ----------
226 cube: Cube
227 Cube of variable for which the colorbar information is desired.
228 axis: "x", "y", optional
229 Select the levels for just this axis of a line plot. The min and max
230 can be set by xmin/xmax or ymin/ymax respectively. For variables where
231 setting a universal range is not desirable (e.g. temperature), users
232 can set ymin/ymax values to "auto" in the colorbar definitions file.
233 Where no additional xmin/xmax or ymin/ymax values are provided, the
234 axis bounds default to use the vmin/vmax values provided.
236 Returns
237 -------
238 cmap:
239 Matplotlib colormap.
240 levels:
241 List of levels to use for plotting. For continuous plots the min and max
242 should be taken as the range.
243 norm:
244 BoundaryNorm information.
245 """
246 # Grab the colorbar file from the recipe global metadata.
247 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
248 colorbar = _load_colorbar_map(user_colorbar_file)
249 cmap = None
251 try:
252 # We assume that pressure is a scalar coordinate here.
253 pressure_level_raw = cube.coord("pressure").points[0]
254 # Ensure pressure_level is a string, as it is used as a JSON key.
255 pressure_level = str(int(pressure_level_raw))
256 except iris.exceptions.CoordinateNotFoundError:
257 pressure_level = None
259 # First try long name, then standard name, then var name. This order is used
260 # as long name is the one we correct between models, so it most likely to be
261 # consistent.
262 varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name]))
263 for varname in varnames:
264 # Get the colormap for this variable.
265 try:
266 var_colorbar = colorbar[varname]
267 cmap = plt.get_cmap(colorbar[varname]["cmap"], 51)
268 varname_key = varname
269 break
270 except KeyError:
271 logging.debug("Cube name %s has no colorbar definition.", varname)
273 # Get colormap if it is a mask.
274 if any("mask_for_" in name for name in varnames):
275 cmap, levels, norm = _custom_colormap_mask(cube, axis=axis)
276 return cmap, levels, norm
277 # If winds on Beaufort Scale use custom colorbar and levels
278 if any("Beaufort_Scale" in name for name in varnames):
279 cmap, levels, norm = _custom_beaufort_scale(cube, axis=axis)
280 return cmap, levels, norm
281 # If probability is plotted use custom colorbar and levels
282 if any("probability_of_" in name for name in varnames):
283 cmap, levels, norm = _custom_colormap_probability(cube, axis=axis)
284 return cmap, levels, norm
285 # If aviation colour state use custom colorbar and levels
286 if any("aviation_colour_state" in name for name in varnames):
287 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)
288 return cmap, levels, norm
290 # If no valid colormap has been defined, use defaults and return.
291 if not cmap:
292 logging.warning("No colorbar definition exists for %s.", cube.name())
293 cmap, levels, norm = mpl.colormaps["viridis"], None, None
294 return cmap, levels, norm
296 # Test if pressure-level specific settings are provided for cube.
297 if pressure_level:
298 try:
299 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]
300 except KeyError:
301 logging.debug(
302 "%s has no colorbar definition for pressure level %s.",
303 varname,
304 pressure_level,
305 )
307 # Check for availability of x-axis or y-axis user-specific overrides
308 # for setting level bounds for line plot types and return just levels.
309 # Line plots do not need a colormap, and just use the data range.
310 if axis:
311 if axis == "x":
312 try:
313 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]
314 except KeyError:
315 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
316 if axis == "y":
317 try:
318 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]
319 except KeyError:
320 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
321 # Check if user-specified auto-scaling for this variable
322 if vmin == "auto" or vmax == "auto":
323 levels = None
324 else:
325 levels = [vmin, vmax]
326 return None, levels, None
327 # Get and use the colorbar levels for this variable if spatial or histogram.
328 else:
329 try:
330 levels = var_colorbar["levels"]
331 # Use discrete bins when levels are specified, rather
332 # than a smooth range.
333 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
334 logging.debug("Using levels for %s colorbar.", varname)
335 logging.info("Using levels: %s", levels)
336 except KeyError:
337 # Get the range for this variable.
338 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
339 logging.debug("Using min and max for %s colorbar.", varname)
340 # Calculate levels from range.
341 levels = np.linspace(vmin, vmax, 101)
342 norm = None
344 # Overwrite cmap, levels and norm for specific variables that
345 # require custom colorbar_map as these can not be defined in the
346 # JSON file.
347 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)
348 cmap, levels, norm = _custom_colourmap_visibility_in_air(
349 cube, cmap, levels, norm
350 )
351 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)
352 return cmap, levels, norm
355def _setup_spatial_map(
356 cube: iris.cube.Cube,
357 figure,
358 cmap,
359 grid_size: int | None = None,
360 subplot: int | None = None,
361):
362 """Define map projections, extent and add coastlines and borderlines for spatial plots.
364 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
365 is specified, and map extent defined based on the input data.
367 Parameters
368 ----------
369 cube: Cube
370 2 dimensional (lat and lon) Cube of the data to plot.
371 figure:
372 Matplotlib Figure object holding all plot elements.
373 cmap:
374 Matplotlib colormap.
375 grid_size: int, optional
376 Size of grid for subplots if multiple spatial subplots in figure.
377 subplot: int, optional
378 Subplot index if multiple spatial subplots in figure.
380 Returns
381 -------
382 axes:
383 Matplotlib GeoAxes definition.
384 """
385 # Identify min/max plot bounds.
386 try:
387 lat_axis, lon_axis = get_cube_yxcoordname(cube)
388 x1 = np.min(cube.coord(lon_axis).points)
389 x2 = np.max(cube.coord(lon_axis).points)
390 y1 = np.min(cube.coord(lat_axis).points)
391 y2 = np.max(cube.coord(lat_axis).points)
393 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
394 if np.abs(x2 - x1) > 180.0:
395 x1 = x1 - 180.0
396 x2 = x2 - 180.0
397 logging.debug("Adjusting plot bounds to fit global extent.")
399 # Consider map projection orientation.
400 # Adapting orientation enables plotting across international dateline.
401 # Users can adapt the default central_longitude if alternative projections views.
402 if x2 > 180.0 or x1 < -180.0:
403 central_longitude = 180.0
404 else:
405 central_longitude = 0.0
407 # Define spatial map projection.
408 coord_system = cube.coord(lat_axis).coord_system
409 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
410 # Define rotated pole map projection for rotated pole inputs.
411 projection = ccrs.RotatedPole(
412 pole_longitude=coord_system.grid_north_pole_longitude,
413 pole_latitude=coord_system.grid_north_pole_latitude,
414 central_rotated_longitude=central_longitude,
415 )
416 crs = projection
417 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 417 ↛ 419line 417 didn't jump to line 419 because the condition on line 417 was never true
418 # Define Transverse Mercator projection for TM inputs.
419 projection = ccrs.TransverseMercator(
420 central_longitude=coord_system.longitude_of_central_meridian,
421 central_latitude=coord_system.latitude_of_projection_origin,
422 false_easting=coord_system.false_easting,
423 false_northing=coord_system.false_northing,
424 scale_factor=coord_system.scale_factor_at_central_meridian,
425 )
426 crs = projection
427 else:
428 # Define regular map projection for non-rotated pole inputs.
429 # Alternatives might include e.g. for global model outputs:
430 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
431 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
432 projection = ccrs.PlateCarree(central_longitude=central_longitude)
433 crs = ccrs.PlateCarree()
435 # Define axes for plot (or subplot) with required map projection.
436 if subplot is not None:
437 axes = figure.add_subplot(
438 grid_size, grid_size, subplot, projection=projection
439 )
440 else:
441 axes = figure.add_subplot(projection=projection)
443 # Add coastlines and borderlines if cube contains x and y map coordinates.
444 if cmap.name in ["viridis", "Greys"]:
445 coastcol = "magenta"
446 else:
447 coastcol = "black"
448 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
449 axes.coastlines(resolution="10m", color=coastcol)
450 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
452 # Add gridlines.
453 if subplot is None:
454 draw_labels = True
455 else:
456 draw_labels = False
457 gl = axes.gridlines(
458 alpha=0.3,
459 draw_labels=draw_labels,
460 dms=False,
461 x_inline=False,
462 y_inline=False,
463 )
464 gl.top_labels = False
465 gl.right_labels = False
467 # If is lat/lon spatial map, fix extent to keep plot tight.
468 # Specifying crs within set_extent helps ensure only data region is shown.
469 if isinstance(coord_system, iris.coord_systems.GeogCS):
470 axes.set_extent([x1, x2, y1, y2], crs=crs)
472 except ValueError:
473 # Skip if not both x and y map coordinates.
474 axes = figure.gca()
475 pass
477 return axes
480def _get_plot_resolution() -> int:
481 """Get resolution of rasterised plots in pixels per inch."""
482 return get_recipe_metadata().get("plot_resolution", 100)
485def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
486 """Return title and filename based on start and end points or bounds."""
487 if use_bounds and seq_coord.has_bounds():
488 vals = seq_coord.bounds.flatten()
489 else:
490 vals = seq_coord.points
491 start = seq_coord.units.title(vals[0])
492 end = seq_coord.units.title(vals[-1])
494 if start == end:
495 sequence_title = f"\n [{start}]"
496 sequence_fname = f"_{filename_slugify(start)}"
497 else:
498 sequence_title = f"\n [{start} to {end}]"
499 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
501 return sequence_title, sequence_fname
504def _set_title_and_filename(
505 seq_coord: iris.coords.Coord,
506 nplot: int,
507 recipe_title: str,
508 filename: str,
509):
510 """Set plot title and filename based on cube coordinate.
512 Parameters
513 ----------
514 sequence_coordinate: iris.coords.Coord
515 Coordinate about which to make a plot sequence.
516 nplot: int
517 Number of output plots to generate - controls title/naming.
518 recipe_title: str
519 Default plot title, potentially to update.
520 filename: str
521 Input plot filename, potentially to update.
523 Returns
524 -------
525 plot_title: str
526 Output formatted plot title string, based on plotted data.
527 plot_filename: str
528 Output formatted plot filename string.
529 """
530 ndim = seq_coord.ndim
531 npoints = np.size(seq_coord.points)
532 sequence_title = ""
533 sequence_fname = ""
535 # Case 1: Multiple dimension sequence input - list number of aggregated cases
536 # (e.g. aggregation histogram plots)
537 if ndim > 1:
538 ncase = np.shape(seq_coord)[0]
539 sequence_title = f"\n [{ncase} cases]"
540 sequence_fname = f"_{ncase}cases"
542 # Case 2: Single dimension input
543 else:
544 # Single sequence point
545 if npoints == 1:
546 if nplot > 1:
547 # Default labels for sequence inputs
548 sequence_value = seq_coord.units.title(seq_coord.points[0])
549 sequence_title = f"\n [{sequence_value}]"
550 sequence_fname = f"_{filename_slugify(sequence_value)}"
551 else:
552 # Aggregated attribute available where input collapsed over aggregation
553 try:
554 ncase = seq_coord.attributes["number_reference_times"]
555 sequence_title = f"\n [{ncase} cases]"
556 sequence_fname = f"_{ncase}cases"
557 except KeyError:
558 sequence_title, sequence_fname = _get_start_end_strings(
559 seq_coord, use_bounds=seq_coord.has_bounds()
560 )
561 # Multiple sequence (e.g. time) points
562 else:
563 sequence_title, sequence_fname = _get_start_end_strings(
564 seq_coord, use_bounds=False
565 )
567 # Set plot title and filename
568 plot_title = f"{recipe_title}{sequence_title}"
570 # Set plot filename, defaulting to user input if provided.
571 if filename is None:
572 filename = slugify(recipe_title)
573 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
574 else:
575 if nplot > 1:
576 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
577 else:
578 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
580 return plot_title, plot_filename
583def select_series_coord(cube, series_coordinate):
584 """Determine the grid coordinates to use to calculate grid spacing."""
585 try:
586 # Try the requested coordinate first
587 return cube.coord(series_coordinate)
589 except iris.exceptions.CoordinateNotFoundError:
590 # Fallback logic
591 if series_coordinate == "frequency":
592 fallbacks = ("physical_wavenumber", "wavelength")
594 elif series_coordinate == "physical_wavenumber":
595 fallbacks = ("frequency", "wavelength")
597 elif series_coordinate == "wavelength": 597 ↛ 602line 597 didn't jump to line 602 because the condition on line 597 was always true
598 fallbacks = ("frequency", "physical_wavenumber")
600 else:
601 # Unknown coordinate type -> re-raise the original exception
602 raise
604 # Try the fallbacks in order
605 for fb in fallbacks: 605 ↛ 613line 605 didn't jump to line 613 because the loop on line 605 didn't complete
606 try:
607 return cube.coord(fb)
608 except iris.exceptions.CoordinateNotFoundError:
609 continue
611 # If we get here, none of the fallback options were found
613 raise iris.exceptions.CoordinateNotFoundError(
614 f"No valid coordinate found for '{series_coordinate}' "
615 f"or fallback options {fallbacks}"
616 ) from None
619def _plot_and_save_spatial_plot(
620 cube: iris.cube.Cube,
621 filename: str,
622 title: str,
623 method: Literal["contourf", "pcolormesh"],
624 overlay_cube: iris.cube.Cube | None = None,
625 contour_cube: iris.cube.Cube | None = None,
626 **kwargs,
627):
628 """Plot and save a spatial plot.
630 Parameters
631 ----------
632 cube: Cube
633 2 dimensional (lat and lon) Cube of the data to plot.
634 filename: str
635 Filename of the plot to write.
636 title: str
637 Plot title.
638 method: "contourf" | "pcolormesh"
639 The plotting method to use.
640 overlay_cube: Cube, optional
641 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
642 contour_cube: Cube, optional
643 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
644 """
645 # Setup plot details, size, resolution, etc.
646 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
648 # Specify the color bar
649 cmap, levels, norm = _colorbar_map_levels(cube)
651 # If overplotting, set required colorbars
652 if overlay_cube:
653 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
654 if contour_cube:
655 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
657 # Setup plot map projection, extent and coastlines and borderlines.
658 axes = _setup_spatial_map(cube, fig, cmap)
660 # Plot the field.
661 if method == "contourf":
662 # Filled contour plot of the field.
663 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
664 elif method == "pcolormesh":
665 try:
666 vmin = min(levels)
667 vmax = max(levels)
668 except TypeError:
669 vmin, vmax = None, None
670 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
671 # if levels are defined.
672 if norm is not None:
673 vmin = None
674 vmax = None
675 logging.debug("Plotting using defined levels.")
676 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
677 else:
678 raise ValueError(f"Unknown plotting method: {method}")
680 # Overplot overlay field, if required
681 if overlay_cube:
682 try:
683 over_vmin = min(over_levels)
684 over_vmax = max(over_levels)
685 except TypeError:
686 over_vmin, over_vmax = None, None
687 if over_norm is not None: 687 ↛ 688line 687 didn't jump to line 688 because the condition on line 687 was never true
688 over_vmin = None
689 over_vmax = None
690 overlay = iplt.pcolormesh(
691 overlay_cube,
692 cmap=over_cmap,
693 norm=over_norm,
694 alpha=0.8,
695 vmin=over_vmin,
696 vmax=over_vmax,
697 )
698 # Overplot contour field, if required, with contour labelling.
699 if contour_cube:
700 contour = iplt.contour(
701 contour_cube,
702 colors="darkgray",
703 levels=cntr_levels,
704 norm=cntr_norm,
705 alpha=0.5,
706 linestyles="--",
707 linewidths=1,
708 )
709 plt.clabel(contour)
711 # Check to see if transect, and if so, adjust y axis.
712 if is_transect(cube):
713 if "pressure" in [coord.name() for coord in cube.coords()]:
714 axes.invert_yaxis()
715 axes.set_yscale("log")
716 axes.set_ylim(1100, 100)
717 # If both model_level_number and level_height exists, iplt can construct
718 # plot as a function of height above orography (NOT sea level).
719 elif {"model_level_number", "level_height"}.issubset( 719 ↛ 724line 719 didn't jump to line 724 because the condition on line 719 was always true
720 {coord.name() for coord in cube.coords()}
721 ):
722 axes.set_yscale("log")
724 axes.set_title(
725 f"{title}\n"
726 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
727 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
728 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
729 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
730 fontsize=16,
731 )
733 # Inset code
734 axins = inset_axes(
735 axes,
736 width="20%",
737 height="20%",
738 loc="upper right",
739 axes_class=GeoAxes,
740 axes_kwargs={"map_projection": ccrs.PlateCarree()},
741 )
743 # Slightly transparent to reduce plot blocking.
744 axins.patch.set_alpha(0.4)
746 axins.coastlines(resolution="50m")
747 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
749 SLat, SLon, ELat, ELon = (
750 float(coord) for coord in cube.attributes["transect_coords"].split("_")
751 )
753 # Draw line between them
754 axins.plot(
755 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
756 )
758 # Plot points (note: lon, lat order for Cartopy)
759 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
760 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
762 lon_min, lon_max = sorted([SLon, ELon])
763 lat_min, lat_max = sorted([SLat, ELat])
765 # Midpoints
766 lon_mid = (lon_min + lon_max) / 2
767 lat_mid = (lat_min + lat_max) / 2
769 # Maximum half-range
770 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
771 if half_range == 0: # points identical → provide small default 771 ↛ 775line 771 didn't jump to line 775 because the condition on line 771 was always true
772 half_range = 1
774 # Set square extent
775 axins.set_extent(
776 [
777 lon_mid - half_range,
778 lon_mid + half_range,
779 lat_mid - half_range,
780 lat_mid + half_range,
781 ],
782 crs=ccrs.PlateCarree(),
783 )
785 # Ensure square aspect
786 axins.set_aspect("equal")
788 else:
789 # Add title.
790 axes.set_title(title, fontsize=16)
792 # Adjust padding if spatial plot or transect
793 if is_transect(cube):
794 yinfopad = -0.1
795 ycbarpad = 0.1
796 else:
797 yinfopad = 0.01
798 ycbarpad = 0.042
800 # Add watermark with min/max/mean. Currently not user togglable.
801 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
802 axes.annotate(
803 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
804 xy=(1, yinfopad),
805 xycoords="axes fraction",
806 xytext=(-5, 5),
807 textcoords="offset points",
808 ha="right",
809 va="bottom",
810 size=11,
811 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
812 )
814 # Add secondary colour bar for overlay_cube field if required.
815 if overlay_cube:
816 cbarB = fig.colorbar(
817 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
818 )
819 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
820 # add ticks and tick_labels for every levels if less than 20 levels exist
821 if over_levels is not None and len(over_levels) < 20: 821 ↛ 822line 821 didn't jump to line 822 because the condition on line 821 was never true
822 cbarB.set_ticks(over_levels)
823 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
824 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
825 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
826 logging.debug("Set secondary colorbar ticks and labels.")
828 # Add main colour bar.
829 cbar = fig.colorbar(
830 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
831 )
833 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
834 # add ticks and tick_labels for every levels if less than 20 levels exist
835 if levels is not None and len(levels) < 20:
836 cbar.set_ticks(levels)
837 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
838 if "rainfall" or "snowfall" or "visibility" in cube.name(): 838 ↛ 840line 838 didn't jump to line 840 because the condition on line 838 was always true
839 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
840 logging.debug("Set colorbar ticks and labels.")
842 # Save plot.
843 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
844 logging.info("Saved spatial plot to %s", filename)
845 plt.close(fig)
848def _plot_and_save_postage_stamp_spatial_plot(
849 cube: iris.cube.Cube,
850 filename: str,
851 stamp_coordinate: str,
852 title: str,
853 method: Literal["contourf", "pcolormesh"],
854 overlay_cube: iris.cube.Cube | None = None,
855 contour_cube: iris.cube.Cube | None = None,
856 **kwargs,
857):
858 """Plot postage stamp spatial plots from an ensemble.
860 Parameters
861 ----------
862 cube: Cube
863 Iris cube of data to be plotted. It must have the stamp coordinate.
864 filename: str
865 Filename of the plot to write.
866 stamp_coordinate: str
867 Coordinate that becomes different plots.
868 method: "contourf" | "pcolormesh"
869 The plotting method to use.
870 overlay_cube: Cube, optional
871 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
872 contour_cube: Cube, optional
873 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
875 Raises
876 ------
877 ValueError
878 If the cube doesn't have the right dimensions.
879 """
880 # Use the smallest square grid that will fit the members.
881 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
883 fig = plt.figure(figsize=(10, 10))
885 # Specify the color bar
886 cmap, levels, norm = _colorbar_map_levels(cube)
887 # If overplotting, set required colorbars
888 if overlay_cube: 888 ↛ 889line 888 didn't jump to line 889 because the condition on line 888 was never true
889 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
890 if contour_cube: 890 ↛ 891line 890 didn't jump to line 891 because the condition on line 890 was never true
891 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
893 # Make a subplot for each member.
894 for member, subplot in zip(
895 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
896 ):
897 # Setup subplot map projection, extent and coastlines and borderlines.
898 axes = _setup_spatial_map(
899 member, fig, cmap, grid_size=grid_size, subplot=subplot
900 )
901 if method == "contourf":
902 # Filled contour plot of the field.
903 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
904 elif method == "pcolormesh":
905 if levels is not None:
906 vmin = min(levels)
907 vmax = max(levels)
908 else:
909 raise TypeError("Unknown vmin and vmax range.")
910 vmin, vmax = None, None
911 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
912 # if levels are defined.
913 if norm is not None: 913 ↛ 914line 913 didn't jump to line 914 because the condition on line 913 was never true
914 vmin = None
915 vmax = None
916 # pcolormesh plot of the field.
917 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
918 else:
919 raise ValueError(f"Unknown plotting method: {method}")
921 # Overplot overlay field, if required
922 if overlay_cube: 922 ↛ 923line 922 didn't jump to line 923 because the condition on line 922 was never true
923 try:
924 over_vmin = min(over_levels)
925 over_vmax = max(over_levels)
926 except TypeError:
927 over_vmin, over_vmax = None, None
928 if over_norm is not None:
929 over_vmin = None
930 over_vmax = None
931 iplt.pcolormesh(
932 overlay_cube[member.coord(stamp_coordinate).points[0]],
933 cmap=over_cmap,
934 norm=over_norm,
935 alpha=0.6,
936 vmin=over_vmin,
937 vmax=over_vmax,
938 )
939 # Overplot contour field, if required
940 if contour_cube: 940 ↛ 941line 940 didn't jump to line 941 because the condition on line 940 was never true
941 iplt.contour(
942 contour_cube[member.coord(stamp_coordinate).points[0]],
943 colors="darkgray",
944 levels=cntr_levels,
945 norm=cntr_norm,
946 alpha=0.6,
947 linestyles="--",
948 linewidths=1,
949 )
950 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
951 axes.set_axis_off()
953 # Put the shared colorbar in its own axes.
954 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])
955 colorbar = fig.colorbar(
956 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
957 )
958 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
960 # Overall figure title.
961 fig.suptitle(title, fontsize=16)
963 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
964 logging.info("Saved contour postage stamp plot to %s", filename)
965 plt.close(fig)
968def _plot_and_save_line_series(
969 cubes: iris.cube.CubeList,
970 coords: list[iris.coords.Coord],
971 ensemble_coord: str,
972 filename: str,
973 title: str,
974 **kwargs,
975):
976 """Plot and save a 1D line series.
978 Parameters
979 ----------
980 cubes: Cube or CubeList
981 Cube or CubeList containing the cubes to plot on the y-axis.
982 coords: list[Coord]
983 Coordinates to plot on the x-axis, one per cube.
984 ensemble_coord: str
985 Ensemble coordinate in the cube.
986 filename: str
987 Filename of the plot to write.
988 title: str
989 Plot title.
990 """
991 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
993 model_colors_map = _get_model_colors_map(cubes)
995 # Store min/max ranges.
996 y_levels = []
998 # Check match-up across sequence coords gives consistent sizes
999 _validate_cubes_coords(cubes, coords)
1001 for cube, coord in zip(cubes, coords, strict=True):
1002 label = None
1003 color = "black"
1004 if model_colors_map:
1005 label = cube.attributes.get("model_name")
1006 color = model_colors_map.get(label)
1007 for cube_slice in cube.slices_over(ensemble_coord):
1008 # Label with (control) if part of an ensemble or not otherwise.
1009 if cube_slice.coord(ensemble_coord).points == [0]:
1010 iplt.plot(
1011 coord,
1012 cube_slice,
1013 color=color,
1014 marker="o",
1015 ls="-",
1016 lw=3,
1017 label=f"{label} (control)"
1018 if len(cube.coord(ensemble_coord).points) > 1
1019 else label,
1020 )
1021 # Label with (perturbed) if part of an ensemble and not the control.
1022 else:
1023 iplt.plot(
1024 coord,
1025 cube_slice,
1026 color=color,
1027 ls="-",
1028 lw=1.5,
1029 alpha=0.75,
1030 label=f"{label} (member)",
1031 )
1033 # Calculate the global min/max if multiple cubes are given.
1034 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1035 if levels is not None: 1035 ↛ 1036line 1035 didn't jump to line 1036 because the condition on line 1035 was never true
1036 y_levels.append(min(levels))
1037 y_levels.append(max(levels))
1039 # Get the current axes.
1040 ax = plt.gca()
1042 # Add some labels and tweak the style.
1043 # check if cubes[0] works for single cube if not CubeList
1044 if coords[0].name() == "time":
1045 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1046 else:
1047 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1048 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1049 ax.set_title(title, fontsize=16)
1051 ax.ticklabel_format(axis="y", useOffset=False)
1052 ax.tick_params(axis="x", labelrotation=15)
1053 ax.tick_params(axis="both", labelsize=12)
1055 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1056 if y_levels: 1056 ↛ 1057line 1056 didn't jump to line 1057 because the condition on line 1056 was never true
1057 ax.set_ylim(min(y_levels), max(y_levels))
1058 # Add zero line.
1059 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1060 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1061 logging.debug(
1062 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1063 )
1064 else:
1065 ax.autoscale()
1067 # Add gridlines
1068 ax.grid(linestyle="--", color="grey", linewidth=1)
1069 # Ientify unique labels for legend
1070 handles = list(
1071 {
1072 label: handle
1073 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1074 }.values()
1075 )
1076 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1078 # Save plot.
1079 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1080 logging.info("Saved line plot to %s", filename)
1081 plt.close(fig)
1084def _plot_and_save_line_power_spectrum_series(
1085 cubes: iris.cube.CubeList,
1086 coords: list[iris.coords.Coord],
1087 ensemble_coord: str,
1088 filename: str,
1089 title: str,
1090 series_coordinate: str = None,
1091 **kwargs,
1092):
1093 """Plot and save a 1D line series.
1095 Parameters
1096 ----------
1097 cubes: Cube or CubeList
1098 Cube or CubeList containing the cubes to plot on the y-axis.
1099 coords: list[Coord]
1100 Coordinates to plot on the x-axis, one per cube.
1101 ensemble_coord: str
1102 Ensemble coordinate in the cube.
1103 filename: str
1104 Filename of the plot to write.
1105 title: str
1106 Plot title.
1107 series_coordinate: str, optional
1108 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1109 """
1110 # xn = coords[0].name() # x-axis (e.g. frequency)
1112 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1113 model_colors_map = _get_model_colors_map(cubes)
1114 ax = plt.gca()
1116 # Store min/max ranges.
1117 y_levels = []
1119 line_marker = None
1120 line_width = 1
1122 # for cube, coord in zip(cubes, coords, strict=True):
1123 for cube in iter_maybe(cubes):
1124 # next 2 lines replace chunk of code.
1125 xcoord = select_series_coord(cube, series_coordinate)
1126 xname = xcoord.points
1128 yfield = cube.data # power spectrum
1129 label = None
1130 color = "black"
1131 if model_colors_map: 1131 ↛ 1134line 1131 didn't jump to line 1134 because the condition on line 1131 was always true
1132 label = cube.attributes.get("model_name")
1133 color = model_colors_map.get(label)
1134 for cube_slice in cube.slices_over(ensemble_coord):
1135 # Label with (control) if part of an ensemble or not otherwise.
1136 if cube_slice.coord(ensemble_coord).points == [0]: 1136 ↛ 1150line 1136 didn't jump to line 1150 because the condition on line 1136 was always true
1137 ax.plot(
1138 xname,
1139 yfield,
1140 color=color,
1141 marker=line_marker,
1142 ls="-",
1143 lw=line_width,
1144 label=f"{label} (control)"
1145 if len(cube.coord(ensemble_coord).points) > 1
1146 else label,
1147 )
1148 # Label with (perturbed) if part of an ensemble and not the control.
1149 else:
1150 ax.plot(
1151 xname,
1152 yfield,
1153 color=color,
1154 ls="-",
1155 lw=1.5,
1156 alpha=0.75,
1157 label=f"{label} (member)",
1158 )
1160 # Calculate the global min/max if multiple cubes are given.
1161 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1162 if levels is not None: 1162 ↛ 1163line 1162 didn't jump to line 1163 because the condition on line 1162 was never true
1163 y_levels.append(min(levels))
1164 y_levels.append(max(levels))
1166 # Add some labels and tweak the style.
1167 # check if cubes[0] works for single cube if not CubeList
1169 # title = f"{title}\n [{coords.units.title(coords.points[0])}]"
1170 title = f"{title}"
1171 ax.set_title(title, fontsize=16)
1173 # Set appropriate x-axis label based on coordinate
1174 if series_coordinate == "wavelength" or ( 1174 ↛ 1177line 1174 didn't jump to line 1177 because the condition on line 1174 was never true
1175 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1176 ):
1177 ax.set_xlabel("Wavelength (km)", fontsize=14)
1178 elif series_coordinate == "physical_wavenumber" or ( 1178 ↛ 1181line 1178 didn't jump to line 1181 because the condition on line 1178 was never true
1179 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1180 ):
1181 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1182 else: # frequency or check units
1183 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1183 ↛ 1184line 1183 didn't jump to line 1184 because the condition on line 1183 was never true
1184 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1185 else:
1186 ax.set_xlabel("Wavenumber", fontsize=14)
1188 ax.set_ylabel("Power Spectral Density", fontsize=14)
1189 ax.tick_params(axis="both", labelsize=12)
1191 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1193 # Set log-log scale
1194 ax.set_xscale("log")
1195 ax.set_yscale("log")
1197 # Add gridlines
1198 ax.grid(linestyle="--", color="grey", linewidth=1)
1199 # Ientify unique labels for legend
1200 handles = list(
1201 {
1202 label: handle
1203 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1204 }.values()
1205 )
1206 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1208 # Save plot.
1209 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1210 logging.info("Saved line plot to %s", filename)
1211 plt.close(fig)
1214def _plot_and_save_vertical_line_series(
1215 cubes: iris.cube.CubeList,
1216 coords: list[iris.coords.Coord],
1217 ensemble_coord: str,
1218 filename: str,
1219 series_coordinate: str,
1220 title: str,
1221 vmin: float,
1222 vmax: float,
1223 **kwargs,
1224):
1225 """Plot and save a 1D line series in vertical.
1227 Parameters
1228 ----------
1229 cubes: CubeList
1230 1 dimensional Cube or CubeList of the data to plot on x-axis.
1231 coord: list[Coord]
1232 Coordinates to plot on the y-axis, one per cube.
1233 ensemble_coord: str
1234 Ensemble coordinate in the cube.
1235 filename: str
1236 Filename of the plot to write.
1237 series_coordinate: str
1238 Coordinate to use as vertical axis.
1239 title: str
1240 Plot title.
1241 vmin: float
1242 Minimum value for the x-axis.
1243 vmax: float
1244 Maximum value for the x-axis.
1245 """
1246 # plot the vertical pressure axis using log scale
1247 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1249 model_colors_map = _get_model_colors_map(cubes)
1251 # Check match-up across sequence coords gives consistent sizes
1252 _validate_cubes_coords(cubes, coords)
1254 for cube, coord in zip(cubes, coords, strict=True):
1255 label = None
1256 color = "black"
1257 if model_colors_map: 1257 ↛ 1258line 1257 didn't jump to line 1258 because the condition on line 1257 was never true
1258 label = cube.attributes.get("model_name")
1259 color = model_colors_map.get(label)
1261 for cube_slice in cube.slices_over(ensemble_coord):
1262 # If ensemble data given plot control member with (control)
1263 # unless single forecast.
1264 if cube_slice.coord(ensemble_coord).points == [0]:
1265 iplt.plot(
1266 cube_slice,
1267 coord,
1268 color=color,
1269 marker="o",
1270 ls="-",
1271 lw=3,
1272 label=f"{label} (control)"
1273 if len(cube.coord(ensemble_coord).points) > 1
1274 else label,
1275 )
1276 # If ensemble data given plot perturbed members with (perturbed).
1277 else:
1278 iplt.plot(
1279 cube_slice,
1280 coord,
1281 color=color,
1282 ls="-",
1283 lw=1.5,
1284 alpha=0.75,
1285 label=f"{label} (member)",
1286 )
1288 # Get the current axis
1289 ax = plt.gca()
1291 # Special handling for pressure level data.
1292 if series_coordinate == "pressure": 1292 ↛ 1314line 1292 didn't jump to line 1314 because the condition on line 1292 was always true
1293 # Invert y-axis and set to log scale.
1294 ax.invert_yaxis()
1295 ax.set_yscale("log")
1297 # Define y-ticks and labels for pressure log axis.
1298 y_tick_labels = [
1299 "1000",
1300 "850",
1301 "700",
1302 "500",
1303 "300",
1304 "200",
1305 "100",
1306 ]
1307 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1309 # Set y-axis limits and ticks.
1310 ax.set_ylim(1100, 100)
1312 # Test if series_coordinate is model level data. The UM data uses
1313 # model_level_number and lfric uses full_levels as coordinate.
1314 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1315 # Define y-ticks and labels for vertical axis.
1316 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1317 y_tick_labels = [str(int(i)) for i in y_ticks]
1318 ax.set_ylim(min(y_ticks), max(y_ticks))
1320 ax.set_yticks(y_ticks)
1321 ax.set_yticklabels(y_tick_labels)
1323 # Set x-axis limits.
1324 ax.set_xlim(vmin, vmax)
1325 # Mark y=0 if present in plot.
1326 if vmin < 0.0 and vmax > 0.0: 1326 ↛ 1327line 1326 didn't jump to line 1327 because the condition on line 1326 was never true
1327 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1329 # Add some labels and tweak the style.
1330 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1331 ax.set_xlabel(
1332 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1333 )
1334 ax.set_title(title, fontsize=16)
1335 ax.ticklabel_format(axis="x")
1336 ax.tick_params(axis="y")
1337 ax.tick_params(axis="both", labelsize=12)
1339 # Add gridlines
1340 ax.grid(linestyle="--", color="grey", linewidth=1)
1341 # Ientify unique labels for legend
1342 handles = list(
1343 {
1344 label: handle
1345 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1346 }.values()
1347 )
1348 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1350 # Save plot.
1351 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1352 logging.info("Saved line plot to %s", filename)
1353 plt.close(fig)
1356def _plot_and_save_scatter_plot(
1357 cube_x: iris.cube.Cube | iris.cube.CubeList,
1358 cube_y: iris.cube.Cube | iris.cube.CubeList,
1359 filename: str,
1360 title: str,
1361 one_to_one: bool,
1362 model_names: list[str] = None,
1363 **kwargs,
1364):
1365 """Plot and save a 2D scatter plot.
1367 Parameters
1368 ----------
1369 cube_x: Cube | CubeList
1370 1 dimensional Cube or CubeList of the data to plot on x-axis.
1371 cube_y: Cube | CubeList
1372 1 dimensional Cube or CubeList of the data to plot on y-axis.
1373 filename: str
1374 Filename of the plot to write.
1375 title: str
1376 Plot title.
1377 one_to_one: bool
1378 Whether a 1:1 line is plotted.
1379 """
1380 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1381 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1382 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1383 # over the pairs simultaneously.
1385 # Ensure cube_x and cube_y are iterable
1386 cube_x_iterable = iter_maybe(cube_x)
1387 cube_y_iterable = iter_maybe(cube_y)
1389 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1390 iplt.scatter(cube_x_iter, cube_y_iter)
1391 if one_to_one is True:
1392 plt.plot(
1393 [
1394 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1395 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1396 ],
1397 [
1398 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1399 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1400 ],
1401 "k",
1402 linestyle="--",
1403 )
1404 ax = plt.gca()
1406 # Add some labels and tweak the style.
1407 if model_names is None:
1408 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1409 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1410 else:
1411 # Add the model names, these should be order of base (x) and other (y).
1412 ax.set_xlabel(
1413 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1414 )
1415 ax.set_ylabel(
1416 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1417 )
1418 ax.set_title(title, fontsize=16)
1419 ax.ticklabel_format(axis="y", useOffset=False)
1420 ax.tick_params(axis="x", labelrotation=15)
1421 ax.tick_params(axis="both", labelsize=12)
1422 ax.autoscale()
1424 # Save plot.
1425 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1426 logging.info("Saved scatter plot to %s", filename)
1427 plt.close(fig)
1430def _plot_and_save_vector_plot(
1431 cube_u: iris.cube.Cube,
1432 cube_v: iris.cube.Cube,
1433 filename: str,
1434 title: str,
1435 method: Literal["contourf", "pcolormesh"],
1436 **kwargs,
1437):
1438 """Plot and save a 2D vector plot.
1440 Parameters
1441 ----------
1442 cube_u: Cube
1443 2 dimensional Cube of u component of the data.
1444 cube_v: Cube
1445 2 dimensional Cube of v component of the data.
1446 filename: str
1447 Filename of the plot to write.
1448 title: str
1449 Plot title.
1450 """
1451 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1453 # Create a cube containing the magnitude of the vector field.
1454 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1455 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1457 # Specify the color bar
1458 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1460 # Setup plot map projection, extent and coastlines and borderlines.
1461 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1463 if method == "contourf":
1464 # Filled contour plot of the field.
1465 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1466 elif method == "pcolormesh":
1467 try:
1468 vmin = min(levels)
1469 vmax = max(levels)
1470 except TypeError:
1471 vmin, vmax = None, None
1472 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1473 # if levels are defined.
1474 if norm is not None:
1475 vmin = None
1476 vmax = None
1477 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1478 else:
1479 raise ValueError(f"Unknown plotting method: {method}")
1481 # Check to see if transect, and if so, adjust y axis.
1482 if is_transect(cube_vec_mag):
1483 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1484 axes.invert_yaxis()
1485 axes.set_yscale("log")
1486 axes.set_ylim(1100, 100)
1487 # If both model_level_number and level_height exists, iplt can construct
1488 # plot as a function of height above orography (NOT sea level).
1489 elif {"model_level_number", "level_height"}.issubset(
1490 {coord.name() for coord in cube_vec_mag.coords()}
1491 ):
1492 axes.set_yscale("log")
1494 axes.set_title(
1495 f"{title}\n"
1496 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1497 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1498 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1499 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1500 fontsize=16,
1501 )
1503 else:
1504 # Add title.
1505 axes.set_title(title, fontsize=16)
1507 # Add watermark with min/max/mean. Currently not user togglable.
1508 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1509 axes.annotate(
1510 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}",
1511 xy=(1, -0.05),
1512 xycoords="axes fraction",
1513 xytext=(-5, 5),
1514 textcoords="offset points",
1515 ha="right",
1516 va="bottom",
1517 size=11,
1518 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1519 )
1521 # Add colour bar.
1522 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1523 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1524 # add ticks and tick_labels for every levels if less than 20 levels exist
1525 if levels is not None and len(levels) < 20:
1526 cbar.set_ticks(levels)
1527 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1529 # 30 barbs along the longest axis of the plot, or a barb per point for data
1530 # with less than 30 points.
1531 step = max(max(cube_u.shape) // 30, 1)
1532 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1534 # Save plot.
1535 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1536 logging.info("Saved vector plot to %s", filename)
1537 plt.close(fig)
1540def _plot_and_save_histogram_series(
1541 cubes: iris.cube.Cube | iris.cube.CubeList,
1542 filename: str,
1543 title: str,
1544 vmin: float,
1545 vmax: float,
1546 **kwargs,
1547):
1548 """Plot and save a histogram series.
1550 Parameters
1551 ----------
1552 cubes: Cube or CubeList
1553 2 dimensional Cube or CubeList of the data to plot as histogram.
1554 filename: str
1555 Filename of the plot to write.
1556 title: str
1557 Plot title.
1558 vmin: float
1559 minimum for colorbar
1560 vmax: float
1561 maximum for colorbar
1562 """
1563 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1564 ax = plt.gca()
1566 model_colors_map = _get_model_colors_map(cubes)
1568 # Set default that histograms will produce probability density function
1569 # at each bin (integral over range sums to 1).
1570 density = True
1572 for cube in iter_maybe(cubes):
1573 # Easier to check title (where var name originates)
1574 # than seeing if long names exist etc.
1575 # Exception case, where distribution better fits log scales/bins.
1576 if "surface_microphysical" in title:
1577 if "amount" in title: 1577 ↛ 1579line 1577 didn't jump to line 1579 because the condition on line 1577 was never true
1578 # Compute histogram following Klingaman et al. (2017): ASoP
1579 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1580 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1581 density = False
1582 else:
1583 bins = 10.0 ** (
1584 np.arange(-10, 27, 1) / 10.0
1585 ) # Suggestion from RMED toolbox.
1586 bins = np.insert(bins, 0, 0)
1587 ax.set_yscale("log")
1588 vmin = bins[1]
1589 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1590 ax.set_xscale("log")
1591 elif "lightning" in title:
1592 bins = [0, 1, 2, 3, 4, 5]
1593 else:
1594 bins = np.linspace(vmin, vmax, 51)
1595 logging.debug(
1596 "Plotting histogram with %s bins %s - %s.",
1597 np.size(bins),
1598 np.min(bins),
1599 np.max(bins),
1600 )
1602 # Reshape cube data into a single array to allow for a single histogram.
1603 # Otherwise we plot xdim histograms stacked.
1604 cube_data_1d = (cube.data).flatten()
1606 label = None
1607 color = "black"
1608 if model_colors_map: 1608 ↛ 1609line 1608 didn't jump to line 1609 because the condition on line 1608 was never true
1609 label = cube.attributes.get("model_name")
1610 color = model_colors_map[label]
1611 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1613 # Compute area under curve.
1614 if "surface_microphysical" in title and "amount" in title: 1614 ↛ 1615line 1614 didn't jump to line 1615 because the condition on line 1614 was never true
1615 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1616 x = x * bin_mean / x.sum()
1617 x = x[1:]
1618 y = y[1:]
1620 ax.plot(
1621 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1622 )
1624 # Add some labels and tweak the style.
1625 ax.set_title(title, fontsize=16)
1626 ax.set_xlabel(
1627 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1628 )
1629 ax.set_ylabel("Normalised probability density", fontsize=14)
1630 if "surface_microphysical" in title and "amount" in title: 1630 ↛ 1631line 1630 didn't jump to line 1631 because the condition on line 1630 was never true
1631 ax.set_ylabel(
1632 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1633 )
1634 ax.set_xlim(vmin, vmax)
1635 ax.tick_params(axis="both", labelsize=12)
1637 # Overlay grid-lines onto histogram plot.
1638 ax.grid(linestyle="--", color="grey", linewidth=1)
1639 if model_colors_map: 1639 ↛ 1640line 1639 didn't jump to line 1640 because the condition on line 1639 was never true
1640 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1642 # Save plot.
1643 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1644 logging.info("Saved histogram plot to %s", filename)
1645 plt.close(fig)
1648def _plot_and_save_postage_stamp_histogram_series(
1649 cube: iris.cube.Cube,
1650 filename: str,
1651 title: str,
1652 stamp_coordinate: str,
1653 vmin: float,
1654 vmax: float,
1655 **kwargs,
1656):
1657 """Plot and save postage (ensemble members) stamps for a histogram series.
1659 Parameters
1660 ----------
1661 cube: Cube
1662 2 dimensional Cube of the data to plot as histogram.
1663 filename: str
1664 Filename of the plot to write.
1665 title: str
1666 Plot title.
1667 stamp_coordinate: str
1668 Coordinate that becomes different plots.
1669 vmin: float
1670 minimum for pdf x-axis
1671 vmax: float
1672 maximum for pdf x-axis
1673 """
1674 # Use the smallest square grid that will fit the members.
1675 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1677 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1678 # Make a subplot for each member.
1679 for member, subplot in zip(
1680 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1681 ):
1682 # Implicit interface is much easier here, due to needing to have the
1683 # cartopy GeoAxes generated.
1684 plt.subplot(grid_size, grid_size, subplot)
1685 # Reshape cube data into a single array to allow for a single histogram.
1686 # Otherwise we plot xdim histograms stacked.
1687 member_data_1d = (member.data).flatten()
1688 plt.hist(member_data_1d, density=True, stacked=True)
1689 ax = plt.gca()
1690 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1691 ax.set_xlim(vmin, vmax)
1693 # Overall figure title.
1694 fig.suptitle(title, fontsize=16)
1696 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1697 logging.info("Saved histogram postage stamp plot to %s", filename)
1698 plt.close(fig)
1701def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1702 cube: iris.cube.Cube,
1703 filename: str,
1704 title: str,
1705 stamp_coordinate: str,
1706 vmin: float,
1707 vmax: float,
1708 **kwargs,
1709):
1710 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1711 ax.set_title(title, fontsize=16)
1712 ax.set_xlim(vmin, vmax)
1713 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1714 ax.set_ylabel("normalised probability density", fontsize=14)
1715 # Loop over all slices along the stamp_coordinate
1716 for member in cube.slices_over(stamp_coordinate):
1717 # Flatten the member data to 1D
1718 member_data_1d = member.data.flatten()
1719 # Plot the histogram using plt.hist
1720 plt.hist(
1721 member_data_1d,
1722 density=True,
1723 stacked=True,
1724 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1725 )
1727 # Add a legend
1728 ax.legend(fontsize=16)
1730 # Save the figure to a file
1731 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1732 logging.info("Saved histogram postage stamp plot to %s", filename)
1734 # Close the figure
1735 plt.close(fig)
1738def _plot_and_save_scattermap_plot(
1739 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1740):
1741 """Plot and save a geographical scatter plot.
1743 Parameters
1744 ----------
1745 cube: Cube
1746 1 dimensional Cube of the data points with auxiliary latitude and
1747 longitude coordinates,
1748 filename: str
1749 Filename of the plot to write.
1750 title: str
1751 Plot title.
1752 projection: str
1753 Mapping projection to be used by cartopy.
1754 """
1755 # Setup plot details, size, resolution, etc.
1756 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1757 if projection is not None:
1758 # Apart from the default, the only projection we currently support is
1759 # a stereographic projection over the North Pole.
1760 if projection == "NP_Stereo":
1761 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1762 else:
1763 raise ValueError(f"Unknown projection: {projection}")
1764 else:
1765 axes = plt.axes(projection=ccrs.PlateCarree())
1767 # Scatter plot of the field. The marker size is chosen to give
1768 # symbols that decrease in size as the number of observations
1769 # increases, although the fraction of the figure covered by
1770 # symbols increases roughly as N^(1/2), disregarding overlaps,
1771 # and has been selected for the default figure size of (10, 10).
1772 # Should this be changed, the marker size should be adjusted in
1773 # proportion to the area of the figure.
1774 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1775 klon = None
1776 klat = None
1777 for kc in range(len(cube.aux_coords)):
1778 if cube.aux_coords[kc].standard_name == "latitude":
1779 klat = kc
1780 elif cube.aux_coords[kc].standard_name == "longitude":
1781 klon = kc
1782 scatter_map = iplt.scatter(
1783 cube.aux_coords[klon],
1784 cube.aux_coords[klat],
1785 c=cube.data[:],
1786 s=mrk_size,
1787 cmap="jet",
1788 edgecolors="k",
1789 )
1791 # Add coastlines and borderlines.
1792 try:
1793 axes.coastlines(resolution="10m")
1794 axes.add_feature(cfeature.BORDERS)
1795 except AttributeError:
1796 pass
1798 # Add title.
1799 axes.set_title(title, fontsize=16)
1801 # Add colour bar.
1802 cbar = fig.colorbar(scatter_map)
1803 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1805 # Save plot.
1806 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1807 logging.info("Saved geographical scatter plot to %s", filename)
1808 plt.close(fig)
1811def _spatial_plot(
1812 method: Literal["contourf", "pcolormesh"],
1813 cube: iris.cube.Cube,
1814 filename: str | None,
1815 sequence_coordinate: str,
1816 stamp_coordinate: str,
1817 overlay_cube: iris.cube.Cube | None = None,
1818 contour_cube: iris.cube.Cube | None = None,
1819 **kwargs,
1820):
1821 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1823 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1824 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1825 is present then postage stamp plots will be produced.
1827 If an overlay_cube and/or contour_cube are specified, multiple variables can
1828 be overplotted on the same figure.
1830 Parameters
1831 ----------
1832 method: "contourf" | "pcolormesh"
1833 The plotting method to use.
1834 cube: Cube
1835 Iris cube of the data to plot. It should have two spatial dimensions,
1836 such as lat and lon, and may also have a another two dimension to be
1837 plotted sequentially and/or as postage stamp plots.
1838 filename: str | None
1839 Name of the plot to write, used as a prefix for plot sequences. If None
1840 uses the recipe name.
1841 sequence_coordinate: str
1842 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1843 This coordinate must exist in the cube.
1844 stamp_coordinate: str
1845 Coordinate about which to plot postage stamp plots. Defaults to
1846 ``"realization"``.
1847 overlay_cube: Cube | None, optional
1848 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1849 contour_cube: Cube | None, optional
1850 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1852 Raises
1853 ------
1854 ValueError
1855 If the cube doesn't have the right dimensions.
1856 TypeError
1857 If the cube isn't a single cube.
1858 """
1859 recipe_title = get_recipe_metadata().get("title", "Untitled")
1861 # Ensure we've got a single cube.
1862 cube = _check_single_cube(cube)
1864 # Make postage stamp plots if stamp_coordinate exists and has more than a
1865 # single point.
1866 plotting_func = _plot_and_save_spatial_plot
1867 try:
1868 if cube.coord(stamp_coordinate).shape[0] > 1:
1869 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1870 except iris.exceptions.CoordinateNotFoundError:
1871 pass
1873 # Produce a geographical scatter plot if the data have a
1874 # dimension called observation or model_obs_error
1875 if any( 1875 ↛ 1879line 1875 didn't jump to line 1879 because the condition on line 1875 was never true
1876 crd.var_name == "station" or crd.var_name == "model_obs_error"
1877 for crd in cube.coords()
1878 ):
1879 plotting_func = _plot_and_save_scattermap_plot
1881 # Must have a sequence coordinate.
1882 try:
1883 cube.coord(sequence_coordinate)
1884 except iris.exceptions.CoordinateNotFoundError as err:
1885 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1887 # Create a plot for each value of the sequence coordinate.
1888 plot_index = []
1889 nplot = np.size(cube.coord(sequence_coordinate).points)
1891 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1892 # Set plot titles and filename
1893 seq_coord = cube_slice.coord(sequence_coordinate)
1894 plot_title, plot_filename = _set_title_and_filename(
1895 seq_coord, nplot, recipe_title, filename
1896 )
1898 # Extract sequence slice for overlay_cube and contour_cube if required.
1899 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1900 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1902 # Do the actual plotting.
1903 plotting_func(
1904 cube_slice,
1905 filename=plot_filename,
1906 stamp_coordinate=stamp_coordinate,
1907 title=plot_title,
1908 method=method,
1909 overlay_cube=overlay_slice,
1910 contour_cube=contour_slice,
1911 **kwargs,
1912 )
1913 plot_index.append(plot_filename)
1915 # Add list of plots to plot metadata.
1916 complete_plot_index = _append_to_plot_index(plot_index)
1918 # Make a page to display the plots.
1919 _make_plot_html_page(complete_plot_index)
1922def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1923 """Get colourmap for mask.
1925 If "mask_for_" appears anywhere in the name of a cube this function will be called
1926 regardless of the name of the variable to ensure a consistent plot.
1928 Parameters
1929 ----------
1930 cube: Cube
1931 Cube of variable for which the colorbar information is desired.
1932 axis: "x", "y", optional
1933 Select the levels for just this axis of a line plot. The min and max
1934 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1935 setting a universal range is not desirable (e.g. temperature), users
1936 can set ymin/ymax values to "auto" in the colorbar definitions file.
1937 Where no additional xmin/xmax or ymin/ymax values are provided, the
1938 axis bounds default to use the vmin/vmax values provided.
1940 Returns
1941 -------
1942 cmap:
1943 Matplotlib colormap.
1944 levels:
1945 List of levels to use for plotting. For continuous plots the min and max
1946 should be taken as the range.
1947 norm:
1948 BoundaryNorm information.
1949 """
1950 if "difference" not in cube.long_name:
1951 if axis:
1952 levels = [0, 1]
1953 # Complete settings based on levels.
1954 return None, levels, None
1955 else:
1956 # Define the levels and colors.
1957 levels = [0, 1, 2]
1958 colors = ["white", "dodgerblue"]
1959 # Create a custom color map.
1960 cmap = mcolors.ListedColormap(colors)
1961 # Normalize the levels.
1962 norm = mcolors.BoundaryNorm(levels, cmap.N)
1963 logging.debug("Colourmap for %s.", cube.long_name)
1964 return cmap, levels, norm
1965 else:
1966 if axis:
1967 levels = [-1, 1]
1968 return None, levels, None
1969 else:
1970 # Search for if mask difference, set to +/- 0.5 as values plotted <
1971 # not <=.
1972 levels = [-2, -0.5, 0.5, 2]
1973 colors = ["goldenrod", "white", "teal"]
1974 cmap = mcolors.ListedColormap(colors)
1975 norm = mcolors.BoundaryNorm(levels, cmap.N)
1976 logging.debug("Colourmap for %s.", cube.long_name)
1977 return cmap, levels, norm
1980def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1981 """Get a custom colorbar for a cube in the Beaufort Scale.
1983 Specific variable ranges can be separately set in user-supplied definition
1984 for x- or y-axis limits, or indicate where automated range preferred.
1986 Parameters
1987 ----------
1988 cube: Cube
1989 Cube of variable with Beaufort Scale in name.
1990 axis: "x", "y", optional
1991 Select the levels for just this axis of a line plot. The min and max
1992 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1993 setting a universal range is not desirable (e.g. temperature), users
1994 can set ymin/ymax values to "auto" in the colorbar definitions file.
1995 Where no additional xmin/xmax or ymin/ymax values are provided, the
1996 axis bounds default to use the vmin/vmax values provided.
1998 Returns
1999 -------
2000 cmap:
2001 Matplotlib colormap.
2002 levels:
2003 List of levels to use for plotting. For continuous plots the min and max
2004 should be taken as the range.
2005 norm:
2006 BoundaryNorm information.
2007 """
2008 if "difference" not in cube.long_name:
2009 if axis:
2010 levels = [0, 12]
2011 return None, levels, None
2012 else:
2013 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
2014 colors = [
2015 "black",
2016 (0, 0, 0.6),
2017 "blue",
2018 "cyan",
2019 "green",
2020 "yellow",
2021 (1, 0.5, 0),
2022 "red",
2023 "pink",
2024 "magenta",
2025 "purple",
2026 "maroon",
2027 "white",
2028 ]
2029 cmap = mcolors.ListedColormap(colors)
2030 norm = mcolors.BoundaryNorm(levels, cmap.N)
2031 logging.info("change colormap for Beaufort Scale colorbar.")
2032 return cmap, levels, norm
2033 else:
2034 if axis:
2035 levels = [-4, 4]
2036 return None, levels, None
2037 else:
2038 levels = [
2039 -3.5,
2040 -2.5,
2041 -1.5,
2042 -0.5,
2043 0.5,
2044 1.5,
2045 2.5,
2046 3.5,
2047 ]
2048 cmap = plt.get_cmap("bwr", 8)
2049 norm = mcolors.BoundaryNorm(levels, cmap.N)
2050 return cmap, levels, norm
2053def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
2054 """Return altered colourmap for temperature with change in units to Celsius.
2056 If "Celsius" appears anywhere in the name of a cube this function will be called.
2058 Parameters
2059 ----------
2060 cube: Cube
2061 Cube of variable for which the colorbar information is desired.
2062 cmap: Matplotlib colormap.
2063 levels: List
2064 List of levels to use for plotting. For continuous plots the min and max
2065 should be taken as the range.
2066 norm: BoundaryNorm.
2068 Returns
2069 -------
2070 cmap: Matplotlib colormap.
2071 levels: List
2072 List of levels to use for plotting. For continuous plots the min and max
2073 should be taken as the range.
2074 norm: BoundaryNorm.
2075 """
2076 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2077 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
2078 levels = np.array(levels)
2079 levels -= 273
2080 levels = levels.tolist()
2081 else:
2082 # Do nothing keep the existing colourbar attributes
2083 levels = levels
2084 cmap = cmap
2085 norm = norm
2086 return cmap, levels, norm
2089def _custom_colormap_probability(
2090 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
2091):
2092 """Get a custom colorbar for a probability cube.
2094 Specific variable ranges can be separately set in user-supplied definition
2095 for x- or y-axis limits, or indicate where automated range preferred.
2097 Parameters
2098 ----------
2099 cube: Cube
2100 Cube of variable with probability in name.
2101 axis: "x", "y", optional
2102 Select the levels for just this axis of a line plot. The min and max
2103 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2104 setting a universal range is not desirable (e.g. temperature), users
2105 can set ymin/ymax values to "auto" in the colorbar definitions file.
2106 Where no additional xmin/xmax or ymin/ymax values are provided, the
2107 axis bounds default to use the vmin/vmax values provided.
2109 Returns
2110 -------
2111 cmap:
2112 Matplotlib colormap.
2113 levels:
2114 List of levels to use for plotting. For continuous plots the min and max
2115 should be taken as the range.
2116 norm:
2117 BoundaryNorm information.
2118 """
2119 if axis:
2120 levels = [0, 1]
2121 return None, levels, None
2122 else:
2123 cmap = mcolors.ListedColormap(
2124 [
2125 "#FFFFFF",
2126 "#636363",
2127 "#e1dada",
2128 "#B5CAFF",
2129 "#8FB3FF",
2130 "#7F97FF",
2131 "#ABCF63",
2132 "#E8F59E",
2133 "#FFFA14",
2134 "#FFD121",
2135 "#FFA30A",
2136 ]
2137 )
2138 levels = [0.0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
2139 norm = mcolors.BoundaryNorm(levels, cmap.N)
2140 return cmap, levels, norm
2143def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
2144 """Return a custom colourmap for the current recipe."""
2145 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2146 if (
2147 any("surface_microphysical" in name for name in varnames)
2148 and "difference" not in cube.long_name
2149 and "mask" not in cube.long_name
2150 ):
2151 # Define the levels and colors
2152 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
2153 colors = [
2154 "w",
2155 (0, 0, 0.6),
2156 "b",
2157 "c",
2158 "g",
2159 "y",
2160 (1, 0.5, 0),
2161 "r",
2162 "pink",
2163 "m",
2164 "purple",
2165 "maroon",
2166 "gray",
2167 ]
2168 # Create a custom colormap
2169 cmap = mcolors.ListedColormap(colors)
2170 # Normalize the levels
2171 norm = mcolors.BoundaryNorm(levels, cmap.N)
2172 logging.info("change colormap for surface_microphysical variable colorbar.")
2173 else:
2174 # do nothing and keep existing colorbar attributes
2175 cmap = cmap
2176 levels = levels
2177 norm = norm
2178 return cmap, levels, norm
2181def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
2182 """Return custom colourmap for aviation colour state.
2184 If "aviation_colour_state" appears anywhere in the name of a cube
2185 this function will be called.
2187 Parameters
2188 ----------
2189 cube: Cube
2190 Cube of variable for which the colorbar information is desired.
2192 Returns
2193 -------
2194 cmap: Matplotlib colormap.
2195 levels: List
2196 List of levels to use for plotting. For continuous plots the min and max
2197 should be taken as the range.
2198 norm: BoundaryNorm.
2199 """
2200 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
2201 colors = [
2202 "#87ceeb",
2203 "#ffffff",
2204 "#8ced69",
2205 "#ffff00",
2206 "#ffd700",
2207 "#ffa500",
2208 "#fe3620",
2209 ]
2210 # Create a custom colormap
2211 cmap = mcolors.ListedColormap(colors)
2212 # Normalise the levels
2213 norm = mcolors.BoundaryNorm(levels, cmap.N)
2214 return cmap, levels, norm
2217def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
2218 """Return a custom colourmap for the current recipe."""
2219 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2220 if (
2221 any("visibility_in_air" in name for name in varnames)
2222 and "difference" not in cube.long_name
2223 and "mask" not in cube.long_name
2224 ):
2225 # Define the levels and colors (in km)
2226 levels = [0, 0.05, 0.1, 0.2, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 50.0, 70.0, 100.0]
2227 norm = mcolors.BoundaryNorm(levels, cmap.N)
2228 colours = [
2229 "#8f00d6",
2230 "#d10000",
2231 "#ff9700",
2232 "#ffff00",
2233 "#00007f",
2234 "#6c9ccd",
2235 "#aae8ff",
2236 "#37a648",
2237 "#8edc64",
2238 "#c5ffc5",
2239 "#dcdcdc",
2240 "#ffffff",
2241 ]
2242 # Create a custom colormap
2243 cmap = mcolors.ListedColormap(colours)
2244 # Normalize the levels
2245 norm = mcolors.BoundaryNorm(levels, cmap.N)
2246 logging.info("change colormap for visibility_in_air variable colorbar.")
2247 else:
2248 # do nothing and keep existing colorbar attributes
2249 cmap = cmap
2250 levels = levels
2251 norm = norm
2252 return cmap, levels, norm
2255def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
2256 """Return number of models based on cube attributes."""
2257 model_names = list(
2258 filter(
2259 lambda x: x is not None,
2260 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
2261 )
2262 )
2263 if not model_names:
2264 logging.debug("Missing model names. Will assume single model.")
2265 return 1
2266 else:
2267 return len(model_names)
2270def _validate_cube_shape(
2271 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2272) -> None:
2273 """Check all cubes have a model name."""
2274 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2274 ↛ 2275line 2274 didn't jump to line 2275 because the condition on line 2274 was never true
2275 raise ValueError(
2276 f"The number of model names ({num_models}) should equal the number "
2277 f"of cubes ({len(cube)})."
2278 )
2281def _validate_cubes_coords(
2282 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2283) -> None:
2284 """Check same number of cubes as sequence coordinate for zip functions."""
2285 if len(cubes) != len(coords): 2285 ↛ 2286line 2285 didn't jump to line 2286 because the condition on line 2285 was never true
2286 raise ValueError(
2287 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2288 f"of sequence coordinates ({len(coords)})."
2289 f"Check that number of time entries in input data are consistent if "
2290 f"performing time-averaging steps prior to plotting outputs."
2291 )
2294####################
2295# Public functions #
2296####################
2299def spatial_contour_plot(
2300 cube: iris.cube.Cube,
2301 filename: str = None,
2302 sequence_coordinate: str = "time",
2303 stamp_coordinate: str = "realization",
2304 **kwargs,
2305) -> iris.cube.Cube:
2306 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2308 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2309 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2310 is present then postage stamp plots will be produced.
2312 Parameters
2313 ----------
2314 cube: Cube
2315 Iris cube of the data to plot. It should have two spatial dimensions,
2316 such as lat and lon, and may also have a another two dimension to be
2317 plotted sequentially and/or as postage stamp plots.
2318 filename: str, optional
2319 Name of the plot to write, used as a prefix for plot sequences. Defaults
2320 to the recipe name.
2321 sequence_coordinate: str, optional
2322 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2323 This coordinate must exist in the cube.
2324 stamp_coordinate: str, optional
2325 Coordinate about which to plot postage stamp plots. Defaults to
2326 ``"realization"``.
2328 Returns
2329 -------
2330 Cube
2331 The original cube (so further operations can be applied).
2333 Raises
2334 ------
2335 ValueError
2336 If the cube doesn't have the right dimensions.
2337 TypeError
2338 If the cube isn't a single cube.
2339 """
2340 _spatial_plot(
2341 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2342 )
2343 return cube
2346def spatial_pcolormesh_plot(
2347 cube: iris.cube.Cube,
2348 filename: str = None,
2349 sequence_coordinate: str = "time",
2350 stamp_coordinate: str = "realization",
2351 **kwargs,
2352) -> iris.cube.Cube:
2353 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2355 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2356 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2357 is present then postage stamp plots will be produced.
2359 This function is significantly faster than ``spatial_contour_plot``,
2360 especially at high resolutions, and should be preferred unless contiguous
2361 contour areas are important.
2363 Parameters
2364 ----------
2365 cube: Cube
2366 Iris cube of the data to plot. It should have two spatial dimensions,
2367 such as lat and lon, and may also have a another two dimension to be
2368 plotted sequentially and/or as postage stamp plots.
2369 filename: str, optional
2370 Name of the plot to write, used as a prefix for plot sequences. Defaults
2371 to the recipe name.
2372 sequence_coordinate: str, optional
2373 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2374 This coordinate must exist in the cube.
2375 stamp_coordinate: str, optional
2376 Coordinate about which to plot postage stamp plots. Defaults to
2377 ``"realization"``.
2379 Returns
2380 -------
2381 Cube
2382 The original cube (so further operations can be applied).
2384 Raises
2385 ------
2386 ValueError
2387 If the cube doesn't have the right dimensions.
2388 TypeError
2389 If the cube isn't a single cube.
2390 """
2391 _spatial_plot(
2392 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2393 )
2394 return cube
2397def spatial_multi_pcolormesh_plot(
2398 cube: iris.cube.Cube,
2399 overlay_cube: iris.cube.Cube,
2400 contour_cube: iris.cube.Cube,
2401 filename: str = None,
2402 sequence_coordinate: str = "time",
2403 stamp_coordinate: str = "realization",
2404 **kwargs,
2405) -> iris.cube.Cube:
2406 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2408 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2409 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2410 is present then postage stamp plots will be produced.
2412 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2414 If specified, contours of a contour_cube can be overplotted on top of those.
2416 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2418 This function is significantly faster than ``spatial_contour_plot``,
2419 especially at high resolutions, and should be preferred unless contiguous
2420 contour areas are important.
2422 Parameters
2423 ----------
2424 cube: Cube
2425 Iris cube of the data to plot. It should have two spatial dimensions,
2426 such as lat and lon, and may also have a another two dimension to be
2427 plotted sequentially and/or as postage stamp plots.
2428 overlay_cube: Cube
2429 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2430 such as lat and lon, and may also have a another two dimension to be
2431 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.
2432 contour_cube: Cube
2433 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,
2434 such as lat and lon, and may also have a another two dimension to be
2435 plotted sequentially and/or as postage stamp plots.
2436 filename: str, optional
2437 Name of the plot to write, used as a prefix for plot sequences. Defaults
2438 to the recipe name.
2439 sequence_coordinate: str, optional
2440 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2441 This coordinate must exist in the cube.
2442 stamp_coordinate: str, optional
2443 Coordinate about which to plot postage stamp plots. Defaults to
2444 ``"realization"``.
2446 Returns
2447 -------
2448 Cube
2449 The original cube (so further operations can be applied).
2451 Raises
2452 ------
2453 ValueError
2454 If the cube doesn't have the right dimensions.
2455 TypeError
2456 If the cube isn't a single cube.
2457 """
2458 _spatial_plot(
2459 "pcolormesh",
2460 cube,
2461 filename,
2462 sequence_coordinate,
2463 stamp_coordinate,
2464 overlay_cube=overlay_cube,
2465 contour_cube=contour_cube,
2466 )
2467 return cube, overlay_cube, contour_cube
2470# TODO: Expand function to handle ensemble data.
2471# line_coordinate: str, optional
2472# Coordinate about which to plot multiple lines. Defaults to
2473# ``"realization"``.
2474def plot_line_series(
2475 cube: iris.cube.Cube | iris.cube.CubeList,
2476 filename: str = None,
2477 series_coordinate: str = "time",
2478 sequence_coordinate: str = "time",
2479 # add the following for ensembles
2480 stamp_coordinate: str = "realization",
2481 single_plot: bool = False,
2482 **kwargs,
2483) -> iris.cube.Cube | iris.cube.CubeList:
2484 """Plot a line plot for the specified coordinate.
2486 The Cube or CubeList must be 1D.
2488 Parameters
2489 ----------
2490 iris.cube | iris.cube.CubeList
2491 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2492 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2493 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2494 filename: str, optional
2495 Name of the plot to write, used as a prefix for plot sequences. Defaults
2496 to the recipe name.
2497 series_coordinate: str, optional
2498 Coordinate about which to make a series. Defaults to ``"time"``. This
2499 coordinate must exist in the cube.
2501 Returns
2502 -------
2503 iris.cube.Cube | iris.cube.CubeList
2504 The original Cube or CubeList (so further operations can be applied).
2505 plotted data.
2507 Raises
2508 ------
2509 ValueError
2510 If the cubes don't have the right dimensions.
2511 TypeError
2512 If the cube isn't a Cube or CubeList.
2513 """
2514 stamp_coordinate = "realization"
2515 # Ensure we have a name for the plot file.
2516 recipe_title = get_recipe_metadata().get("title", "Untitled")
2518 num_models = _get_num_models(cube)
2520 _validate_cube_shape(cube, num_models)
2522 # Iterate over all cubes and extract coordinate to plot.
2523 cubes = iter_maybe(cube)
2525 coords = []
2526 for cube in cubes:
2527 try:
2528 coords.append(cube.coord(series_coordinate))
2529 except iris.exceptions.CoordinateNotFoundError as err:
2530 raise ValueError(
2531 f"Cube must have a {series_coordinate} coordinate."
2532 ) from err
2533 if cube.coords("realization"): 2533 ↛ 2537line 2533 didn't jump to line 2537 because the condition on line 2533 was always true
2534 if cube.ndim > 3: 2534 ↛ 2535line 2534 didn't jump to line 2535 because the condition on line 2534 was never true
2535 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2536 else:
2537 raise ValueError("Cube must have a realization coordinate.")
2539 plot_index = []
2541 # Check if this is a spectral plot by looking for spectral coordinates
2542 is_spectral_plot = series_coordinate in [
2543 "frequency",
2544 "physical_wavenumber",
2545 "wavelength",
2546 ]
2548 if is_spectral_plot:
2549 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2550 # coordinate frequency/wavenumber.
2551 # If several power spectra are plotted with time as sequence_coordinate for the
2552 # time slider option.
2554 # Internal plotting function.
2555 plotting_func = _plot_and_save_line_power_spectrum_series
2557 for cube in cubes:
2558 try:
2559 cube.coord(sequence_coordinate)
2560 except iris.exceptions.CoordinateNotFoundError as err:
2561 raise ValueError(
2562 f"Cube must have a {sequence_coordinate} coordinate."
2563 ) from err
2565 if num_models == 1: 2565 ↛ 2579line 2565 didn't jump to line 2579 because the condition on line 2565 was always true
2566 # check for ensembles
2567 if ( 2567 ↛ 2571line 2567 didn't jump to line 2571 because the condition on line 2567 was never true
2568 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2569 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2570 ):
2571 if single_plot:
2572 # Plot spectra, mean and ensemble spread on 1 plot
2573 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2574 else:
2575 # Plot postage stamps
2576 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2577 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2578 else:
2579 all_points = sorted(
2580 set(
2581 itertools.chain.from_iterable(
2582 cb.coord(sequence_coordinate).points for cb in cubes
2583 )
2584 )
2585 )
2586 all_slices = list(
2587 itertools.chain.from_iterable(
2588 cb.slices_over(sequence_coordinate) for cb in cubes
2589 )
2590 )
2591 # Matched slices (matched by seq coord point; it may happen that
2592 # evaluated models do not cover the same seq coord range, hence matching
2593 # necessary)
2594 cube_iterables = [
2595 iris.cube.CubeList(
2596 s
2597 for s in all_slices
2598 if s.coord(sequence_coordinate).points[0] == point
2599 )
2600 for point in all_points
2601 ]
2603 nplot = np.size(cube.coord(sequence_coordinate).points)
2605 # Create a plot for each value of the sequence coordinate. Allowing for
2606 # multiple cubes in a CubeList to be plotted in the same plot for similar
2607 # sequence values. Passing a CubeList into the internal plotting function
2608 # for similar values of the sequence coordinate. cube_slice can be an
2609 # iris.cube.Cube or an iris.cube.CubeList.
2611 for cube_slice in cube_iterables:
2612 # Normalize cube_slice to a list of cubes
2613 if isinstance(cube_slice, iris.cube.CubeList): 2613 ↛ 2614line 2613 didn't jump to line 2614 because the condition on line 2613 was never true
2614 cubes = list(cube_slice)
2615 elif isinstance(cube_slice, iris.cube.Cube): 2615 ↛ 2618line 2615 didn't jump to line 2618 because the condition on line 2615 was always true
2616 cubes = [cube_slice]
2617 else:
2618 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2620 # Use sequence value so multiple sequences can merge.
2621 seq_coord = cube_slice[0].coord(sequence_coordinate)
2622 plot_title, plot_filename = _set_title_and_filename(
2623 seq_coord, nplot, recipe_title, filename
2624 )
2626 # Format the coordinate value in a unit appropriate way.
2627 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2629 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2630 if nplot == 1 and seq_coord.has_bounds: 2630 ↛ 2635line 2630 didn't jump to line 2635 because the condition on line 2630 was always true
2631 if np.size(seq_coord.bounds) > 1: 2631 ↛ 2632line 2631 didn't jump to line 2632 because the condition on line 2631 was never true
2632 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2634 # Do the actual plotting.
2635 plotting_func(
2636 cube_slice,
2637 coords,
2638 stamp_coordinate,
2639 plot_filename,
2640 title,
2641 series_coordinate,
2642 )
2644 plot_index.append(plot_filename)
2645 else:
2646 # Format the title and filename using plotted series coordinate
2647 nplot = 1
2648 seq_coord = coords[0]
2649 plot_title, plot_filename = _set_title_and_filename(
2650 seq_coord, nplot, recipe_title, filename
2651 )
2652 # Do the actual plotting for all other series coordinate options.
2653 _plot_and_save_line_series(
2654 cubes, coords, stamp_coordinate, plot_filename, plot_title
2655 )
2657 plot_index.append(plot_filename)
2659 # append plot to list of plots
2660 complete_plot_index = _append_to_plot_index(plot_index)
2662 # Make a page to display the plots.
2663 _make_plot_html_page(complete_plot_index)
2665 return cube
2668def plot_vertical_line_series(
2669 cubes: iris.cube.Cube | iris.cube.CubeList,
2670 filename: str = None,
2671 series_coordinate: str = "model_level_number",
2672 sequence_coordinate: str = "time",
2673 # line_coordinate: str = "realization",
2674 **kwargs,
2675) -> iris.cube.Cube | iris.cube.CubeList:
2676 """Plot a line plot against a type of vertical coordinate.
2678 The Cube or CubeList must be 1D.
2680 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2681 then a sequence of plots will be produced.
2683 Parameters
2684 ----------
2685 iris.cube | iris.cube.CubeList
2686 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2687 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2688 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2689 filename: str, optional
2690 Name of the plot to write, used as a prefix for plot sequences. Defaults
2691 to the recipe name.
2692 series_coordinate: str, optional
2693 Coordinate to plot on the y-axis. Can be ``pressure`` or
2694 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2695 for LFRic. Defaults to ``model_level_number``.
2696 This coordinate must exist in the cube.
2697 sequence_coordinate: str, optional
2698 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2699 This coordinate must exist in the cube.
2701 Returns
2702 -------
2703 iris.cube.Cube | iris.cube.CubeList
2704 The original Cube or CubeList (so further operations can be applied).
2705 Plotted data.
2707 Raises
2708 ------
2709 ValueError
2710 If the cubes doesn't have the right dimensions.
2711 TypeError
2712 If the cube isn't a Cube or CubeList.
2713 """
2714 # Ensure we have a name for the plot file.
2715 recipe_title = get_recipe_metadata().get("title", "Untitled")
2717 cubes = iter_maybe(cubes)
2718 # Initialise empty list to hold all data from all cubes in a CubeList
2719 all_data = []
2721 # Store min/max ranges for x range.
2722 x_levels = []
2724 num_models = _get_num_models(cubes)
2726 _validate_cube_shape(cubes, num_models)
2728 # Iterate over all cubes in cube or CubeList and plot.
2729 coords = []
2730 for cube in cubes:
2731 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2732 try:
2733 coords.append(cube.coord(series_coordinate))
2734 except iris.exceptions.CoordinateNotFoundError as err:
2735 raise ValueError(
2736 f"Cube must have a {series_coordinate} coordinate."
2737 ) from err
2739 try:
2740 if cube.ndim > 1 or not cube.coords("realization"): 2740 ↛ 2748line 2740 didn't jump to line 2748 because the condition on line 2740 was always true
2741 cube.coord(sequence_coordinate)
2742 except iris.exceptions.CoordinateNotFoundError as err:
2743 raise ValueError(
2744 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2745 ) from err
2747 # Get minimum and maximum from levels information.
2748 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2749 if levels is not None: 2749 ↛ 2753line 2749 didn't jump to line 2753 because the condition on line 2749 was always true
2750 x_levels.append(min(levels))
2751 x_levels.append(max(levels))
2752 else:
2753 all_data.append(cube.data)
2755 if len(x_levels) == 0: 2755 ↛ 2757line 2755 didn't jump to line 2757 because the condition on line 2755 was never true
2756 # Combine all data into a single NumPy array
2757 combined_data = np.concatenate(all_data)
2759 # Set the lower and upper limit for the x-axis to ensure all plots have
2760 # same range. This needs to read the whole cube over the range of the
2761 # sequence and if applicable postage stamp coordinate.
2762 vmin = np.floor(combined_data.min())
2763 vmax = np.ceil(combined_data.max())
2764 else:
2765 vmin = min(x_levels)
2766 vmax = max(x_levels)
2768 # Matching the slices (matching by seq coord point; it may happen that
2769 # evaluated models do not cover the same seq coord range, hence matching
2770 # necessary)
2771 def filter_cube_iterables(cube_iterables) -> bool:
2772 return len(cube_iterables) == len(coords)
2774 cube_iterables = filter(
2775 filter_cube_iterables,
2776 (
2777 iris.cube.CubeList(
2778 s
2779 for s in itertools.chain.from_iterable(
2780 cb.slices_over(sequence_coordinate) for cb in cubes
2781 )
2782 if s.coord(sequence_coordinate).points[0] == point
2783 )
2784 for point in sorted(
2785 set(
2786 itertools.chain.from_iterable(
2787 cb.coord(sequence_coordinate).points for cb in cubes
2788 )
2789 )
2790 )
2791 ),
2792 )
2794 # Create a plot for each value of the sequence coordinate.
2795 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2796 # similar sequence values. Passing a CubeList into the internal plotting function
2797 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2798 # or an iris.cube.CubeList.
2799 plot_index = []
2800 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2801 for cubes_slice in cube_iterables:
2802 # Format the coordinate value in a unit appropriate way.
2803 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2804 plot_title, plot_filename = _set_title_and_filename(
2805 seq_coord, nplot, recipe_title, filename
2806 )
2808 # Do the actual plotting.
2809 _plot_and_save_vertical_line_series(
2810 cubes_slice,
2811 coords,
2812 "realization",
2813 plot_filename,
2814 series_coordinate,
2815 title=plot_title,
2816 vmin=vmin,
2817 vmax=vmax,
2818 )
2819 plot_index.append(plot_filename)
2821 # Add list of plots to plot metadata.
2822 complete_plot_index = _append_to_plot_index(plot_index)
2824 # Make a page to display the plots.
2825 _make_plot_html_page(complete_plot_index)
2827 return cubes
2830def qq_plot(
2831 cubes: iris.cube.CubeList,
2832 coordinates: list[str],
2833 percentiles: list[float],
2834 model_names: list[str],
2835 filename: str = None,
2836 one_to_one: bool = True,
2837 **kwargs,
2838) -> iris.cube.CubeList:
2839 """Plot a Quantile-Quantile plot between two models for common time points.
2841 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2842 collapsed within the operator over all specified coordinates such as
2843 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2844 ensemble members to ensure a 1D cube (array).
2846 Parameters
2847 ----------
2848 cubes: iris.cube.CubeList
2849 Two cubes of the same variable with different models.
2850 coordinate: list[str]
2851 The list of coordinates to collapse over. This list should be
2852 every coordinate within the cube to result in a 1D cube around
2853 the percentile coordinate.
2854 percent: list[float]
2855 A list of percentiles to appear in the plot.
2856 model_names: list[str]
2857 A list of model names to appear on the axis of the plot.
2858 filename: str, optional
2859 Filename of the plot to write.
2860 one_to_one: bool, optional
2861 If True a 1:1 line is plotted; if False it is not. Default is True.
2863 Raises
2864 ------
2865 ValueError
2866 When the cubes are not compatible.
2868 Notes
2869 -----
2870 The quantile-quantile plot is a variant on the scatter plot representing
2871 two datasets by their quantiles (percentiles) for common time points.
2872 This plot does not use a theoretical distribution to compare against, but
2873 compares percentiles of two datasets. This plot does
2874 not use all raw data points, but plots the selected percentiles (quantiles) of
2875 each variable instead for the two datasets, thereby normalising the data for a
2876 direct comparison between the selected percentiles of the two dataset distributions.
2878 Quantile-quantile plots are valuable for comparing against
2879 observations and other models. Identical percentiles between the variables
2880 will lie on the one-to-one line implying the values correspond well to each
2881 other. Where there is a deviation from the one-to-one line a range of
2882 possibilities exist depending on how and where the data is shifted (e.g.,
2883 Wilks 2011 [Wilks2011]_).
2885 For distributions above the one-to-one line the distribution is left-skewed;
2886 below is right-skewed. A distinct break implies a bimodal distribution, and
2887 closer values/values further apart at the tails imply poor representation of
2888 the extremes.
2890 References
2891 ----------
2892 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2893 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2894 """
2895 # Check cubes using same functionality as the difference operator.
2896 if len(cubes) != 2:
2897 raise ValueError("cubes should contain exactly 2 cubes.")
2898 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2899 other: Cube = cubes.extract_cube(
2900 iris.Constraint(
2901 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2902 )
2903 )
2905 # Get spatial coord names.
2906 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2907 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2909 # Ensure cubes to compare are on common differencing grid.
2910 # This is triggered if either
2911 # i) latitude and longitude shapes are not the same. Note grid points
2912 # are not compared directly as these can differ through rounding
2913 # errors.
2914 # ii) or variables are known to often sit on different grid staggering
2915 # in different models (e.g. cell center vs cell edge), as is the case
2916 # for UM and LFRic comparisons.
2917 # In future greater choice of regridding method might be applied depending
2918 # on variable type. Linear regridding can in general be appropriate for smooth
2919 # variables. Care should be taken with interpretation of differences
2920 # given this dependency on regridding.
2921 if (
2922 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2923 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2924 ) or (
2925 base.long_name
2926 in [
2927 "eastward_wind_at_10m",
2928 "northward_wind_at_10m",
2929 "northward_wind_at_cell_centres",
2930 "eastward_wind_at_cell_centres",
2931 "zonal_wind_at_pressure_levels",
2932 "meridional_wind_at_pressure_levels",
2933 "potential_vorticity_at_pressure_levels",
2934 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2935 ]
2936 ):
2937 logging.debug(
2938 "Linear regridding base cube to other grid to compute differences"
2939 )
2940 base = regrid_onto_cube(base, other, method="Linear")
2942 # Extract just common time points.
2943 base, other = _extract_common_time_points(base, other)
2945 # Equalise attributes so we can merge.
2946 fully_equalise_attributes([base, other])
2947 logging.debug("Base: %s\nOther: %s", base, other)
2949 # Collapse cubes.
2950 base = collapse(
2951 base,
2952 coordinate=coordinates,
2953 method="PERCENTILE",
2954 additional_percent=percentiles,
2955 )
2956 other = collapse(
2957 other,
2958 coordinate=coordinates,
2959 method="PERCENTILE",
2960 additional_percent=percentiles,
2961 )
2963 # Ensure we have a name for the plot file.
2964 recipe_title = get_recipe_metadata().get("title", "Untitled")
2965 title = f"{recipe_title}"
2967 if filename is None:
2968 filename = slugify(recipe_title)
2970 # Add file extension.
2971 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2973 # Do the actual plotting on a scatter plot
2974 _plot_and_save_scatter_plot(
2975 base, other, plot_filename, title, one_to_one, model_names
2976 )
2978 # Add list of plots to plot metadata.
2979 plot_index = _append_to_plot_index([plot_filename])
2981 # Make a page to display the plots.
2982 _make_plot_html_page(plot_index)
2984 return iris.cube.CubeList([base, other])
2987def scatter_plot(
2988 cube_x: iris.cube.Cube | iris.cube.CubeList,
2989 cube_y: iris.cube.Cube | iris.cube.CubeList,
2990 filename: str = None,
2991 one_to_one: bool = True,
2992 **kwargs,
2993) -> iris.cube.CubeList:
2994 """Plot a scatter plot between two variables.
2996 Both cubes must be 1D.
2998 Parameters
2999 ----------
3000 cube_x: Cube | CubeList
3001 1 dimensional Cube of the data to plot on y-axis.
3002 cube_y: Cube | CubeList
3003 1 dimensional Cube of the data to plot on x-axis.
3004 filename: str, optional
3005 Filename of the plot to write.
3006 one_to_one: bool, optional
3007 If True a 1:1 line is plotted; if False it is not. Default is True.
3009 Returns
3010 -------
3011 cubes: CubeList
3012 CubeList of the original x and y cubes for further processing.
3014 Raises
3015 ------
3016 ValueError
3017 If the cube doesn't have the right dimensions and cubes not the same
3018 size.
3019 TypeError
3020 If the cube isn't a single cube.
3022 Notes
3023 -----
3024 Scatter plots are used for determining if there is a relationship between
3025 two variables. Positive relations have a slope going from bottom left to top
3026 right; Negative relations have a slope going from top left to bottom right.
3027 """
3028 # Iterate over all cubes in cube or CubeList and plot.
3029 for cube_iter in iter_maybe(cube_x):
3030 # Check cubes are correct shape.
3031 cube_iter = _check_single_cube(cube_iter)
3032 if cube_iter.ndim > 1:
3033 raise ValueError("cube_x must be 1D.")
3035 # Iterate over all cubes in cube or CubeList and plot.
3036 for cube_iter in iter_maybe(cube_y):
3037 # Check cubes are correct shape.
3038 cube_iter = _check_single_cube(cube_iter)
3039 if cube_iter.ndim > 1:
3040 raise ValueError("cube_y must be 1D.")
3042 # Ensure we have a name for the plot file.
3043 recipe_title = get_recipe_metadata().get("title", "Untitled")
3044 title = f"{recipe_title}"
3046 if filename is None:
3047 filename = slugify(recipe_title)
3049 # Add file extension.
3050 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3052 # Do the actual plotting.
3053 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3055 # Add list of plots to plot metadata.
3056 plot_index = _append_to_plot_index([plot_filename])
3058 # Make a page to display the plots.
3059 _make_plot_html_page(plot_index)
3061 return iris.cube.CubeList([cube_x, cube_y])
3064def vector_plot(
3065 cube_u: iris.cube.Cube,
3066 cube_v: iris.cube.Cube,
3067 filename: str = None,
3068 sequence_coordinate: str = "time",
3069 **kwargs,
3070) -> iris.cube.CubeList:
3071 """Plot a vector plot based on the input u and v components."""
3072 recipe_title = get_recipe_metadata().get("title", "Untitled")
3074 # Cubes must have a matching sequence coordinate.
3075 try:
3076 # Check that the u and v cubes have the same sequence coordinate.
3077 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3077 ↛ anywhereline 3077 didn't jump anywhere: it always raised an exception.
3078 raise ValueError("Coordinates do not match.")
3079 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3080 raise ValueError(
3081 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3082 ) from err
3084 # Create a plot for each value of the sequence coordinate.
3085 plot_index = []
3086 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3087 for cube_u_slice, cube_v_slice in zip(
3088 cube_u.slices_over(sequence_coordinate),
3089 cube_v.slices_over(sequence_coordinate),
3090 strict=True,
3091 ):
3092 # Format the coordinate value in a unit appropriate way.
3093 seq_coord = cube_u_slice.coord(sequence_coordinate)
3094 plot_title, plot_filename = _set_title_and_filename(
3095 seq_coord, nplot, recipe_title, filename
3096 )
3098 # Do the actual plotting.
3099 _plot_and_save_vector_plot(
3100 cube_u_slice,
3101 cube_v_slice,
3102 filename=plot_filename,
3103 title=plot_title,
3104 method="contourf",
3105 )
3106 plot_index.append(plot_filename)
3108 # Add list of plots to plot metadata.
3109 complete_plot_index = _append_to_plot_index(plot_index)
3111 # Make a page to display the plots.
3112 _make_plot_html_page(complete_plot_index)
3114 return iris.cube.CubeList([cube_u, cube_v])
3117def plot_histogram_series(
3118 cubes: iris.cube.Cube | iris.cube.CubeList,
3119 filename: str = None,
3120 sequence_coordinate: str = "time",
3121 stamp_coordinate: str = "realization",
3122 single_plot: bool = False,
3123 **kwargs,
3124) -> iris.cube.Cube | iris.cube.CubeList:
3125 """Plot a histogram plot for each vertical level provided.
3127 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3128 is present then a sequence of plots will be produced using the time slider
3129 functionality to scroll through histograms against time. If a
3130 stamp_coordinate is present then postage stamp plots will be produced. If
3131 stamp_coordinate and single_plot is True, all postage stamp plots will be
3132 plotted in a single plot instead of separate postage stamp plots.
3134 Parameters
3135 ----------
3136 cubes: Cube | iris.cube.CubeList
3137 Iris cube or CubeList of the data to plot. It should have a single dimension other
3138 than the stamp coordinate.
3139 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3140 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3141 filename: str, optional
3142 Name of the plot to write, used as a prefix for plot sequences. Defaults
3143 to the recipe name.
3144 sequence_coordinate: str, optional
3145 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3146 This coordinate must exist in the cube and will be used for the time
3147 slider.
3148 stamp_coordinate: str, optional
3149 Coordinate about which to plot postage stamp plots. Defaults to
3150 ``"realization"``.
3151 single_plot: bool, optional
3152 If True, all postage stamp plots will be plotted in a single plot. If
3153 False, each postage stamp plot will be plotted separately. Is only valid
3154 if stamp_coordinate exists and has more than a single point.
3156 Returns
3157 -------
3158 iris.cube.Cube | iris.cube.CubeList
3159 The original Cube or CubeList (so further operations can be applied).
3160 Plotted data.
3162 Raises
3163 ------
3164 ValueError
3165 If the cube doesn't have the right dimensions.
3166 TypeError
3167 If the cube isn't a Cube or CubeList.
3168 """
3169 recipe_title = get_recipe_metadata().get("title", "Untitled")
3171 cubes = iter_maybe(cubes)
3172 # Ensure we have a name for the plot file.
3173 if filename is None:
3174 filename = slugify(recipe_title)
3176 # Internal plotting function.
3177 plotting_func = _plot_and_save_histogram_series
3179 num_models = _get_num_models(cubes)
3181 _validate_cube_shape(cubes, num_models)
3183 # If several histograms are plotted with time as sequence_coordinate for the
3184 # time slider option.
3185 for cube in cubes:
3186 try:
3187 cube.coord(sequence_coordinate)
3188 except iris.exceptions.CoordinateNotFoundError as err:
3189 raise ValueError(
3190 f"Cube must have a {sequence_coordinate} coordinate."
3191 ) from err
3193 # Get minimum and maximum from levels information.
3194 levels = None
3195 for cube in cubes: 3195 ↛ 3211line 3195 didn't jump to line 3211 because the loop on line 3195 didn't complete
3196 # First check if user-specified "auto" range variable.
3197 # This maintains the value of levels as None, so proceed.
3198 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3199 if levels is None:
3200 break
3201 # If levels is changed, recheck to use the vmin,vmax or
3202 # levels-based ranges for histogram plots.
3203 _, levels, _ = _colorbar_map_levels(cube)
3204 logging.debug("levels: %s", levels)
3205 if levels is not None: 3205 ↛ 3195line 3205 didn't jump to line 3195 because the condition on line 3205 was always true
3206 vmin = min(levels)
3207 vmax = max(levels)
3208 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
3209 break
3211 if levels is None:
3212 vmin = min(cb.data.min() for cb in cubes)
3213 vmax = max(cb.data.max() for cb in cubes)
3215 # Make postage stamp plots if stamp_coordinate exists and has more than a
3216 # single point. If single_plot is True:
3217 # -- all postage stamp plots will be plotted in a single plot instead of
3218 # separate postage stamp plots.
3219 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3220 # produced per single model only
3222 if num_models == 1: 3222 ↛ 3235line 3222 didn't jump to line 3235 because the condition on line 3222 was always true
3223 if ( 3223 ↛ 3227line 3223 didn't jump to line 3227 because the condition on line 3223 was never true
3224 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3225 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3226 ):
3227 if single_plot:
3228 plotting_func = (
3229 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3230 )
3231 else:
3232 plotting_func = _plot_and_save_postage_stamp_histogram_series
3233 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3234 else:
3235 all_points = sorted(
3236 set(
3237 itertools.chain.from_iterable(
3238 cb.coord(sequence_coordinate).points for cb in cubes
3239 )
3240 )
3241 )
3242 all_slices = list(
3243 itertools.chain.from_iterable(
3244 cb.slices_over(sequence_coordinate) for cb in cubes
3245 )
3246 )
3247 # Matched slices (matched by seq coord point; it may happen that
3248 # evaluated models do not cover the same seq coord range, hence matching
3249 # necessary)
3250 cube_iterables = [
3251 iris.cube.CubeList(
3252 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3253 )
3254 for point in all_points
3255 ]
3257 plot_index = []
3258 nplot = np.size(cube.coord(sequence_coordinate).points)
3259 # Create a plot for each value of the sequence coordinate. Allowing for
3260 # multiple cubes in a CubeList to be plotted in the same plot for similar
3261 # sequence values. Passing a CubeList into the internal plotting function
3262 # for similar values of the sequence coordinate. cube_slice can be an
3263 # iris.cube.Cube or an iris.cube.CubeList.
3264 for cube_slice in cube_iterables:
3265 single_cube = cube_slice
3266 if isinstance(cube_slice, iris.cube.CubeList): 3266 ↛ 3267line 3266 didn't jump to line 3267 because the condition on line 3266 was never true
3267 single_cube = cube_slice[0]
3269 # Set plot titles and filename, based on sequence coordinate
3270 seq_coord = single_cube.coord(sequence_coordinate)
3271 # Use time coordinate in title and filename if single histogram output.
3272 if sequence_coordinate == "realization" and nplot == 1: 3272 ↛ 3273line 3272 didn't jump to line 3273 because the condition on line 3272 was never true
3273 seq_coord = single_cube.coord("time")
3274 plot_title, plot_filename = _set_title_and_filename(
3275 seq_coord, nplot, recipe_title, filename
3276 )
3278 # Do the actual plotting.
3280 plotting_func(
3281 cube_slice,
3282 filename=plot_filename,
3283 stamp_coordinate=stamp_coordinate,
3284 title=plot_title,
3285 vmin=vmin,
3286 vmax=vmax,
3287 )
3288 plot_index.append(plot_filename)
3290 # Add list of plots to plot metadata.
3291 complete_plot_index = _append_to_plot_index(plot_index)
3293 # Make a page to display the plots.
3294 _make_plot_html_page(complete_plot_index)
3296 return cubes
3299def _plot_and_save_postage_stamp_power_spectrum_series(
3300 cubes: iris.cube.Cube,
3301 coords: list[iris.coords.Coord],
3302 stamp_coordinate: str,
3303 filename: str,
3304 title: str,
3305 series_coordinate: str = None,
3306 **kwargs,
3307):
3308 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3310 Parameters
3311 ----------
3312 cubes: Cube or CubeList
3313 Cube or Cubelist of the power spectrum data.
3314 coords: list[Coord]
3315 Coordinates to plot on the x-axis, one per cube.
3316 stamp_coordinate: str
3317 Coordinate that becomes different plots.
3318 filename: str
3319 Filename of the plot to write.
3320 title: str
3321 Plot title.
3322 series_coordinate: str, optional
3323 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3325 """
3326 # Use the smallest square grid that will fit the members.
3327 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
3329 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3330 model_colors_map = _get_model_colors_map(cubes)
3331 # ax = plt.gca()
3332 # Make a subplot for each member.
3333 for member, subplot in zip(
3334 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3335 ):
3336 ax = plt.subplot(grid_size, grid_size, subplot)
3338 # Store min/max ranges.
3339 y_levels = []
3341 line_marker = None
3342 line_width = 1
3344 for cube in iter_maybe(member):
3345 xcoord = select_series_coord(cube, series_coordinate)
3346 xname = xcoord.points
3348 yfield = cube.data # power spectrum
3349 label = None
3350 color = "black"
3351 if model_colors_map: 3351 ↛ 3352line 3351 didn't jump to line 3352 because the condition on line 3351 was never true
3352 label = cube.attributes.get("model_name")
3353 color = model_colors_map.get(label)
3355 if member.coord(stamp_coordinate).points == [0]:
3356 ax.plot(
3357 xname,
3358 yfield,
3359 color=color,
3360 marker=line_marker,
3361 ls="-",
3362 lw=line_width,
3363 label=f"{label} (control)"
3364 if len(cube.coord(stamp_coordinate).points) > 1
3365 else label,
3366 )
3367 # Label with member if part of an ensemble and not the control.
3368 else:
3369 ax.plot(
3370 xname,
3371 yfield,
3372 color=color,
3373 ls="-",
3374 lw=1.5,
3375 alpha=0.75,
3376 label=f"{label} (member)",
3377 )
3379 # Calculate the global min/max if multiple cubes are given.
3380 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3381 if levels is not None: 3381 ↛ 3382line 3381 didn't jump to line 3382 because the condition on line 3381 was never true
3382 y_levels.append(min(levels))
3383 y_levels.append(max(levels))
3385 # Add some labels and tweak the style.
3386 title = f"{title}"
3387 ax.set_title(title, fontsize=16)
3389 # Set appropriate x-axis label based on coordinate
3390 if series_coordinate == "wavelength" or ( 3390 ↛ 3393line 3390 didn't jump to line 3393 because the condition on line 3390 was never true
3391 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3392 ):
3393 ax.set_xlabel("Wavelength (km)", fontsize=14)
3394 elif series_coordinate == "physical_wavenumber" or ( 3394 ↛ 3399line 3394 didn't jump to line 3399 because the condition on line 3394 was always true
3395 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3396 ):
3397 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3398 else: # frequency or check units
3399 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3400 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3401 else:
3402 ax.set_xlabel("Wavenumber", fontsize=14)
3404 ax.set_ylabel("Power Spectral Density", fontsize=14)
3405 ax.tick_params(axis="both", labelsize=12)
3407 # Set log-log scale
3408 ax.set_xscale("log")
3409 ax.set_yscale("log")
3411 # Add gridlines
3412 ax.grid(linestyle="--", color="grey", linewidth=1)
3413 # Ientify unique labels for legend
3414 handles = list(
3415 {
3416 label: handle
3417 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3418 }.values()
3419 )
3420 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
3422 ax = plt.gca()
3423 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3425 # Overall figure title.
3426 fig.suptitle(title, fontsize=16)
3428 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3429 logging.info("Saved histogram postage stamp plot to %s", filename)
3430 plt.close(fig)
3433def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3434 cubes: iris.cube.Cube,
3435 coords: list[iris.coords.Coord],
3436 stamp_coordinate: str,
3437 filename: str,
3438 title: str,
3439 series_coordinate: str = None,
3440 **kwargs,
3441):
3442 """Plot and save power spectra for ensemble members in single plot.
3444 Parameters
3445 ----------
3446 cubes: Cube or CubeList
3447 Cube or Cubelist of the power spectrum data.
3448 coords: list[Coord]
3449 Coordinates to plot on the x-axis, one per cube.
3450 stamp_coordinate: str
3451 Coordinate that becomes different plots.
3452 filename: str
3453 Filename of the plot to write.
3454 title: str
3455 Plot title.
3456 series_coordinate: str, optional
3457 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3459 """
3460 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3461 model_colors_map = _get_model_colors_map(cubes)
3463 line_marker = None
3464 line_width = 1
3466 # Compute ensemble statistics to show spread
3467 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3468 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3469 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3471 xcoord_global = mean_cube.coord(series_coordinate)
3472 x_global = xcoord_global.points
3474 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3475 xcoord = select_series_coord(member, series_coordinate)
3476 xname = xcoord.points
3478 yfield = member.data # power spectrum
3479 color = "black"
3480 if model_colors_map: 3480 ↛ 3484line 3480 didn't jump to line 3484 because the condition on line 3480 was always true
3481 label = member.attributes.get("model_name") if i == 0 else None
3482 color = model_colors_map.get(label)
3484 if member.coord(stamp_coordinate).points == [0]:
3485 ax.plot(
3486 xname,
3487 yfield,
3488 color=color,
3489 marker=line_marker,
3490 ls="-",
3491 lw=line_width,
3492 label=f"{label} (control)"
3493 if len(member.coord(stamp_coordinate).points) > 1
3494 else label,
3495 )
3496 # Label with member number if part of an ensemble and not the control.
3497 else:
3498 ax.plot(
3499 xname,
3500 yfield,
3501 color=color,
3502 ls="-",
3503 lw=1.5,
3504 alpha=0.75,
3505 label=label,
3506 )
3508 # Set appropriate x-axis label based on coordinate
3509 if series_coordinate == "wavelength" or ( 3509 ↛ 3512line 3509 didn't jump to line 3512 because the condition on line 3509 was never true
3510 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3511 ):
3512 ax.set_xlabel("Wavelength (km)", fontsize=14)
3513 elif series_coordinate == "physical_wavenumber" or ( 3513 ↛ 3518line 3513 didn't jump to line 3518 because the condition on line 3513 was always true
3514 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3515 ):
3516 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3517 else: # frequency or check units
3518 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3519 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3520 else:
3521 ax.set_xlabel("Wavenumber", fontsize=14)
3523 # Add ensemble spread shading
3524 ax.fill_between(
3525 x_global,
3526 min_cube.data,
3527 max_cube.data,
3528 color="grey",
3529 alpha=0.3,
3530 label="Ensemble spread",
3531 )
3533 # Add ensemble mean line
3534 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3536 ax.set_ylabel("Power Spectral Density", fontsize=14)
3537 ax.tick_params(axis="both", labelsize=12)
3539 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3540 # Set log-log scale
3541 ax.set_xscale("log")
3542 ax.set_yscale("log")
3544 # Add gridlines
3545 ax.grid(linestyle="--", color="grey", linewidth=1)
3546 # Identify unique labels for legend
3547 handles = list(
3548 {
3549 label: handle
3550 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3551 }.values()
3552 )
3553 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
3555 # Add a legend
3556 ax.legend(fontsize=16)
3558 # Figure title.
3559 ax.set_title(title, fontsize=16)
3561 # Save the figure to a file
3562 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3564 # Close the figure
3565 plt.close(fig)