Coverage for src/CSET/operators/plot.py: 83%
1092 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-16 13:58 +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
38import scipy.fft as fft
39from cartopy.mpl.geoaxes import GeoAxes
40from iris.cube import Cube
41from markdown_it import MarkdownIt
42from mpl_toolkits.axes_grid1.inset_locator import inset_axes
44from CSET._common import (
45 combine_dicts,
46 filename_slugify,
47 get_recipe_metadata,
48 iter_maybe,
49 render_file,
50 slugify,
51)
52from CSET.operators._utils import (
53 check_stamp_coordinate,
54 fully_equalise_attributes,
55 get_cube_yxcoordname,
56 is_transect,
57 slice_over_maybe,
58)
59from CSET.operators.collapse import collapse
60from CSET.operators.misc import _extract_common_time_points
61from CSET.operators.regrid import regrid_onto_cube
63# Use a non-interactive plotting backend.
64mpl.use("agg")
66DEFAULT_DISCRETE_COLORS = mpl.colormaps["tab10"].colors + mpl.colormaps["Accent"].colors
68############################
69# Private helper functions #
70############################
73def _append_to_plot_index(plot_index: list) -> list:
74 """Add plots into the plot index, returning the complete plot index."""
75 with open("meta.json", "r+t", encoding="UTF-8") as fp:
76 fcntl.flock(fp, fcntl.LOCK_EX)
77 fp.seek(0)
78 meta = json.load(fp)
79 complete_plot_index = meta.get("plots", [])
80 complete_plot_index = complete_plot_index + plot_index
81 meta["plots"] = complete_plot_index
82 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
83 os.getenv("DO_CASE_AGGREGATION")
84 ):
85 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
86 fp.seek(0)
87 fp.truncate()
88 json.dump(meta, fp, indent=2)
89 return complete_plot_index
92def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:
93 """Ensure a single cube is given.
95 If a CubeList of length one is given that the contained cube is returned,
96 otherwise an error is raised.
98 Parameters
99 ----------
100 cube: Cube | CubeList
101 The cube to check.
103 Returns
104 -------
105 cube: Cube
106 The checked cube.
108 Raises
109 ------
110 TypeError
111 If the input cube is not a Cube or CubeList of a single Cube.
112 """
113 if isinstance(cube, iris.cube.Cube):
114 return cube
115 if isinstance(cube, iris.cube.CubeList):
116 if len(cube) == 1:
117 return cube[0]
118 raise TypeError("Must have a single cube", cube)
121def _make_plot_html_page(plots: list):
122 """Create a HTML page to display a plot image."""
123 # Debug check that plots actually contains some strings.
124 assert isinstance(plots[0], str)
126 # Load HTML template file.
127 operator_files = importlib.resources.files()
128 template_file = operator_files.joinpath("_plot_page_template.html")
130 # Get some metadata.
131 meta = get_recipe_metadata()
132 title = meta.get("title", "Untitled")
133 description = MarkdownIt().render(meta.get("description", "*No description.*"))
135 # Prepare template variables.
136 variables = {
137 "title": title,
138 "description": description,
139 "initial_plot": plots[0],
140 "plots": plots,
141 "title_slug": slugify(title),
142 }
144 # Render template.
145 html = render_file(template_file, **variables)
147 # Save completed HTML.
148 with open("index.html", "wt", encoding="UTF-8") as fp:
149 fp.write(html)
152@functools.cache
153def _load_colorbar_map(user_colorbar_file: str = None) -> dict:
154 """Load the colorbar definitions from a file.
156 This is a separate function to make it cacheable.
157 """
158 colorbar_file = importlib.resources.files().joinpath("_colorbar_definition.json")
159 with open(colorbar_file, "rt", encoding="UTF-8") as fp:
160 colorbar = json.load(fp)
162 logging.debug("User colour bar file: %s", user_colorbar_file)
163 override_colorbar = {}
164 if user_colorbar_file:
165 try:
166 with open(user_colorbar_file, "rt", encoding="UTF-8") as fp:
167 override_colorbar = json.load(fp)
168 except FileNotFoundError:
169 logging.warning("Colorbar file does not exist. Using default values.")
171 # Overwrite values with the user supplied colorbar definition.
172 colorbar = combine_dicts(colorbar, override_colorbar)
173 return colorbar
176def _get_model_colors_map(cubes: iris.cube.CubeList | iris.cube.Cube) -> dict:
177 """Get an appropriate colors for model lines in line plots.
179 For each model in the list of cubes colors either from user provided
180 color definition file (so-called style file) or from default colors are mapped
181 to model_name attribute.
183 Parameters
184 ----------
185 cubes: CubeList or Cube
186 Cubes with model_name attribute
188 Returns
189 -------
190 model_colors_map:
191 Dictionary mapping model_name attribute to colors
192 """
193 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
194 colorbar = _load_colorbar_map(user_colorbar_file)
195 model_names = sorted(
196 filter(
197 lambda x: x is not None,
198 (cube.attributes.get("model_name", None) for cube in iter_maybe(cubes)),
199 )
200 )
201 if not model_names:
202 return {}
203 use_user_colors = all(mname in colorbar.keys() for mname in model_names)
204 if use_user_colors: 204 ↛ 205line 204 didn't jump to line 205 because the condition on line 204 was never true
205 return {mname: colorbar[mname] for mname in model_names}
207 color_list = itertools.cycle(DEFAULT_DISCRETE_COLORS)
208 return {mname: color for mname, color in zip(model_names, color_list, strict=False)}
211def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
212 """Get an appropriate colorbar for the given cube.
214 For the given variable the appropriate colorbar is looked up from a
215 combination of the built-in CSET colorbar definitions, and any user supplied
216 definitions. As well as varying on variables, these definitions may also
217 exist for specific pressure levels to account for variables with
218 significantly different ranges at different heights. The colorbars also exist
219 for masks and mask differences for considering variable presence diagnostics.
220 Specific variable ranges can be separately set in user-supplied definition
221 for x- or y-axis limits, or indicate where automated range preferred.
223 Parameters
224 ----------
225 cube: Cube
226 Cube of variable for which the colorbar information is desired.
227 axis: "x", "y", optional
228 Select the levels for just this axis of a line plot. The min and max
229 can be set by xmin/xmax or ymin/ymax respectively. For variables where
230 setting a universal range is not desirable (e.g. temperature), users
231 can set ymin/ymax values to "auto" in the colorbar definitions file.
232 Where no additional xmin/xmax or ymin/ymax values are provided, the
233 axis bounds default to use the vmin/vmax values provided.
235 Returns
236 -------
237 cmap:
238 Matplotlib colormap.
239 levels:
240 List of levels to use for plotting. For continuous plots the min and max
241 should be taken as the range.
242 norm:
243 BoundaryNorm information.
244 """
245 # Grab the colorbar file from the recipe global metadata.
246 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
247 colorbar = _load_colorbar_map(user_colorbar_file)
248 cmap = None
250 try:
251 # We assume that pressure is a scalar coordinate here.
252 pressure_level_raw = cube.coord("pressure").points[0]
253 # Ensure pressure_level is a string, as it is used as a JSON key.
254 pressure_level = str(int(pressure_level_raw))
255 except iris.exceptions.CoordinateNotFoundError:
256 pressure_level = None
258 # First try long name, then standard name, then var name. This order is used
259 # as long name is the one we correct between models, so it most likely to be
260 # consistent.
261 varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name]))
262 for varname in varnames:
263 # Get the colormap for this variable.
264 try:
265 var_colorbar = colorbar[varname]
266 cmap = plt.get_cmap(colorbar[varname]["cmap"], 51)
267 varname_key = varname
268 break
269 except KeyError:
270 logging.debug("Cube name %s has no colorbar definition.", varname)
272 # Get colormap if it is a mask.
273 if any("mask_for_" in name for name in varnames):
274 cmap, levels, norm = _custom_colormap_mask(cube, axis=axis)
275 return cmap, levels, norm
276 # If winds on Beaufort Scale use custom colorbar and levels
277 if any("Beaufort_Scale" in name for name in varnames):
278 cmap, levels, norm = _custom_beaufort_scale(cube, axis=axis)
279 return cmap, levels, norm
280 # If probability is plotted use custom colorbar and levels
281 if any("probability_of_" in name for name in varnames):
282 cmap, levels, norm = _custom_colormap_probability(cube, axis=axis)
283 return cmap, levels, norm
284 # If aviation colour state use custom colorbar and levels
285 if any("aviation_colour_state" in name for name in varnames): 285 ↛ 286line 285 didn't jump to line 286 because the condition on line 285 was never true
286 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)
287 return cmap, levels, norm
288 if any("RMSE_" in name for name in varnames):
289 # Variables levels and norm set to None to use specific colorbar routine.
290 levels = None
291 norm = None
292 cmap, levels, norm = _custom_colormap_scores(cube, cmap, levels, norm)
293 return cmap, levels, norm
294 # If no valid colormap has been defined, use defaults and return.
295 if not cmap:
296 logging.warning("No colorbar definition exists for %s.", cube.name())
297 cmap, levels, norm = mpl.colormaps["viridis"], None, None
298 return cmap, levels, norm
300 # Test if pressure-level specific settings are provided for cube.
301 if pressure_level:
302 try:
303 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]
304 except KeyError:
305 logging.debug(
306 "%s has no colorbar definition for pressure level %s.",
307 varname,
308 pressure_level,
309 )
311 # Check for availability of x-axis or y-axis user-specific overrides
312 # for setting level bounds for line plot types and return just levels.
313 # Line plots do not need a colormap, and just use the data range.
314 if axis:
315 if axis == "x":
316 try:
317 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]
318 except KeyError:
319 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
320 if axis == "y":
321 try:
322 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]
323 except KeyError:
324 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
325 # Check if user-specified auto-scaling for this variable
326 if vmin == "auto" or vmax == "auto":
327 levels = None
328 else:
329 levels = [vmin, vmax]
330 return None, levels, None
331 # Get and use the colorbar levels for this variable if spatial or histogram.
332 else:
333 try:
334 levels = var_colorbar["levels"]
335 # Use discrete bins when levels are specified, rather
336 # than a smooth range.
337 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
338 logging.debug("Using levels for %s colorbar.", varname)
339 logging.info("Using levels: %s", levels)
340 except KeyError:
341 # Get the range for this variable.
342 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
343 logging.debug("Using min and max for %s colorbar.", varname)
344 # Calculate levels from range.
345 if vmin == "auto" or vmax == "auto": 345 ↛ 346line 345 didn't jump to line 346 because the condition on line 345 was never true
346 levels = None
347 else:
348 levels = np.linspace(vmin, vmax, 101)
349 norm = None
351 # Overwrite cmap, levels and norm for specific variables that
352 # require custom colorbar_map as these can not be defined in the
353 # JSON file.
354 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)
355 cmap, levels, norm = _custom_colourmap_visibility_in_air(
356 cube, cmap, levels, norm
357 )
358 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)
359 return cmap, levels, norm
362def _custom_colormap_scores(cube: iris.cube.Cube, cmap, levels, norm):
363 """Return altered colourmap for statistical metrics.
365 Parameters
366 ----------
367 cube: Cube
368 Cube of variable for which the colorbar information is desired.
369 cmap: Matplotlib colormap.
370 levels: List
371 List of levels to use for plotting. For continuous plots the min and max
372 should be taken as the range.
373 norm: BoundaryNorm.
375 Returns
376 -------
377 cmap: Matplotlib colormap.
378 levels: List
379 List of levels to use for plotting. For continuous plots the min and max
380 should be taken as the range.
381 norm: BoundaryNorm.
382 """
383 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
384 if any("RMSE_" in name for name in varnames): 384 ↛ 386line 384 didn't jump to line 386 because the condition on line 384 was always true
385 cmap = plt.get_cmap("PuRd", 51)
386 return cmap, levels, norm
389def _setup_spatial_map(
390 cube: iris.cube.Cube,
391 figure,
392 cmap,
393 grid_size: tuple[int, int] | None = None,
394 subplot: int | None = None,
395):
396 """Define map projections, extent and add coastlines and borderlines for spatial plots.
398 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
399 is specified, and map extent defined based on the input data.
401 Parameters
402 ----------
403 cube: Cube
404 2 dimensional (lat and lon) Cube of the data to plot.
405 figure:
406 Matplotlib Figure object holding all plot elements.
407 cmap:
408 Matplotlib colormap.
409 grid_size: (int, int), optional
410 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
411 subplot: int, optional
412 Subplot index if multiple spatial subplots in figure.
414 Returns
415 -------
416 axes:
417 Matplotlib GeoAxes definition.
418 """
419 # Identify min/max plot bounds.
420 try:
421 lat_axis, lon_axis = get_cube_yxcoordname(cube)
422 x1 = np.min(cube.coord(lon_axis).points)
423 x2 = np.max(cube.coord(lon_axis).points)
424 y1 = np.min(cube.coord(lat_axis).points)
425 y2 = np.max(cube.coord(lat_axis).points)
427 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
428 if np.abs(x2 - x1) > 180.0:
429 x1 = x1 - 180.0
430 x2 = x2 - 180.0
431 logging.debug("Adjusting plot bounds to fit global extent.")
433 # Consider map projection orientation.
434 # Adapting orientation enables plotting across international dateline.
435 # Users can adapt the default central_longitude if alternative projections views.
436 if x2 > 180.0 or x1 < -180.0:
437 central_longitude = 180.0
438 else:
439 central_longitude = 0.0
441 # Define spatial map projection.
442 coord_system = cube.coord(lat_axis).coord_system
443 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
444 # Define rotated pole map projection for rotated pole inputs.
445 projection = ccrs.RotatedPole(
446 pole_longitude=coord_system.grid_north_pole_longitude,
447 pole_latitude=coord_system.grid_north_pole_latitude,
448 central_rotated_longitude=central_longitude,
449 )
450 crs = projection
451 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 451 ↛ 453line 451 didn't jump to line 453 because the condition on line 451 was never true
452 # Define Transverse Mercator projection for TM inputs.
453 projection = ccrs.TransverseMercator(
454 central_longitude=coord_system.longitude_of_central_meridian,
455 central_latitude=coord_system.latitude_of_projection_origin,
456 false_easting=coord_system.false_easting,
457 false_northing=coord_system.false_northing,
458 scale_factor=coord_system.scale_factor_at_central_meridian,
459 )
460 crs = projection
461 else:
462 # Define regular map projection for non-rotated pole inputs.
463 # Alternatives might include e.g. for global model outputs:
464 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
465 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
466 projection = ccrs.PlateCarree(central_longitude=central_longitude)
467 crs = ccrs.PlateCarree()
469 # Define axes for plot (or subplot) with required map projection.
470 if subplot is not None:
471 axes = figure.add_subplot(
472 grid_size[0], grid_size[1], subplot, projection=projection
473 )
474 else:
475 axes = figure.add_subplot(projection=projection)
477 # Add coastlines and borderlines if cube contains x and y map coordinates.
478 # Avoid adding lines for masked data or specific fixed ancillary spatial plots.
479 if iris.util.is_masked(cube.data) or any( 479 ↛ 482line 479 didn't jump to line 482 because the condition on line 479 was never true
480 name in cube.name() for name in ["land_", "orography", "altitude"]
481 ):
482 pass
483 else:
484 if cmap.name in ["viridis", "Greys"]:
485 coastcol = "magenta"
486 else:
487 coastcol = "black"
488 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
489 axes.coastlines(resolution="10m", color=coastcol)
490 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
492 # Add gridlines.
493 gl = axes.gridlines(
494 alpha=0.3,
495 draw_labels=True,
496 dms=False,
497 x_inline=False,
498 y_inline=False,
499 )
500 gl.top_labels = False
501 gl.right_labels = False
502 if subplot:
503 gl.bottom_labels = False
504 gl.left_labels = False
505 if subplot % grid_size[1] == 1:
506 gl.left_labels = True
507 if subplot > ((grid_size[0] - 1) * grid_size[1]): 507 ↛ 512line 507 didn't jump to line 512 because the condition on line 507 was always true
508 gl.bottom_labels = True
510 # If is lat/lon spatial map, fix extent to keep plot tight.
511 # Specifying crs within set_extent helps ensure only data region is shown.
512 if isinstance(coord_system, iris.coord_systems.GeogCS):
513 axes.set_extent([x1, x2, y1, y2], crs=crs)
515 except ValueError:
516 # Skip if not both x and y map coordinates.
517 axes = figure.gca()
518 pass
520 return axes
523def _get_plot_resolution() -> int:
524 """Get resolution of rasterised plots in pixels per inch."""
525 return get_recipe_metadata().get("plot_resolution", 100)
528def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
529 """Return title and filename based on start and end points or bounds."""
530 if use_bounds and seq_coord.has_bounds():
531 vals = seq_coord.bounds.flatten()
532 else:
533 vals = seq_coord.points
534 start = seq_coord.units.title(vals[0])
535 end = seq_coord.units.title(vals[-1])
537 if start == end:
538 sequence_title = f"\n [{start}]"
539 sequence_fname = f"_{filename_slugify(start)}"
540 else:
541 sequence_title = f"\n [{start} to {end}]"
542 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
544 # Do not include time if coord set to zero.
545 if (
546 seq_coord.units == "hours since 0001-01-01 00:00:00"
547 and vals[0] == 0
548 and vals[-1] == 0
549 ):
550 sequence_title = ""
551 sequence_fname = ""
553 return sequence_title, sequence_fname
556def _set_title_and_filename(
557 seq_coord: iris.coords.Coord,
558 nplot: int,
559 recipe_title: str,
560 filename: str,
561):
562 """Set plot title and filename based on cube coordinate.
564 Parameters
565 ----------
566 sequence_coordinate: iris.coords.Coord
567 Coordinate about which to make a plot sequence.
568 nplot: int
569 Number of output plots to generate - controls title/naming.
570 recipe_title: str
571 Default plot title, potentially to update.
572 filename: str
573 Input plot filename, potentially to update.
575 Returns
576 -------
577 plot_title: str
578 Output formatted plot title string, based on plotted data.
579 plot_filename: str
580 Output formatted plot filename string.
581 """
582 ndim = seq_coord.ndim
583 npoints = np.size(seq_coord.points)
584 sequence_title = ""
585 sequence_fname = ""
587 # Case 1: Multiple dimension sequence input - list number of aggregated cases
588 # (e.g. aggregation histogram plots)
589 if ndim > 1:
590 ncase = np.shape(seq_coord)[0]
591 sequence_title = f"\n [{ncase} cases]"
592 sequence_fname = f"_{ncase}cases"
594 # Case 2: Single dimension input
595 else:
596 # Single sequence point
597 if npoints == 1:
598 if nplot > 1:
599 # Default labels for sequence inputs
600 sequence_value = seq_coord.units.title(seq_coord.points[0])
601 sequence_title = f"\n [{sequence_value}]"
602 sequence_fname = f"_{filename_slugify(sequence_value)}"
603 else:
604 # Aggregated attribute available where input collapsed over aggregation
605 try:
606 ncase = seq_coord.attributes["number_reference_times"]
607 sequence_title = f"\n [{ncase} cases]"
608 sequence_fname = f"_{ncase}cases"
609 except KeyError:
610 sequence_title, sequence_fname = _get_start_end_strings(
611 seq_coord, use_bounds=seq_coord.has_bounds()
612 )
613 # Multiple sequence (e.g. time) points
614 else:
615 sequence_title, sequence_fname = _get_start_end_strings(
616 seq_coord, use_bounds=False
617 )
619 # Set plot title and filename
620 plot_title = f"{recipe_title}{sequence_title}"
622 # Set plot filename, defaulting to user input if provided.
623 if filename is None:
624 filename = slugify(recipe_title)
625 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
626 else:
627 if nplot > 1:
628 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
629 else:
630 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
632 return plot_title, plot_filename
635def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
636 """Control postage stamp plot output titles based on stamp coordinate."""
637 if stamp_coord.name() == "realization":
638 mtitle = "Member"
639 else:
640 mtitle = stamp_coord.name().capitalize()
642 if stamp_coord.name() == "time":
643 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
644 else:
645 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
647 return mtitle
650def _plot_and_save_spatial_plot(
651 cube: iris.cube.Cube,
652 filename: str,
653 title: str,
654 method: Literal["contourf", "pcolormesh"],
655 overlay_cube: iris.cube.Cube | None = None,
656 contour_cube: iris.cube.Cube | None = None,
657 **kwargs,
658):
659 """Plot and save a spatial plot.
661 Parameters
662 ----------
663 cube: Cube
664 2 dimensional (lat and lon) Cube of the data to plot.
665 filename: str
666 Filename of the plot to write.
667 title: str
668 Plot title.
669 method: "contourf" | "pcolormesh"
670 The plotting method to use.
671 overlay_cube: Cube, optional
672 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
673 contour_cube: Cube, optional
674 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
675 """
676 # Setup plot details, size, resolution, etc.
677 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
679 # Specify the color bar
680 cmap, levels, norm = _colorbar_map_levels(cube)
682 # If overplotting, set required colorbars
683 if overlay_cube:
684 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
685 if contour_cube:
686 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
688 # Setup plot map projection, extent and coastlines and borderlines.
689 axes = _setup_spatial_map(cube, fig, cmap)
691 # Plot the field.
692 if method == "contourf":
693 # Filled contour plot of the field.
694 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
695 elif method == "pcolormesh":
696 try:
697 vmin = min(levels)
698 vmax = max(levels)
699 except TypeError:
700 vmin, vmax = None, None
701 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
702 # if levels are defined.
703 if norm is not None:
704 vmin = None
705 vmax = None
706 logging.debug("Plotting using defined levels.")
707 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
708 else:
709 raise ValueError(f"Unknown plotting method: {method}")
711 # Overplot overlay field, if required
712 if overlay_cube:
713 try:
714 over_vmin = min(over_levels)
715 over_vmax = max(over_levels)
716 except TypeError:
717 over_vmin, over_vmax = None, None
718 if over_norm is not None: 718 ↛ 719line 718 didn't jump to line 719 because the condition on line 718 was never true
719 over_vmin = None
720 over_vmax = None
721 overlay = iplt.pcolormesh(
722 overlay_cube,
723 cmap=over_cmap,
724 norm=over_norm,
725 alpha=0.8,
726 vmin=over_vmin,
727 vmax=over_vmax,
728 )
729 # Overplot contour field, if required, with contour labelling.
730 if contour_cube:
731 contour = iplt.contour(
732 contour_cube,
733 colors="darkgray",
734 levels=cntr_levels,
735 norm=cntr_norm,
736 alpha=0.5,
737 linestyles="--",
738 linewidths=1,
739 )
740 plt.clabel(contour)
742 # Check to see if transect, and if so, adjust y axis.
743 if is_transect(cube):
744 if "pressure" in [coord.name() for coord in cube.coords()]:
745 axes.invert_yaxis()
746 axes.set_yscale("log")
747 axes.set_ylim(1100, 100)
748 # If both model_level_number and level_height exists, iplt can construct
749 # plot as a function of height above orography (NOT sea level).
750 elif {"model_level_number", "level_height"}.issubset( 750 ↛ 755line 750 didn't jump to line 755 because the condition on line 750 was always true
751 {coord.name() for coord in cube.coords()}
752 ):
753 axes.set_yscale("log")
755 axes.set_title(
756 f"{title}\n"
757 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
758 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
759 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
760 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
761 fontsize=16,
762 )
764 # Inset code
765 axins = inset_axes(
766 axes,
767 width="20%",
768 height="20%",
769 loc="upper right",
770 axes_class=GeoAxes,
771 axes_kwargs={"map_projection": ccrs.PlateCarree()},
772 )
774 # Slightly transparent to reduce plot blocking.
775 axins.patch.set_alpha(0.4)
777 axins.coastlines(resolution="50m")
778 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
780 SLat, SLon, ELat, ELon = (
781 float(coord) for coord in cube.attributes["transect_coords"].split("_")
782 )
784 # Draw line between them
785 axins.plot(
786 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
787 )
789 # Plot points (note: lon, lat order for Cartopy)
790 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
791 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
793 lon_min, lon_max = sorted([SLon, ELon])
794 lat_min, lat_max = sorted([SLat, ELat])
796 # Midpoints
797 lon_mid = (lon_min + lon_max) / 2
798 lat_mid = (lat_min + lat_max) / 2
800 # Maximum half-range
801 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
802 if half_range == 0: # points identical → provide small default 802 ↛ 806line 802 didn't jump to line 806 because the condition on line 802 was always true
803 half_range = 1
805 # Set square extent
806 axins.set_extent(
807 [
808 lon_mid - half_range,
809 lon_mid + half_range,
810 lat_mid - half_range,
811 lat_mid + half_range,
812 ],
813 crs=ccrs.PlateCarree(),
814 )
816 # Ensure square aspect
817 axins.set_aspect("equal")
819 else:
820 # Add title.
821 axes.set_title(title, fontsize=16)
823 # Adjust padding if spatial plot or transect
824 if is_transect(cube):
825 yinfopad = -0.1
826 ycbarpad = 0.1
827 else:
828 yinfopad = 0.01
829 ycbarpad = 0.042
831 # Add watermark with min/max/mean. Currently not user togglable.
832 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
833 axes.annotate(
834 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
835 xy=(0.025, yinfopad),
836 xycoords="axes fraction",
837 xytext=(-5, 5),
838 textcoords="offset points",
839 ha="left",
840 va="bottom",
841 size=11,
842 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
843 )
845 # Add secondary colour bar for overlay_cube field if required.
846 if overlay_cube:
847 cbarB = fig.colorbar(
848 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
849 )
850 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
851 # add ticks and tick_labels for every levels if less than 20 levels exist
852 if over_levels is not None and len(over_levels) < 20: 852 ↛ 853line 852 didn't jump to line 853 because the condition on line 852 was never true
853 cbarB.set_ticks(over_levels)
854 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
855 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
856 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
857 logging.debug("Set secondary colorbar ticks and labels.")
859 # Add main colour bar.
860 cbar = fig.colorbar(
861 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
862 )
864 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
865 # add ticks and tick_labels for every levels if less than 20 levels exist
866 if levels is not None and len(levels) < 20:
867 cbar.set_ticks(levels)
868 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
869 if "rainfall" or "snowfall" or "visibility" in cube.name(): 869 ↛ 871line 869 didn't jump to line 871 because the condition on line 869 was always true
870 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
871 logging.debug("Set colorbar ticks and labels.")
873 # Save plot.
874 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
875 logging.info("Saved spatial plot to %s", filename)
876 plt.close(fig)
879def _plot_and_save_postage_stamp_spatial_plot(
880 cube: iris.cube.Cube,
881 filename: str,
882 stamp_coordinate: str,
883 title: str,
884 method: Literal["contourf", "pcolormesh"],
885 overlay_cube: iris.cube.Cube | None = None,
886 contour_cube: iris.cube.Cube | None = None,
887 **kwargs,
888):
889 """Plot postage stamp spatial plots from an ensemble.
891 Parameters
892 ----------
893 cube: Cube
894 Iris cube of data to be plotted. It must have the stamp coordinate.
895 filename: str
896 Filename of the plot to write.
897 stamp_coordinate: str
898 Coordinate that becomes different plots.
899 method: "contourf" | "pcolormesh"
900 The plotting method to use.
901 overlay_cube: Cube, optional
902 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
903 contour_cube: Cube, optional
904 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
906 Raises
907 ------
908 ValueError
909 If the cube doesn't have the right dimensions.
910 """
911 # Use the smallest square grid that will fit the members.
912 nmember = len(cube.coord(stamp_coordinate).points)
913 grid_rows = int(math.sqrt(nmember))
914 grid_size = math.ceil(nmember / grid_rows)
916 fig = plt.figure(
917 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
918 )
920 # Specify the color bar
921 cmap, levels, norm = _colorbar_map_levels(cube)
922 # If overplotting, set required colorbars
923 if overlay_cube: 923 ↛ 924line 923 didn't jump to line 924 because the condition on line 923 was never true
924 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
925 if contour_cube: 925 ↛ 926line 925 didn't jump to line 926 because the condition on line 925 was never true
926 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
928 # Make a subplot for each member.
929 for member, subplot in zip(
930 cube.slices_over(stamp_coordinate),
931 range(1, grid_size * grid_rows + 1),
932 strict=False,
933 ):
934 # Setup subplot map projection, extent and coastlines and borderlines.
935 axes = _setup_spatial_map(
936 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
937 )
938 if method == "contourf":
939 # Filled contour plot of the field.
940 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
941 elif method == "pcolormesh":
942 if levels is not None:
943 vmin = min(levels)
944 vmax = max(levels)
945 else:
946 raise TypeError("Unknown vmin and vmax range.")
947 vmin, vmax = None, None
948 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
949 # if levels are defined.
950 if norm is not None: 950 ↛ 951line 950 didn't jump to line 951 because the condition on line 950 was never true
951 vmin = None
952 vmax = None
953 # pcolormesh plot of the field.
954 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
955 else:
956 raise ValueError(f"Unknown plotting method: {method}")
958 # Overplot overlay field, if required
959 if overlay_cube: 959 ↛ 960line 959 didn't jump to line 960 because the condition on line 959 was never true
960 try:
961 over_vmin = min(over_levels)
962 over_vmax = max(over_levels)
963 except TypeError:
964 over_vmin, over_vmax = None, None
965 if over_norm is not None:
966 over_vmin = None
967 over_vmax = None
968 iplt.pcolormesh(
969 overlay_cube[member.coord(stamp_coordinate).points[0]],
970 cmap=over_cmap,
971 norm=over_norm,
972 alpha=0.6,
973 vmin=over_vmin,
974 vmax=over_vmax,
975 )
976 # Overplot contour field, if required
977 if contour_cube: 977 ↛ 978line 977 didn't jump to line 978 because the condition on line 977 was never true
978 iplt.contour(
979 contour_cube[member.coord(stamp_coordinate).points[0]],
980 colors="darkgray",
981 levels=cntr_levels,
982 norm=cntr_norm,
983 alpha=0.6,
984 linestyles="--",
985 linewidths=1,
986 )
987 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
988 axes.set_title(f"{mtitle}")
990 # Put the shared colorbar in its own axes.
991 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
992 colorbar = fig.colorbar(
993 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
994 )
995 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
997 # Overall figure title.
998 fig.suptitle(title, fontsize=16)
1000 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1001 logging.info("Saved contour postage stamp plot to %s", filename)
1002 plt.close(fig)
1005def _plot_and_save_line_series(
1006 cubes: iris.cube.CubeList,
1007 coords: list[iris.coords.Coord],
1008 ensemble_coord: str,
1009 filename: str,
1010 title: str,
1011 **kwargs,
1012):
1013 """Plot and save a 1D line series.
1015 Parameters
1016 ----------
1017 cubes: Cube or CubeList
1018 Cube or CubeList containing the cubes to plot on the y-axis.
1019 coords: list[Coord]
1020 Coordinates to plot on the x-axis, one per cube.
1021 ensemble_coord: str
1022 Ensemble coordinate in the cube.
1023 filename: str
1024 Filename of the plot to write.
1025 title: str
1026 Plot title.
1027 """
1028 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1030 model_colors_map = _get_model_colors_map(cubes)
1032 # Store min/max ranges.
1033 y_levels = []
1035 # Check match-up across sequence coords gives consistent sizes
1036 _validate_cubes_coords(cubes, coords)
1038 for cube, coord in zip(cubes, coords, strict=True):
1039 label = None
1040 color = "black"
1041 if model_colors_map:
1042 label = cube.attributes.get("model_name")
1043 color = model_colors_map.get(label)
1044 for cube_slice in cube.slices_over(ensemble_coord):
1045 # Label with (control) if part of an ensemble or not otherwise.
1046 if cube_slice.coord(ensemble_coord).points == [0]:
1047 iplt.plot(
1048 coord,
1049 cube_slice,
1050 color=color,
1051 marker="o",
1052 ls="-",
1053 lw=3,
1054 label=f"{label} (control)"
1055 if len(cube.coord(ensemble_coord).points) > 1
1056 else label,
1057 )
1058 # Label with (perturbed) if part of an ensemble and not the control.
1059 else:
1060 iplt.plot(
1061 coord,
1062 cube_slice,
1063 color=color,
1064 ls="-",
1065 lw=1.5,
1066 alpha=0.75,
1067 label=f"{label} (member)",
1068 )
1070 # Calculate the global min/max if multiple cubes are given.
1071 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1072 if levels is not None: 1072 ↛ 1073line 1072 didn't jump to line 1073 because the condition on line 1072 was never true
1073 y_levels.append(min(levels))
1074 y_levels.append(max(levels))
1076 # Get the current axes.
1077 ax = plt.gca()
1079 # Add some labels and tweak the style.
1080 # check if cubes[0] works for single cube if not CubeList
1081 if coords[0].name() == "time":
1082 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1083 else:
1084 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1085 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1086 ax.set_title(title, fontsize=16)
1088 ax.ticklabel_format(axis="y", useOffset=False)
1089 ax.tick_params(axis="x", labelrotation=15)
1090 ax.tick_params(axis="both", labelsize=12)
1092 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1093 if y_levels: 1093 ↛ 1094line 1093 didn't jump to line 1094 because the condition on line 1093 was never true
1094 ax.set_ylim(min(y_levels), max(y_levels))
1095 # Add zero line.
1096 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1097 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1098 logging.debug(
1099 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1100 )
1101 else:
1102 ax.autoscale()
1104 # Add gridlines
1105 ax.grid(linestyle="--", color="grey", linewidth=1)
1106 # Ientify unique labels for legend
1107 handles = list(
1108 {
1109 label: handle
1110 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1111 }.values()
1112 )
1113 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1115 # Save plot.
1116 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1117 logging.info("Saved line plot to %s", filename)
1118 plt.close(fig)
1121def _plot_and_save_vertical_line_series(
1122 cubes: iris.cube.CubeList,
1123 coords: list[iris.coords.Coord],
1124 ensemble_coord: str,
1125 filename: str,
1126 series_coordinate: str,
1127 title: str,
1128 vmin: float,
1129 vmax: float,
1130 **kwargs,
1131):
1132 """Plot and save a 1D line series in vertical.
1134 Parameters
1135 ----------
1136 cubes: CubeList
1137 1 dimensional Cube or CubeList of the data to plot on x-axis.
1138 coord: list[Coord]
1139 Coordinates to plot on the y-axis, one per cube.
1140 ensemble_coord: str
1141 Ensemble coordinate in the cube.
1142 filename: str
1143 Filename of the plot to write.
1144 series_coordinate: str
1145 Coordinate to use as vertical axis.
1146 title: str
1147 Plot title.
1148 vmin: float
1149 Minimum value for the x-axis.
1150 vmax: float
1151 Maximum value for the x-axis.
1152 """
1153 # plot the vertical pressure axis using log scale
1154 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1156 model_colors_map = _get_model_colors_map(cubes)
1158 # Check match-up across sequence coords gives consistent sizes
1159 _validate_cubes_coords(cubes, coords)
1161 for cube, coord in zip(cubes, coords, strict=True):
1162 label = None
1163 color = "black"
1164 if model_colors_map: 1164 ↛ 1165line 1164 didn't jump to line 1165 because the condition on line 1164 was never true
1165 label = cube.attributes.get("model_name")
1166 color = model_colors_map.get(label)
1168 for cube_slice in cube.slices_over(ensemble_coord):
1169 # If ensemble data given plot control member with (control)
1170 # unless single forecast.
1171 if cube_slice.coord(ensemble_coord).points == [0]:
1172 iplt.plot(
1173 cube_slice,
1174 coord,
1175 color=color,
1176 marker="o",
1177 ls="-",
1178 lw=3,
1179 label=f"{label} (control)"
1180 if len(cube.coord(ensemble_coord).points) > 1
1181 else label,
1182 )
1183 # If ensemble data given plot perturbed members with (perturbed).
1184 else:
1185 iplt.plot(
1186 cube_slice,
1187 coord,
1188 color=color,
1189 ls="-",
1190 lw=1.5,
1191 alpha=0.75,
1192 label=f"{label} (member)",
1193 )
1195 # Get the current axis
1196 ax = plt.gca()
1198 # Special handling for pressure level data.
1199 if series_coordinate == "pressure": 1199 ↛ 1221line 1199 didn't jump to line 1221 because the condition on line 1199 was always true
1200 # Invert y-axis and set to log scale.
1201 ax.invert_yaxis()
1202 ax.set_yscale("log")
1204 # Define y-ticks and labels for pressure log axis.
1205 y_tick_labels = [
1206 "1000",
1207 "850",
1208 "700",
1209 "500",
1210 "300",
1211 "200",
1212 "100",
1213 ]
1214 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1216 # Set y-axis limits and ticks.
1217 ax.set_ylim(1100, 100)
1219 # Test if series_coordinate is model level data. The UM data uses
1220 # model_level_number and lfric uses full_levels as coordinate.
1221 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1222 # Define y-ticks and labels for vertical axis.
1223 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1224 y_tick_labels = [str(int(i)) for i in y_ticks]
1225 ax.set_ylim(min(y_ticks), max(y_ticks))
1227 ax.set_yticks(y_ticks)
1228 ax.set_yticklabels(y_tick_labels)
1230 # Set x-axis limits.
1231 ax.set_xlim(vmin, vmax)
1232 # Mark y=0 if present in plot.
1233 if vmin < 0.0 and vmax > 0.0: 1233 ↛ 1234line 1233 didn't jump to line 1234 because the condition on line 1233 was never true
1234 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1236 # Add some labels and tweak the style.
1237 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1238 ax.set_xlabel(
1239 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1240 )
1241 ax.set_title(title, fontsize=16)
1242 ax.ticklabel_format(axis="x")
1243 ax.tick_params(axis="y")
1244 ax.tick_params(axis="both", labelsize=12)
1246 # Add gridlines
1247 ax.grid(linestyle="--", color="grey", linewidth=1)
1248 # Ientify unique labels for legend
1249 handles = list(
1250 {
1251 label: handle
1252 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1253 }.values()
1254 )
1255 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1257 # Save plot.
1258 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1259 logging.info("Saved line plot to %s", filename)
1260 plt.close(fig)
1263def _plot_and_save_scatter_plot(
1264 cube_x: iris.cube.Cube | iris.cube.CubeList,
1265 cube_y: iris.cube.Cube | iris.cube.CubeList,
1266 filename: str,
1267 title: str,
1268 one_to_one: bool,
1269 model_names: list[str] = None,
1270 **kwargs,
1271):
1272 """Plot and save a 2D scatter plot.
1274 Parameters
1275 ----------
1276 cube_x: Cube | CubeList
1277 1 dimensional Cube or CubeList of the data to plot on x-axis.
1278 cube_y: Cube | CubeList
1279 1 dimensional Cube or CubeList of the data to plot on y-axis.
1280 filename: str
1281 Filename of the plot to write.
1282 title: str
1283 Plot title.
1284 one_to_one: bool
1285 Whether a 1:1 line is plotted.
1286 """
1287 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1288 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1289 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1290 # over the pairs simultaneously.
1292 # Ensure cube_x and cube_y are iterable
1293 cube_x_iterable = iter_maybe(cube_x)
1294 cube_y_iterable = iter_maybe(cube_y)
1296 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1297 iplt.scatter(cube_x_iter, cube_y_iter)
1298 if one_to_one is True:
1299 plt.plot(
1300 [
1301 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1302 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1303 ],
1304 [
1305 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1306 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1307 ],
1308 "k",
1309 linestyle="--",
1310 )
1311 ax = plt.gca()
1313 # Add some labels and tweak the style.
1314 if model_names is None:
1315 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1316 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1317 else:
1318 # Add the model names, these should be order of base (x) and other (y).
1319 ax.set_xlabel(
1320 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1321 )
1322 ax.set_ylabel(
1323 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1324 )
1325 ax.set_title(title, fontsize=16)
1326 ax.ticklabel_format(axis="y", useOffset=False)
1327 ax.tick_params(axis="x", labelrotation=15)
1328 ax.tick_params(axis="both", labelsize=12)
1329 ax.autoscale()
1331 # Save plot.
1332 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1333 logging.info("Saved scatter plot to %s", filename)
1334 plt.close(fig)
1337def _plot_and_save_vector_plot(
1338 cube_u: iris.cube.Cube,
1339 cube_v: iris.cube.Cube,
1340 filename: str,
1341 title: str,
1342 method: Literal["contourf", "pcolormesh"],
1343 **kwargs,
1344):
1345 """Plot and save a 2D vector plot.
1347 Parameters
1348 ----------
1349 cube_u: Cube
1350 2 dimensional Cube of u component of the data.
1351 cube_v: Cube
1352 2 dimensional Cube of v component of the data.
1353 filename: str
1354 Filename of the plot to write.
1355 title: str
1356 Plot title.
1357 """
1358 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1360 # Create a cube containing the magnitude of the vector field.
1361 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1362 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1364 # Specify the color bar
1365 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1367 # Setup plot map projection, extent and coastlines and borderlines.
1368 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1370 if method == "contourf":
1371 # Filled contour plot of the field.
1372 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1373 elif method == "pcolormesh":
1374 try:
1375 vmin = min(levels)
1376 vmax = max(levels)
1377 except TypeError:
1378 vmin, vmax = None, None
1379 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1380 # if levels are defined.
1381 if norm is not None:
1382 vmin = None
1383 vmax = None
1384 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1385 else:
1386 raise ValueError(f"Unknown plotting method: {method}")
1388 # Check to see if transect, and if so, adjust y axis.
1389 if is_transect(cube_vec_mag):
1390 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1391 axes.invert_yaxis()
1392 axes.set_yscale("log")
1393 axes.set_ylim(1100, 100)
1394 # If both model_level_number and level_height exists, iplt can construct
1395 # plot as a function of height above orography (NOT sea level).
1396 elif {"model_level_number", "level_height"}.issubset(
1397 {coord.name() for coord in cube_vec_mag.coords()}
1398 ):
1399 axes.set_yscale("log")
1401 axes.set_title(
1402 f"{title}\n"
1403 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1404 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1405 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1406 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1407 fontsize=16,
1408 )
1410 else:
1411 # Add title.
1412 axes.set_title(title, fontsize=16)
1414 # Add watermark with min/max/mean. Currently not user togglable.
1415 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1416 axes.annotate(
1417 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}",
1418 xy=(0.05, -0.05),
1419 xycoords="axes fraction",
1420 xytext=(-5, 5),
1421 textcoords="offset points",
1422 ha="right",
1423 va="bottom",
1424 size=11,
1425 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1426 )
1428 # Add colour bar.
1429 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1430 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1431 # add ticks and tick_labels for every levels if less than 20 levels exist
1432 if levels is not None and len(levels) < 20:
1433 cbar.set_ticks(levels)
1434 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1436 # 30 barbs along the longest axis of the plot, or a barb per point for data
1437 # with less than 30 points.
1438 step = max(max(cube_u.shape) // 30, 1)
1439 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1441 # Save plot.
1442 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1443 logging.info("Saved vector plot to %s", filename)
1444 plt.close(fig)
1447def _plot_and_save_histogram_series(
1448 cubes: iris.cube.Cube | iris.cube.CubeList,
1449 filename: str,
1450 title: str,
1451 vmin: float,
1452 vmax: float,
1453 **kwargs,
1454):
1455 """Plot and save a histogram series.
1457 Parameters
1458 ----------
1459 cubes: Cube or CubeList
1460 2 dimensional Cube or CubeList of the data to plot as histogram.
1461 filename: str
1462 Filename of the plot to write.
1463 title: str
1464 Plot title.
1465 vmin: float
1466 minimum for colorbar
1467 vmax: float
1468 maximum for colorbar
1469 """
1470 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1471 ax = plt.gca()
1473 model_colors_map = _get_model_colors_map(cubes)
1475 # Set default that histograms will produce probability density function
1476 # at each bin (integral over range sums to 1).
1477 density = True
1479 for cube in iter_maybe(cubes):
1480 # Easier to check title (where var name originates)
1481 # than seeing if long names exist etc.
1482 # Exception case, where distribution better fits log scales/bins.
1483 if "surface_microphysical" in title:
1484 if "amount" in title: 1484 ↛ 1486line 1484 didn't jump to line 1486 because the condition on line 1484 was never true
1485 # Compute histogram following Klingaman et al. (2017): ASoP
1486 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1487 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1488 density = False
1489 else:
1490 bins = 10.0 ** (
1491 np.arange(-10, 27, 1) / 10.0
1492 ) # Suggestion from RMED toolbox.
1493 bins = np.insert(bins, 0, 0)
1494 ax.set_yscale("log")
1495 vmin = bins[1]
1496 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1497 ax.set_xscale("log")
1498 elif "lightning" in title:
1499 bins = [0, 1, 2, 3, 4, 5]
1500 else:
1501 bins = np.linspace(vmin, vmax, 51)
1502 logging.debug(
1503 "Plotting histogram with %s bins %s - %s.",
1504 np.size(bins),
1505 np.min(bins),
1506 np.max(bins),
1507 )
1509 # Reshape cube data into a single array to allow for a single histogram.
1510 # Otherwise we plot xdim histograms stacked.
1511 cube_data_1d = (cube.data).flatten()
1513 label = None
1514 color = "black"
1515 if model_colors_map: 1515 ↛ 1516line 1515 didn't jump to line 1516 because the condition on line 1515 was never true
1516 label = cube.attributes.get("model_name")
1517 color = model_colors_map[label]
1518 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1520 # Compute area under curve.
1521 if "surface_microphysical" in title and "amount" in title: 1521 ↛ 1522line 1521 didn't jump to line 1522 because the condition on line 1521 was never true
1522 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1523 x = x * bin_mean / x.sum()
1524 x = x[1:]
1525 y = y[1:]
1527 ax.plot(
1528 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1529 )
1531 # Add some labels and tweak the style.
1532 ax.set_title(title, fontsize=16)
1533 ax.set_xlabel(
1534 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1535 )
1536 ax.set_ylabel("Normalised probability density", fontsize=14)
1537 if "surface_microphysical" in title and "amount" in title: 1537 ↛ 1538line 1537 didn't jump to line 1538 because the condition on line 1537 was never true
1538 ax.set_ylabel(
1539 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1540 )
1541 ax.set_xlim(vmin, vmax)
1542 ax.tick_params(axis="both", labelsize=12)
1544 # Overlay grid-lines onto histogram plot.
1545 ax.grid(linestyle="--", color="grey", linewidth=1)
1546 if model_colors_map: 1546 ↛ 1547line 1546 didn't jump to line 1547 because the condition on line 1546 was never true
1547 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1549 # Save plot.
1550 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1551 logging.info("Saved histogram plot to %s", filename)
1552 plt.close(fig)
1555def _plot_and_save_postage_stamp_histogram_series(
1556 cube: iris.cube.Cube,
1557 filename: str,
1558 title: str,
1559 stamp_coordinate: str,
1560 vmin: float,
1561 vmax: float,
1562 **kwargs,
1563):
1564 """Plot and save postage (ensemble members) stamps for a histogram series.
1566 Parameters
1567 ----------
1568 cube: Cube
1569 2 dimensional Cube of the data to plot as histogram.
1570 filename: str
1571 Filename of the plot to write.
1572 title: str
1573 Plot title.
1574 stamp_coordinate: str
1575 Coordinate that becomes different plots.
1576 vmin: float
1577 minimum for pdf x-axis
1578 vmax: float
1579 maximum for pdf x-axis
1580 """
1581 # Use the smallest square grid that will fit the members.
1582 nmember = len(cube.coord(stamp_coordinate).points)
1583 grid_rows = int(math.sqrt(nmember))
1584 grid_size = math.ceil(nmember / grid_rows)
1586 fig = plt.figure(
1587 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1588 )
1589 # Make a subplot for each member.
1590 for member, subplot in zip(
1591 cube.slices_over(stamp_coordinate),
1592 range(1, grid_size * grid_rows + 1),
1593 strict=False,
1594 ):
1595 # Implicit interface is much easier here, due to needing to have the
1596 # cartopy GeoAxes generated.
1597 plt.subplot(grid_rows, grid_size, subplot)
1598 # Reshape cube data into a single array to allow for a single histogram.
1599 # Otherwise we plot xdim histograms stacked.
1600 member_data_1d = (member.data).flatten()
1601 plt.hist(member_data_1d, density=True, stacked=True)
1602 axes = plt.gca()
1603 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1604 axes.set_title(f"{mtitle}")
1605 axes.set_xlim(vmin, vmax)
1607 # Overall figure title.
1608 fig.suptitle(title, fontsize=16)
1610 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1611 logging.info("Saved histogram postage stamp plot to %s", filename)
1612 plt.close(fig)
1615def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1616 cube: iris.cube.Cube,
1617 filename: str,
1618 title: str,
1619 stamp_coordinate: str,
1620 vmin: float,
1621 vmax: float,
1622 **kwargs,
1623):
1624 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1625 ax.set_title(title, fontsize=16)
1626 ax.set_xlim(vmin, vmax)
1627 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1628 ax.set_ylabel("normalised probability density", fontsize=14)
1629 # Loop over all slices along the stamp_coordinate
1630 for member in cube.slices_over(stamp_coordinate):
1631 # Flatten the member data to 1D
1632 member_data_1d = member.data.flatten()
1633 # Plot the histogram using plt.hist
1634 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1635 plt.hist(
1636 member_data_1d,
1637 density=True,
1638 stacked=True,
1639 label=f"{mtitle}",
1640 )
1642 # Add a legend
1643 ax.legend(fontsize=16)
1645 # Save the figure to a file
1646 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1647 logging.info("Saved histogram postage stamp plot to %s", filename)
1649 # Close the figure
1650 plt.close(fig)
1653def _plot_and_save_scattermap_plot(
1654 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1655):
1656 """Plot and save a geographical scatter plot.
1658 Parameters
1659 ----------
1660 cube: Cube
1661 1 dimensional Cube of the data points with auxiliary latitude and
1662 longitude coordinates,
1663 filename: str
1664 Filename of the plot to write.
1665 title: str
1666 Plot title.
1667 projection: str
1668 Mapping projection to be used by cartopy.
1669 """
1670 # Setup plot details, size, resolution, etc.
1671 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1672 if projection is not None:
1673 # Apart from the default, the only projection we currently support is
1674 # a stereographic projection over the North Pole.
1675 if projection == "NP_Stereo":
1676 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1677 else:
1678 raise ValueError(f"Unknown projection: {projection}")
1679 else:
1680 axes = plt.axes(projection=ccrs.PlateCarree())
1682 # Scatter plot of the field. The marker size is chosen to give
1683 # symbols that decrease in size as the number of observations
1684 # increases, although the fraction of the figure covered by
1685 # symbols increases roughly as N^(1/2), disregarding overlaps,
1686 # and has been selected for the default figure size of (10, 10).
1687 # Should this be changed, the marker size should be adjusted in
1688 # proportion to the area of the figure.
1689 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1690 klon = None
1691 klat = None
1692 for kc in range(len(cube.aux_coords)):
1693 if cube.aux_coords[kc].standard_name == "latitude":
1694 klat = kc
1695 elif cube.aux_coords[kc].standard_name == "longitude":
1696 klon = kc
1697 scatter_map = iplt.scatter(
1698 cube.aux_coords[klon],
1699 cube.aux_coords[klat],
1700 c=cube.data[:],
1701 s=mrk_size,
1702 cmap="jet",
1703 edgecolors="k",
1704 )
1706 # Add coastlines and borderlines.
1707 try:
1708 axes.coastlines(resolution="10m")
1709 axes.add_feature(cfeature.BORDERS)
1710 except AttributeError:
1711 pass
1713 # Add title.
1714 axes.set_title(title, fontsize=16)
1716 # Add colour bar.
1717 cbar = fig.colorbar(scatter_map)
1718 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1720 # Save plot.
1721 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1722 logging.info("Saved geographical scatter plot to %s", filename)
1723 plt.close(fig)
1726def _plot_and_save_power_spectrum_series(
1727 cubes: iris.cube.Cube | iris.cube.CubeList,
1728 filename: str,
1729 title: str,
1730 **kwargs,
1731):
1732 """Plot and save a power spectrum series.
1734 Parameters
1735 ----------
1736 cubes: Cube or CubeList
1737 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1738 filename: str
1739 Filename of the plot to write.
1740 title: str
1741 Plot title.
1742 """
1743 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1744 ax = plt.gca()
1746 model_colors_map = _get_model_colors_map(cubes)
1748 for cube in iter_maybe(cubes):
1749 # Calculate power spectrum
1751 # Extract time coordinate and convert to datetime
1752 time_coord = cube.coord("time")
1753 time_points = time_coord.units.num2date(time_coord.points)
1755 # Choose one time point (e.g., the first one)
1756 target_time = time_points[0]
1758 # Bind target_time inside the lambda using a default argument
1759 time_constraint = iris.Constraint(
1760 time=lambda cell, target_time=target_time: cell.point == target_time
1761 )
1763 cube = cube.extract(time_constraint)
1765 if cube.ndim == 2:
1766 cube_3d = cube.data[np.newaxis, :, :]
1767 logging.debug("Adding in new axis for a 2 dimensional cube.")
1768 elif cube.ndim == 3: 1768 ↛ 1769line 1768 didn't jump to line 1769 because the condition on line 1768 was never true
1769 cube_3d = cube.data
1770 else:
1771 raise ValueError("Cube dimensions unsuitable for power spectra code")
1772 raise ValueError(
1773 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1774 )
1776 # Calculate spectra
1777 ps_array = _DCT_ps(cube_3d)
1779 ps_cube = iris.cube.Cube(
1780 ps_array,
1781 long_name="power_spectra",
1782 )
1784 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1786 # Create a frequency/wavelength array for coordinate
1787 ps_len = ps_cube.data.shape[1]
1788 freqs = np.arange(1, ps_len + 1)
1789 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1791 # Convert datetime to numeric time using original units
1792 numeric_time = time_coord.units.date2num(time_points)
1793 # Create a new DimCoord with numeric time
1794 new_time_coord = iris.coords.DimCoord(
1795 numeric_time, standard_name="time", units=time_coord.units
1796 )
1798 # Add time and frequency coordinate to spectra cube.
1799 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1800 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1802 # Extract data from the cube
1803 frequency = ps_cube.coord("frequency").points
1804 power_spectrum = ps_cube.data
1806 label = None
1807 color = "black"
1808 if model_colors_map: 1808 ↛ 1809line 1808 didn't jump to line 1809 because the condition on line 1808 was never true
1809 label = ps_cube.attributes.get("model_name")
1810 color = model_colors_map[label]
1811 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1813 # Add some labels and tweak the style.
1814 ax.set_title(title, fontsize=16)
1815 ax.set_xlabel("Wavenumber", fontsize=14)
1816 ax.set_ylabel("Power", fontsize=14)
1817 ax.tick_params(axis="both", labelsize=12)
1819 # Set log-log scale
1820 ax.set_xscale("log")
1821 ax.set_yscale("log")
1823 # Overlay grid-lines onto power spectrum plot.
1824 ax.grid(linestyle="--", color="grey", linewidth=1)
1825 if model_colors_map: 1825 ↛ 1826line 1825 didn't jump to line 1826 because the condition on line 1825 was never true
1826 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1828 # Save plot.
1829 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1830 logging.info("Saved power spectrum plot to %s", filename)
1831 plt.close(fig)
1834def _plot_and_save_postage_stamp_power_spectrum_series(
1835 cube: iris.cube.Cube,
1836 filename: str,
1837 title: str,
1838 stamp_coordinate: str,
1839 **kwargs,
1840):
1841 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1843 Parameters
1844 ----------
1845 cube: Cube
1846 2 dimensional Cube of the data to plot as power spectrum.
1847 filename: str
1848 Filename of the plot to write.
1849 title: str
1850 Plot title.
1851 stamp_coordinate: str
1852 Coordinate that becomes different plots.
1853 """
1854 # Use the smallest square grid that will fit the members.
1855 nmember = len(cube.coord(stamp_coordinate).points)
1856 grid_rows = int(math.sqrt(nmember))
1857 grid_size = math.ceil(nmember / grid_rows)
1859 fig = plt.figure(
1860 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1861 )
1863 # Make a subplot for each member.
1864 for member, subplot in zip(
1865 cube.slices_over(stamp_coordinate),
1866 range(1, grid_size * grid_rows + 1),
1867 strict=False,
1868 ):
1869 # Implicit interface is much easier here, due to needing to have the
1870 # cartopy GeoAxes generated.
1871 plt.subplot(grid_rows, grid_size, subplot)
1873 frequency = member.coord("frequency").points
1875 axes = plt.gca()
1876 axes.plot(frequency, member.data)
1877 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1879 # Overall figure title.
1880 fig.suptitle(title, fontsize=16)
1882 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1883 logging.info("Saved power spectra postage stamp plot to %s", filename)
1884 plt.close(fig)
1887def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1888 cube: iris.cube.Cube,
1889 filename: str,
1890 title: str,
1891 stamp_coordinate: str,
1892 **kwargs,
1893):
1894 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1895 ax.set_title(title, fontsize=16)
1896 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1897 ax.set_ylabel("Power", fontsize=14)
1898 # Loop over all slices along the stamp_coordinate
1899 for member in cube.slices_over(stamp_coordinate):
1900 frequency = member.coord("frequency").points
1901 ax.plot(
1902 frequency,
1903 member.data,
1904 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1905 )
1907 # Add a legend
1908 ax.legend(fontsize=16)
1910 # Save the figure to a file
1911 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1912 logging.info("Saved power spectra plot to %s", filename)
1914 # Close the figure
1915 plt.close(fig)
1918def _spatial_plot(
1919 method: Literal["contourf", "pcolormesh"],
1920 cube: iris.cube.Cube,
1921 filename: str | None,
1922 sequence_coordinate: str,
1923 stamp_coordinate: str,
1924 overlay_cube: iris.cube.Cube | None = None,
1925 contour_cube: iris.cube.Cube | None = None,
1926 **kwargs,
1927):
1928 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1930 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1931 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1932 is present then postage stamp plots will be produced.
1934 If an overlay_cube and/or contour_cube are specified, multiple variables can
1935 be overplotted on the same figure.
1937 Parameters
1938 ----------
1939 method: "contourf" | "pcolormesh"
1940 The plotting method to use.
1941 cube: Cube
1942 Iris cube of the data to plot. It should have two spatial dimensions,
1943 such as lat and lon, and may also have a another two dimension to be
1944 plotted sequentially and/or as postage stamp plots.
1945 filename: str | None
1946 Name of the plot to write, used as a prefix for plot sequences. If None
1947 uses the recipe name.
1948 sequence_coordinate: str
1949 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1950 This coordinate must exist in the cube.
1951 stamp_coordinate: str
1952 Coordinate about which to plot postage stamp plots. Defaults to
1953 ``"realization"``.
1954 overlay_cube: Cube | None, optional
1955 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1956 contour_cube: Cube | None, optional
1957 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1959 Raises
1960 ------
1961 ValueError
1962 If the cube doesn't have the right dimensions.
1963 TypeError
1964 If the cube isn't a single cube.
1965 """
1966 recipe_title = get_recipe_metadata().get("title", "Untitled")
1968 # Ensure we've got a single cube.
1969 cube = _check_single_cube(cube)
1971 # Check if there is a valid stamp coordinate in cube dimensions.
1972 if stamp_coordinate == "realization": 1972 ↛ 1977line 1972 didn't jump to line 1977 because the condition on line 1972 was always true
1973 stamp_coordinate = check_stamp_coordinate(cube)
1975 # Make postage stamp plots if stamp_coordinate exists and has more than a
1976 # single point.
1977 plotting_func = _plot_and_save_spatial_plot
1978 try:
1979 if cube.coord(stamp_coordinate).shape[0] > 1:
1980 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1981 except iris.exceptions.CoordinateNotFoundError:
1982 pass
1984 # Produce a geographical scatter plot if the data have a
1985 # dimension called observation or model_obs_error
1986 if any( 1986 ↛ 1990line 1986 didn't jump to line 1990 because the condition on line 1986 was never true
1987 crd.var_name == "station" or crd.var_name == "model_obs_error"
1988 for crd in cube.coords()
1989 ):
1990 plotting_func = _plot_and_save_scattermap_plot
1992 # Must have a sequence coordinate.
1993 try:
1994 cube.coord(sequence_coordinate)
1995 except iris.exceptions.CoordinateNotFoundError as err:
1996 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1998 # Create a plot for each value of the sequence coordinate.
1999 plot_index = []
2000 nplot = np.size(cube.coord(sequence_coordinate).points)
2002 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
2003 # Set plot titles and filename
2004 seq_coord = cube_slice.coord(sequence_coordinate)
2005 plot_title, plot_filename = _set_title_and_filename(
2006 seq_coord, nplot, recipe_title, filename
2007 )
2009 # Extract sequence slice for overlay_cube and contour_cube if required.
2010 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
2011 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
2013 # Do the actual plotting.
2014 plotting_func(
2015 cube_slice,
2016 filename=plot_filename,
2017 stamp_coordinate=stamp_coordinate,
2018 title=plot_title,
2019 method=method,
2020 overlay_cube=overlay_slice,
2021 contour_cube=contour_slice,
2022 **kwargs,
2023 )
2024 plot_index.append(plot_filename)
2026 # Add list of plots to plot metadata.
2027 complete_plot_index = _append_to_plot_index(plot_index)
2029 # Make a page to display the plots.
2030 _make_plot_html_page(complete_plot_index)
2033def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2034 """Get colourmap for mask.
2036 If "mask_for_" appears anywhere in the name of a cube this function will be called
2037 regardless of the name of the variable to ensure a consistent plot.
2039 Parameters
2040 ----------
2041 cube: Cube
2042 Cube of variable for which the colorbar information is desired.
2043 axis: "x", "y", optional
2044 Select the levels for just this axis of a line plot. The min and max
2045 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2046 setting a universal range is not desirable (e.g. temperature), users
2047 can set ymin/ymax values to "auto" in the colorbar definitions file.
2048 Where no additional xmin/xmax or ymin/ymax values are provided, the
2049 axis bounds default to use the vmin/vmax values provided.
2051 Returns
2052 -------
2053 cmap:
2054 Matplotlib colormap.
2055 levels:
2056 List of levels to use for plotting. For continuous plots the min and max
2057 should be taken as the range.
2058 norm:
2059 BoundaryNorm information.
2060 """
2061 if "difference" not in cube.long_name:
2062 if axis:
2063 levels = [0, 1]
2064 # Complete settings based on levels.
2065 return None, levels, None
2066 else:
2067 # Define the levels and colors.
2068 levels = [0, 1, 2]
2069 colors = ["white", "dodgerblue"]
2070 # Create a custom color map.
2071 cmap = mcolors.ListedColormap(colors)
2072 # Normalize the levels.
2073 norm = mcolors.BoundaryNorm(levels, cmap.N)
2074 logging.debug("Colourmap for %s.", cube.long_name)
2075 return cmap, levels, norm
2076 else:
2077 if axis:
2078 levels = [-1, 1]
2079 return None, levels, None
2080 else:
2081 # Search for if mask difference, set to +/- 0.5 as values plotted <
2082 # not <=.
2083 levels = [-2, -0.5, 0.5, 2]
2084 colors = ["goldenrod", "white", "teal"]
2085 cmap = mcolors.ListedColormap(colors)
2086 norm = mcolors.BoundaryNorm(levels, cmap.N)
2087 logging.debug("Colourmap for %s.", cube.long_name)
2088 return cmap, levels, norm
2091def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2092 """Get a custom colorbar for a cube in the Beaufort Scale.
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 Beaufort Scale 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 "difference" not in cube.long_name:
2120 if axis:
2121 levels = [0, 12]
2122 return None, levels, None
2123 else:
2124 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
2125 colors = [
2126 "black",
2127 (0, 0, 0.6),
2128 "blue",
2129 "cyan",
2130 "green",
2131 "yellow",
2132 (1, 0.5, 0),
2133 "red",
2134 "pink",
2135 "magenta",
2136 "purple",
2137 "maroon",
2138 "white",
2139 ]
2140 cmap = mcolors.ListedColormap(colors)
2141 norm = mcolors.BoundaryNorm(levels, cmap.N)
2142 logging.info("change colormap for Beaufort Scale colorbar.")
2143 return cmap, levels, norm
2144 else:
2145 if axis:
2146 levels = [-4, 4]
2147 return None, levels, None
2148 else:
2149 levels = [
2150 -3.5,
2151 -2.5,
2152 -1.5,
2153 -0.5,
2154 0.5,
2155 1.5,
2156 2.5,
2157 3.5,
2158 ]
2159 cmap = plt.get_cmap("bwr", 8)
2160 norm = mcolors.BoundaryNorm(levels, cmap.N)
2161 return cmap, levels, norm
2164def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
2165 """Return altered colourmap for temperature with change in units to Celsius.
2167 If "Celsius" appears anywhere in the name of a cube this function will be called.
2169 Parameters
2170 ----------
2171 cube: Cube
2172 Cube of variable for which the colorbar information is desired.
2173 cmap: Matplotlib colormap.
2174 levels: List
2175 List of levels to use for plotting. For continuous plots the min and max
2176 should be taken as the range.
2177 norm: BoundaryNorm.
2179 Returns
2180 -------
2181 cmap: Matplotlib colormap.
2182 levels: List
2183 List of levels to use for plotting. For continuous plots the min and max
2184 should be taken as the range.
2185 norm: BoundaryNorm.
2186 """
2187 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2188 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
2189 levels = np.array(levels)
2190 levels -= 273
2191 levels = levels.tolist()
2192 else:
2193 # Do nothing keep the existing colourbar attributes
2194 levels = levels
2195 cmap = cmap
2196 norm = norm
2197 return cmap, levels, norm
2200def _custom_colormap_probability(
2201 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
2202):
2203 """Get a custom colorbar for a probability cube.
2205 Specific variable ranges can be separately set in user-supplied definition
2206 for x- or y-axis limits, or indicate where automated range preferred.
2208 Parameters
2209 ----------
2210 cube: Cube
2211 Cube of variable with probability in name.
2212 axis: "x", "y", optional
2213 Select the levels for just this axis of a line plot. The min and max
2214 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2215 setting a universal range is not desirable (e.g. temperature), users
2216 can set ymin/ymax values to "auto" in the colorbar definitions file.
2217 Where no additional xmin/xmax or ymin/ymax values are provided, the
2218 axis bounds default to use the vmin/vmax values provided.
2220 Returns
2221 -------
2222 cmap:
2223 Matplotlib colormap.
2224 levels:
2225 List of levels to use for plotting. For continuous plots the min and max
2226 should be taken as the range.
2227 norm:
2228 BoundaryNorm information.
2229 """
2230 if axis:
2231 levels = [0, 1]
2232 return None, levels, None
2233 else:
2234 cmap = mcolors.ListedColormap(
2235 [
2236 "#FFFFFF",
2237 "#636363",
2238 "#e1dada",
2239 "#B5CAFF",
2240 "#8FB3FF",
2241 "#7F97FF",
2242 "#ABCF63",
2243 "#E8F59E",
2244 "#FFFA14",
2245 "#FFD121",
2246 "#FFA30A",
2247 ]
2248 )
2249 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]
2250 norm = mcolors.BoundaryNorm(levels, cmap.N)
2251 return cmap, levels, norm
2254def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
2255 """Return a custom colourmap for the current recipe."""
2256 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2257 if (
2258 any("surface_microphysical" in name for name in varnames)
2259 and "difference" not in cube.long_name
2260 and "mask" not in cube.long_name
2261 ):
2262 # Define the levels and colors
2263 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
2264 colors = [
2265 "w",
2266 (0, 0, 0.6),
2267 "b",
2268 "c",
2269 "g",
2270 "y",
2271 (1, 0.5, 0),
2272 "r",
2273 "pink",
2274 "m",
2275 "purple",
2276 "maroon",
2277 "gray",
2278 ]
2279 # Create a custom colormap
2280 cmap = mcolors.ListedColormap(colors)
2281 # Normalize the levels
2282 norm = mcolors.BoundaryNorm(levels, cmap.N)
2283 logging.info("change colormap for surface_microphysical variable colorbar.")
2284 else:
2285 # do nothing and keep existing colorbar attributes
2286 cmap = cmap
2287 levels = levels
2288 norm = norm
2289 return cmap, levels, norm
2292def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
2293 """Return custom colourmap for aviation colour state.
2295 If "aviation_colour_state" appears anywhere in the name of a cube
2296 this function will be called.
2298 Parameters
2299 ----------
2300 cube: Cube
2301 Cube of variable for which the colorbar information is desired.
2303 Returns
2304 -------
2305 cmap: Matplotlib colormap.
2306 levels: List
2307 List of levels to use for plotting. For continuous plots the min and max
2308 should be taken as the range.
2309 norm: BoundaryNorm.
2310 """
2311 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
2312 colors = [
2313 "#87ceeb",
2314 "#ffffff",
2315 "#8ced69",
2316 "#ffff00",
2317 "#ffd700",
2318 "#ffa500",
2319 "#fe3620",
2320 ]
2321 # Create a custom colormap
2322 cmap = mcolors.ListedColormap(colors)
2323 # Normalise the levels
2324 norm = mcolors.BoundaryNorm(levels, cmap.N)
2325 return cmap, levels, norm
2328def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
2329 """Return a custom colourmap for the current recipe."""
2330 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2331 if (
2332 any("visibility_in_air" in name for name in varnames)
2333 and "difference" not in cube.long_name
2334 and "mask" not in cube.long_name
2335 ):
2336 # Define the levels and colors (in km)
2337 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]
2338 norm = mcolors.BoundaryNorm(levels, cmap.N)
2339 colours = [
2340 "#8f00d6",
2341 "#d10000",
2342 "#ff9700",
2343 "#ffff00",
2344 "#00007f",
2345 "#6c9ccd",
2346 "#aae8ff",
2347 "#37a648",
2348 "#8edc64",
2349 "#c5ffc5",
2350 "#dcdcdc",
2351 "#ffffff",
2352 ]
2353 # Create a custom colormap
2354 cmap = mcolors.ListedColormap(colours)
2355 # Normalize the levels
2356 norm = mcolors.BoundaryNorm(levels, cmap.N)
2357 logging.info("change colormap for visibility_in_air variable colorbar.")
2358 else:
2359 # do nothing and keep existing colorbar attributes
2360 cmap = cmap
2361 levels = levels
2362 norm = norm
2363 return cmap, levels, norm
2366def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
2367 """Return number of models based on cube attributes."""
2368 model_names = list(
2369 filter(
2370 lambda x: x is not None,
2371 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
2372 )
2373 )
2374 if not model_names:
2375 logging.debug("Missing model names. Will assume single model.")
2376 return 1
2377 else:
2378 return len(model_names)
2381def _validate_cube_shape(
2382 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2383) -> None:
2384 """Check all cubes have a model name."""
2385 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2385 ↛ 2386line 2385 didn't jump to line 2386 because the condition on line 2385 was never true
2386 raise ValueError(
2387 f"The number of model names ({num_models}) should equal the number "
2388 f"of cubes ({len(cube)})."
2389 )
2392def _validate_cubes_coords(
2393 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2394) -> None:
2395 """Check same number of cubes as sequence coordinate for zip functions."""
2396 if len(cubes) != len(coords): 2396 ↛ 2397line 2396 didn't jump to line 2397 because the condition on line 2396 was never true
2397 raise ValueError(
2398 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2399 f"of sequence coordinates ({len(coords)})."
2400 f"Check that number of time entries in input data are consistent if "
2401 f"performing time-averaging steps prior to plotting outputs."
2402 )
2405####################
2406# Public functions #
2407####################
2410def spatial_contour_plot(
2411 cube: iris.cube.Cube,
2412 filename: str = None,
2413 sequence_coordinate: str = "time",
2414 stamp_coordinate: str = "realization",
2415 **kwargs,
2416) -> iris.cube.Cube:
2417 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2419 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2420 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2421 is present then postage stamp plots will be produced.
2423 Parameters
2424 ----------
2425 cube: Cube
2426 Iris cube of the data to plot. It should have two spatial dimensions,
2427 such as lat and lon, and may also have a another two dimension to be
2428 plotted sequentially and/or as postage stamp plots.
2429 filename: str, optional
2430 Name of the plot to write, used as a prefix for plot sequences. Defaults
2431 to the recipe name.
2432 sequence_coordinate: str, optional
2433 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2434 This coordinate must exist in the cube.
2435 stamp_coordinate: str, optional
2436 Coordinate about which to plot postage stamp plots. Defaults to
2437 ``"realization"``.
2439 Returns
2440 -------
2441 Cube
2442 The original cube (so further operations can be applied).
2444 Raises
2445 ------
2446 ValueError
2447 If the cube doesn't have the right dimensions.
2448 TypeError
2449 If the cube isn't a single cube.
2450 """
2451 _spatial_plot(
2452 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2453 )
2454 return cube
2457def spatial_pcolormesh_plot(
2458 cube: iris.cube.Cube,
2459 filename: str = None,
2460 sequence_coordinate: str = "time",
2461 stamp_coordinate: str = "realization",
2462 **kwargs,
2463) -> iris.cube.Cube:
2464 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2466 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2467 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2468 is present then postage stamp plots will be produced.
2470 This function is significantly faster than ``spatial_contour_plot``,
2471 especially at high resolutions, and should be preferred unless contiguous
2472 contour areas are important.
2474 Parameters
2475 ----------
2476 cube: Cube
2477 Iris cube of the data to plot. It should have two spatial dimensions,
2478 such as lat and lon, and may also have a another two dimension to be
2479 plotted sequentially and/or as postage stamp plots.
2480 filename: str, optional
2481 Name of the plot to write, used as a prefix for plot sequences. Defaults
2482 to the recipe name.
2483 sequence_coordinate: str, optional
2484 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2485 This coordinate must exist in the cube.
2486 stamp_coordinate: str, optional
2487 Coordinate about which to plot postage stamp plots. Defaults to
2488 ``"realization"``.
2490 Returns
2491 -------
2492 Cube
2493 The original cube (so further operations can be applied).
2495 Raises
2496 ------
2497 ValueError
2498 If the cube doesn't have the right dimensions.
2499 TypeError
2500 If the cube isn't a single cube.
2501 """
2502 _spatial_plot(
2503 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2504 )
2505 return cube
2508def spatial_multi_pcolormesh_plot(
2509 cube: iris.cube.Cube,
2510 overlay_cube: iris.cube.Cube,
2511 contour_cube: iris.cube.Cube,
2512 filename: str = None,
2513 sequence_coordinate: str = "time",
2514 stamp_coordinate: str = "realization",
2515 **kwargs,
2516) -> iris.cube.Cube:
2517 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2519 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2520 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2521 is present then postage stamp plots will be produced.
2523 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2525 If specified, contours of a contour_cube can be overplotted on top of those.
2527 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2529 This function is significantly faster than ``spatial_contour_plot``,
2530 especially at high resolutions, and should be preferred unless contiguous
2531 contour areas are important.
2533 Parameters
2534 ----------
2535 cube: Cube
2536 Iris cube of the data to plot. It should have two spatial dimensions,
2537 such as lat and lon, and may also have a another two dimension to be
2538 plotted sequentially and/or as postage stamp plots.
2539 overlay_cube: Cube
2540 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2541 such as lat and lon, and may also have a another two dimension to be
2542 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.
2543 contour_cube: Cube
2544 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,
2545 such as lat and lon, and may also have a another two dimension to be
2546 plotted sequentially and/or as postage stamp plots.
2547 filename: str, optional
2548 Name of the plot to write, used as a prefix for plot sequences. Defaults
2549 to the recipe name.
2550 sequence_coordinate: str, optional
2551 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2552 This coordinate must exist in the cube.
2553 stamp_coordinate: str, optional
2554 Coordinate about which to plot postage stamp plots. Defaults to
2555 ``"realization"``.
2557 Returns
2558 -------
2559 Cube
2560 The original cube (so further operations can be applied).
2562 Raises
2563 ------
2564 ValueError
2565 If the cube doesn't have the right dimensions.
2566 TypeError
2567 If the cube isn't a single cube.
2568 """
2569 _spatial_plot(
2570 "pcolormesh",
2571 cube,
2572 filename,
2573 sequence_coordinate,
2574 stamp_coordinate,
2575 overlay_cube=overlay_cube,
2576 contour_cube=contour_cube,
2577 )
2578 return cube, overlay_cube, contour_cube
2581# TODO: Expand function to handle ensemble data.
2582# line_coordinate: str, optional
2583# Coordinate about which to plot multiple lines. Defaults to
2584# ``"realization"``.
2585def plot_line_series(
2586 cube: iris.cube.Cube | iris.cube.CubeList,
2587 filename: str = None,
2588 series_coordinate: str = "time",
2589 # line_coordinate: str = "realization",
2590 **kwargs,
2591) -> iris.cube.Cube | iris.cube.CubeList:
2592 """Plot a line plot for the specified coordinate.
2594 The Cube or CubeList must be 1D.
2596 Parameters
2597 ----------
2598 iris.cube | iris.cube.CubeList
2599 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2600 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2601 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2602 filename: str, optional
2603 Name of the plot to write, used as a prefix for plot sequences. Defaults
2604 to the recipe name.
2605 series_coordinate: str, optional
2606 Coordinate about which to make a series. Defaults to ``"time"``. This
2607 coordinate must exist in the cube.
2609 Returns
2610 -------
2611 iris.cube.Cube | iris.cube.CubeList
2612 The original Cube or CubeList (so further operations can be applied).
2613 plotted data.
2615 Raises
2616 ------
2617 ValueError
2618 If the cubes don't have the right dimensions.
2619 TypeError
2620 If the cube isn't a Cube or CubeList.
2621 """
2622 # Ensure we have a name for the plot file.
2623 recipe_title = get_recipe_metadata().get("title", "Untitled")
2625 num_models = _get_num_models(cube)
2627 _validate_cube_shape(cube, num_models)
2629 # Iterate over all cubes and extract coordinate to plot.
2630 cubes = iter_maybe(cube)
2631 coords = []
2632 for cube in cubes:
2633 try:
2634 coords.append(cube.coord(series_coordinate))
2635 except iris.exceptions.CoordinateNotFoundError as err:
2636 raise ValueError(
2637 f"Cube must have a {series_coordinate} coordinate."
2638 ) from err
2639 if cube.ndim > 2 or not cube.coords("realization"):
2640 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2642 # Format the title and filename using plotted series coordinate
2643 nplot = 1
2644 seq_coord = coords[0]
2645 plot_title, plot_filename = _set_title_and_filename(
2646 seq_coord, nplot, recipe_title, filename
2647 )
2649 # Do the actual plotting.
2650 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, plot_title)
2652 # Add list of plots to plot metadata.
2653 plot_index = _append_to_plot_index([plot_filename])
2655 # Make a page to display the plots.
2656 _make_plot_html_page(plot_index)
2658 return cube
2661def plot_vertical_line_series(
2662 cubes: iris.cube.Cube | iris.cube.CubeList,
2663 filename: str = None,
2664 series_coordinate: str = "model_level_number",
2665 sequence_coordinate: str = "time",
2666 # line_coordinate: str = "realization",
2667 **kwargs,
2668) -> iris.cube.Cube | iris.cube.CubeList:
2669 """Plot a line plot against a type of vertical coordinate.
2671 The Cube or CubeList must be 1D.
2673 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2674 then a sequence of plots will be produced.
2676 Parameters
2677 ----------
2678 iris.cube | iris.cube.CubeList
2679 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2680 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2681 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2682 filename: str, optional
2683 Name of the plot to write, used as a prefix for plot sequences. Defaults
2684 to the recipe name.
2685 series_coordinate: str, optional
2686 Coordinate to plot on the y-axis. Can be ``pressure`` or
2687 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2688 for LFRic. Defaults to ``model_level_number``.
2689 This coordinate must exist in the cube.
2690 sequence_coordinate: str, optional
2691 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2692 This coordinate must exist in the cube.
2694 Returns
2695 -------
2696 iris.cube.Cube | iris.cube.CubeList
2697 The original Cube or CubeList (so further operations can be applied).
2698 Plotted data.
2700 Raises
2701 ------
2702 ValueError
2703 If the cubes doesn't have the right dimensions.
2704 TypeError
2705 If the cube isn't a Cube or CubeList.
2706 """
2707 # Ensure we have a name for the plot file.
2708 recipe_title = get_recipe_metadata().get("title", "Untitled")
2710 cubes = iter_maybe(cubes)
2711 # Initialise empty list to hold all data from all cubes in a CubeList
2712 all_data = []
2714 # Store min/max ranges for x range.
2715 x_levels = []
2717 num_models = _get_num_models(cubes)
2719 _validate_cube_shape(cubes, num_models)
2721 # Iterate over all cubes in cube or CubeList and plot.
2722 coords = []
2723 for cube in cubes:
2724 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2725 try:
2726 coords.append(cube.coord(series_coordinate))
2727 except iris.exceptions.CoordinateNotFoundError as err:
2728 raise ValueError(
2729 f"Cube must have a {series_coordinate} coordinate."
2730 ) from err
2732 try:
2733 if cube.ndim > 1 or not cube.coords("realization"): 2733 ↛ 2741line 2733 didn't jump to line 2741 because the condition on line 2733 was always true
2734 cube.coord(sequence_coordinate)
2735 except iris.exceptions.CoordinateNotFoundError as err:
2736 raise ValueError(
2737 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2738 ) from err
2740 # Get minimum and maximum from levels information.
2741 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2742 if levels is not None: 2742 ↛ 2746line 2742 didn't jump to line 2746 because the condition on line 2742 was always true
2743 x_levels.append(min(levels))
2744 x_levels.append(max(levels))
2745 else:
2746 all_data.append(cube.data)
2748 if len(x_levels) == 0: 2748 ↛ 2750line 2748 didn't jump to line 2750 because the condition on line 2748 was never true
2749 # Combine all data into a single NumPy array
2750 combined_data = np.concatenate(all_data)
2752 # Set the lower and upper limit for the x-axis to ensure all plots have
2753 # same range. This needs to read the whole cube over the range of the
2754 # sequence and if applicable postage stamp coordinate.
2755 vmin = np.floor(combined_data.min())
2756 vmax = np.ceil(combined_data.max())
2757 else:
2758 vmin = min(x_levels)
2759 vmax = max(x_levels)
2761 # Matching the slices (matching by seq coord point; it may happen that
2762 # evaluated models do not cover the same seq coord range, hence matching
2763 # necessary)
2764 def filter_cube_iterables(cube_iterables) -> bool:
2765 return len(cube_iterables) == len(coords)
2767 cube_iterables = filter(
2768 filter_cube_iterables,
2769 (
2770 iris.cube.CubeList(
2771 s
2772 for s in itertools.chain.from_iterable(
2773 cb.slices_over(sequence_coordinate) for cb in cubes
2774 )
2775 if s.coord(sequence_coordinate).points[0] == point
2776 )
2777 for point in sorted(
2778 set(
2779 itertools.chain.from_iterable(
2780 cb.coord(sequence_coordinate).points for cb in cubes
2781 )
2782 )
2783 )
2784 ),
2785 )
2787 # Create a plot for each value of the sequence coordinate.
2788 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2789 # similar sequence values. Passing a CubeList into the internal plotting function
2790 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2791 # or an iris.cube.CubeList.
2792 plot_index = []
2793 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2794 for cubes_slice in cube_iterables:
2795 # Format the coordinate value in a unit appropriate way.
2796 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2797 plot_title, plot_filename = _set_title_and_filename(
2798 seq_coord, nplot, recipe_title, filename
2799 )
2801 # Do the actual plotting.
2802 _plot_and_save_vertical_line_series(
2803 cubes_slice,
2804 coords,
2805 "realization",
2806 plot_filename,
2807 series_coordinate,
2808 title=plot_title,
2809 vmin=vmin,
2810 vmax=vmax,
2811 )
2812 plot_index.append(plot_filename)
2814 # Add list of plots to plot metadata.
2815 complete_plot_index = _append_to_plot_index(plot_index)
2817 # Make a page to display the plots.
2818 _make_plot_html_page(complete_plot_index)
2820 return cubes
2823def qq_plot(
2824 cubes: iris.cube.CubeList,
2825 coordinates: list[str],
2826 percentiles: list[float],
2827 model_names: list[str],
2828 filename: str = None,
2829 one_to_one: bool = True,
2830 **kwargs,
2831) -> iris.cube.CubeList:
2832 """Plot a Quantile-Quantile plot between two models for common time points.
2834 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2835 collapsed within the operator over all specified coordinates such as
2836 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2837 ensemble members to ensure a 1D cube (array).
2839 Parameters
2840 ----------
2841 cubes: iris.cube.CubeList
2842 Two cubes of the same variable with different models.
2843 coordinate: list[str]
2844 The list of coordinates to collapse over. This list should be
2845 every coordinate within the cube to result in a 1D cube around
2846 the percentile coordinate.
2847 percent: list[float]
2848 A list of percentiles to appear in the plot.
2849 model_names: list[str]
2850 A list of model names to appear on the axis of the plot.
2851 filename: str, optional
2852 Filename of the plot to write.
2853 one_to_one: bool, optional
2854 If True a 1:1 line is plotted; if False it is not. Default is True.
2856 Raises
2857 ------
2858 ValueError
2859 When the cubes are not compatible.
2861 Notes
2862 -----
2863 The quantile-quantile plot is a variant on the scatter plot representing
2864 two datasets by their quantiles (percentiles) for common time points.
2865 This plot does not use a theoretical distribution to compare against, but
2866 compares percentiles of two datasets. This plot does
2867 not use all raw data points, but plots the selected percentiles (quantiles) of
2868 each variable instead for the two datasets, thereby normalising the data for a
2869 direct comparison between the selected percentiles of the two dataset distributions.
2871 Quantile-quantile plots are valuable for comparing against
2872 observations and other models. Identical percentiles between the variables
2873 will lie on the one-to-one line implying the values correspond well to each
2874 other. Where there is a deviation from the one-to-one line a range of
2875 possibilities exist depending on how and where the data is shifted (e.g.,
2876 Wilks 2011 [Wilks2011]_).
2878 For distributions above the one-to-one line the distribution is left-skewed;
2879 below is right-skewed. A distinct break implies a bimodal distribution, and
2880 closer values/values further apart at the tails imply poor representation of
2881 the extremes.
2883 References
2884 ----------
2885 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2886 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2887 """
2888 # Check cubes using same functionality as the difference operator.
2889 if len(cubes) != 2:
2890 raise ValueError("cubes should contain exactly 2 cubes.")
2891 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2892 other: Cube = cubes.extract_cube(
2893 iris.Constraint(
2894 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2895 )
2896 )
2898 # Get spatial coord names.
2899 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2900 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2902 # Ensure cubes to compare are on common differencing grid.
2903 # This is triggered if either
2904 # i) latitude and longitude shapes are not the same. Note grid points
2905 # are not compared directly as these can differ through rounding
2906 # errors.
2907 # ii) or variables are known to often sit on different grid staggering
2908 # in different models (e.g. cell center vs cell edge), as is the case
2909 # for UM and LFRic comparisons.
2910 # In future greater choice of regridding method might be applied depending
2911 # on variable type. Linear regridding can in general be appropriate for smooth
2912 # variables. Care should be taken with interpretation of differences
2913 # given this dependency on regridding.
2914 if (
2915 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2916 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2917 ) or (
2918 base.long_name
2919 in [
2920 "eastward_wind_at_10m",
2921 "northward_wind_at_10m",
2922 "northward_wind_at_cell_centres",
2923 "eastward_wind_at_cell_centres",
2924 "zonal_wind_at_pressure_levels",
2925 "meridional_wind_at_pressure_levels",
2926 "potential_vorticity_at_pressure_levels",
2927 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2928 ]
2929 ):
2930 logging.debug(
2931 "Linear regridding base cube to other grid to compute differences"
2932 )
2933 base = regrid_onto_cube(base, other, method="Linear")
2935 # Extract just common time points.
2936 base, other = _extract_common_time_points(base, other)
2938 # Equalise attributes so we can merge.
2939 fully_equalise_attributes([base, other])
2940 logging.debug("Base: %s\nOther: %s", base, other)
2942 # Collapse cubes.
2943 base = collapse(
2944 base,
2945 coordinate=coordinates,
2946 method="PERCENTILE",
2947 additional_percent=percentiles,
2948 )
2949 other = collapse(
2950 other,
2951 coordinate=coordinates,
2952 method="PERCENTILE",
2953 additional_percent=percentiles,
2954 )
2956 # Ensure we have a name for the plot file.
2957 recipe_title = get_recipe_metadata().get("title", "Untitled")
2958 title = f"{recipe_title}"
2960 if filename is None:
2961 filename = slugify(recipe_title)
2963 # Add file extension.
2964 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2966 # Do the actual plotting on a scatter plot
2967 _plot_and_save_scatter_plot(
2968 base, other, plot_filename, title, one_to_one, model_names
2969 )
2971 # Add list of plots to plot metadata.
2972 plot_index = _append_to_plot_index([plot_filename])
2974 # Make a page to display the plots.
2975 _make_plot_html_page(plot_index)
2977 return iris.cube.CubeList([base, other])
2980def scatter_plot(
2981 cube_x: iris.cube.Cube | iris.cube.CubeList,
2982 cube_y: iris.cube.Cube | iris.cube.CubeList,
2983 filename: str = None,
2984 one_to_one: bool = True,
2985 **kwargs,
2986) -> iris.cube.CubeList:
2987 """Plot a scatter plot between two variables.
2989 Both cubes must be 1D.
2991 Parameters
2992 ----------
2993 cube_x: Cube | CubeList
2994 1 dimensional Cube of the data to plot on y-axis.
2995 cube_y: Cube | CubeList
2996 1 dimensional Cube of the data to plot on x-axis.
2997 filename: str, optional
2998 Filename of the plot to write.
2999 one_to_one: bool, optional
3000 If True a 1:1 line is plotted; if False it is not. Default is True.
3002 Returns
3003 -------
3004 cubes: CubeList
3005 CubeList of the original x and y cubes for further processing.
3007 Raises
3008 ------
3009 ValueError
3010 If the cube doesn't have the right dimensions and cubes not the same
3011 size.
3012 TypeError
3013 If the cube isn't a single cube.
3015 Notes
3016 -----
3017 Scatter plots are used for determining if there is a relationship between
3018 two variables. Positive relations have a slope going from bottom left to top
3019 right; Negative relations have a slope going from top left to bottom right.
3020 """
3021 # Iterate over all cubes in cube or CubeList and plot.
3022 for cube_iter in iter_maybe(cube_x):
3023 # Check cubes are correct shape.
3024 cube_iter = _check_single_cube(cube_iter)
3025 if cube_iter.ndim > 1:
3026 raise ValueError("cube_x must be 1D.")
3028 # Iterate over all cubes in cube or CubeList and plot.
3029 for cube_iter in iter_maybe(cube_y):
3030 # Check cubes are correct shape.
3031 cube_iter = _check_single_cube(cube_iter)
3032 if cube_iter.ndim > 1:
3033 raise ValueError("cube_y must be 1D.")
3035 # Ensure we have a name for the plot file.
3036 recipe_title = get_recipe_metadata().get("title", "Untitled")
3037 title = f"{recipe_title}"
3039 if filename is None:
3040 filename = slugify(recipe_title)
3042 # Add file extension.
3043 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3045 # Do the actual plotting.
3046 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3048 # Add list of plots to plot metadata.
3049 plot_index = _append_to_plot_index([plot_filename])
3051 # Make a page to display the plots.
3052 _make_plot_html_page(plot_index)
3054 return iris.cube.CubeList([cube_x, cube_y])
3057def vector_plot(
3058 cube_u: iris.cube.Cube,
3059 cube_v: iris.cube.Cube,
3060 filename: str = None,
3061 sequence_coordinate: str = "time",
3062 **kwargs,
3063) -> iris.cube.CubeList:
3064 """Plot a vector plot based on the input u and v components."""
3065 recipe_title = get_recipe_metadata().get("title", "Untitled")
3067 # Cubes must have a matching sequence coordinate.
3068 try:
3069 # Check that the u and v cubes have the same sequence coordinate.
3070 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3070 ↛ anywhereline 3070 didn't jump anywhere: it always raised an exception.
3071 raise ValueError("Coordinates do not match.")
3072 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3073 raise ValueError(
3074 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3075 ) from err
3077 # Create a plot for each value of the sequence coordinate.
3078 plot_index = []
3079 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3080 for cube_u_slice, cube_v_slice in zip(
3081 cube_u.slices_over(sequence_coordinate),
3082 cube_v.slices_over(sequence_coordinate),
3083 strict=True,
3084 ):
3085 # Format the coordinate value in a unit appropriate way.
3086 seq_coord = cube_u_slice.coord(sequence_coordinate)
3087 plot_title, plot_filename = _set_title_and_filename(
3088 seq_coord, nplot, recipe_title, filename
3089 )
3091 # Do the actual plotting.
3092 _plot_and_save_vector_plot(
3093 cube_u_slice,
3094 cube_v_slice,
3095 filename=plot_filename,
3096 title=plot_title,
3097 method="contourf",
3098 )
3099 plot_index.append(plot_filename)
3101 # Add list of plots to plot metadata.
3102 complete_plot_index = _append_to_plot_index(plot_index)
3104 # Make a page to display the plots.
3105 _make_plot_html_page(complete_plot_index)
3107 return iris.cube.CubeList([cube_u, cube_v])
3110def plot_histogram_series(
3111 cubes: iris.cube.Cube | iris.cube.CubeList,
3112 filename: str = None,
3113 sequence_coordinate: str = "time",
3114 stamp_coordinate: str = "realization",
3115 single_plot: bool = False,
3116 **kwargs,
3117) -> iris.cube.Cube | iris.cube.CubeList:
3118 """Plot a histogram plot for each vertical level provided.
3120 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3121 is present then a sequence of plots will be produced using the time slider
3122 functionality to scroll through histograms against time. If a
3123 stamp_coordinate is present then postage stamp plots will be produced. If
3124 stamp_coordinate and single_plot is True, all postage stamp plots will be
3125 plotted in a single plot instead of separate postage stamp plots.
3127 Parameters
3128 ----------
3129 cubes: Cube | iris.cube.CubeList
3130 Iris cube or CubeList of the data to plot. It should have a single dimension other
3131 than the stamp coordinate.
3132 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3133 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3134 filename: str, optional
3135 Name of the plot to write, used as a prefix for plot sequences. Defaults
3136 to the recipe name.
3137 sequence_coordinate: str, optional
3138 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3139 This coordinate must exist in the cube and will be used for the time
3140 slider.
3141 stamp_coordinate: str, optional
3142 Coordinate about which to plot postage stamp plots. Defaults to
3143 ``"realization"``.
3144 single_plot: bool, optional
3145 If True, all postage stamp plots will be plotted in a single plot. If
3146 False, each postage stamp plot will be plotted separately. Is only valid
3147 if stamp_coordinate exists and has more than a single point.
3149 Returns
3150 -------
3151 iris.cube.Cube | iris.cube.CubeList
3152 The original Cube or CubeList (so further operations can be applied).
3153 Plotted data.
3155 Raises
3156 ------
3157 ValueError
3158 If the cube doesn't have the right dimensions.
3159 TypeError
3160 If the cube isn't a Cube or CubeList.
3161 """
3162 recipe_title = get_recipe_metadata().get("title", "Untitled")
3164 cubes = iter_maybe(cubes)
3166 # Internal plotting function.
3167 plotting_func = _plot_and_save_histogram_series
3169 num_models = _get_num_models(cubes)
3171 _validate_cube_shape(cubes, num_models)
3173 # If several histograms are plotted with time as sequence_coordinate for the
3174 # time slider option.
3175 for cube in cubes:
3176 try:
3177 cube.coord(sequence_coordinate)
3178 except iris.exceptions.CoordinateNotFoundError as err:
3179 raise ValueError(
3180 f"Cube must have a {sequence_coordinate} coordinate."
3181 ) from err
3183 # Get minimum and maximum from levels information.
3184 levels = None
3185 for cube in cubes: 3185 ↛ 3201line 3185 didn't jump to line 3201 because the loop on line 3185 didn't complete
3186 # First check if user-specified "auto" range variable.
3187 # This maintains the value of levels as None, so proceed.
3188 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3189 if levels is None:
3190 break
3191 # If levels is changed, recheck to use the vmin,vmax or
3192 # levels-based ranges for histogram plots.
3193 _, levels, _ = _colorbar_map_levels(cube)
3194 logging.debug("levels: %s", levels)
3195 if levels is not None: 3195 ↛ 3185line 3195 didn't jump to line 3185 because the condition on line 3195 was always true
3196 vmin = min(levels)
3197 vmax = max(levels)
3198 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
3199 break
3201 if levels is None:
3202 vmin = min(cb.data.min() for cb in cubes)
3203 vmax = max(cb.data.max() for cb in cubes)
3205 # Make postage stamp plots if stamp_coordinate exists and has more than a
3206 # single point. If single_plot is True:
3207 # -- all postage stamp plots will be plotted in a single plot instead of
3208 # separate postage stamp plots.
3209 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3210 # produced per single model only
3211 if num_models == 1: 3211 ↛ 3224line 3211 didn't jump to line 3224 because the condition on line 3211 was always true
3212 if ( 3212 ↛ 3216line 3212 didn't jump to line 3216 because the condition on line 3212 was never true
3213 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3214 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3215 ):
3216 if single_plot:
3217 plotting_func = (
3218 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3219 )
3220 else:
3221 plotting_func = _plot_and_save_postage_stamp_histogram_series
3222 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3223 else:
3224 all_points = sorted(
3225 set(
3226 itertools.chain.from_iterable(
3227 cb.coord(sequence_coordinate).points for cb in cubes
3228 )
3229 )
3230 )
3231 all_slices = list(
3232 itertools.chain.from_iterable(
3233 cb.slices_over(sequence_coordinate) for cb in cubes
3234 )
3235 )
3236 # Matched slices (matched by seq coord point; it may happen that
3237 # evaluated models do not cover the same seq coord range, hence matching
3238 # necessary)
3239 cube_iterables = [
3240 iris.cube.CubeList(
3241 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3242 )
3243 for point in all_points
3244 ]
3246 plot_index = []
3247 nplot = np.size(cube.coord(sequence_coordinate).points)
3248 # Create a plot for each value of the sequence coordinate. Allowing for
3249 # multiple cubes in a CubeList to be plotted in the same plot for similar
3250 # sequence values. Passing a CubeList into the internal plotting function
3251 # for similar values of the sequence coordinate. cube_slice can be an
3252 # iris.cube.Cube or an iris.cube.CubeList.
3253 for cube_slice in cube_iterables:
3254 single_cube = cube_slice
3255 if isinstance(cube_slice, iris.cube.CubeList): 3255 ↛ 3256line 3255 didn't jump to line 3256 because the condition on line 3255 was never true
3256 single_cube = cube_slice[0]
3258 # Ensure valid stamp coordinate in cube dimensions
3259 if stamp_coordinate == "realization": 3259 ↛ 3262line 3259 didn't jump to line 3262 because the condition on line 3259 was always true
3260 stamp_coordinate = check_stamp_coordinate(single_cube)
3261 # Set plot titles and filename, based on sequence coordinate
3262 seq_coord = single_cube.coord(sequence_coordinate)
3263 # Use time coordinate in title and filename if single histogram output.
3264 if sequence_coordinate == "realization" and nplot == 1: 3264 ↛ 3265line 3264 didn't jump to line 3265 because the condition on line 3264 was never true
3265 seq_coord = single_cube.coord("time")
3266 plot_title, plot_filename = _set_title_and_filename(
3267 seq_coord, nplot, recipe_title, filename
3268 )
3270 # Do the actual plotting.
3271 plotting_func(
3272 cube_slice,
3273 filename=plot_filename,
3274 stamp_coordinate=stamp_coordinate,
3275 title=plot_title,
3276 vmin=vmin,
3277 vmax=vmax,
3278 )
3279 plot_index.append(plot_filename)
3281 # Add list of plots to plot metadata.
3282 complete_plot_index = _append_to_plot_index(plot_index)
3284 # Make a page to display the plots.
3285 _make_plot_html_page(complete_plot_index)
3287 return cubes
3290def plot_power_spectrum_series(
3291 cubes: iris.cube.Cube | iris.cube.CubeList,
3292 filename: str = None,
3293 sequence_coordinate: str = "time",
3294 stamp_coordinate: str = "realization",
3295 single_plot: bool = False,
3296 **kwargs,
3297) -> iris.cube.Cube | iris.cube.CubeList:
3298 """Plot a power spectrum plot for each vertical level provided.
3300 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
3301 is present then a sequence of plots will be produced using the time slider
3302 functionality to scroll through power spectrum against time. If a
3303 stamp_coordinate is present then postage stamp plots will be produced. If
3304 stamp_coordinate and single_plot is True, all postage stamp plots will be
3305 plotted in a single plot instead of separate postage stamp plots.
3307 Parameters
3308 ----------
3309 cubes: Cube | iris.cube.CubeList
3310 Iris cube or CubeList of the data to plot. It should have a single dimension other
3311 than the stamp coordinate.
3312 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3313 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3314 filename: str, optional
3315 Name of the plot to write, used as a prefix for plot sequences. Defaults
3316 to the recipe name.
3317 sequence_coordinate: str, optional
3318 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3319 This coordinate must exist in the cube and will be used for the time
3320 slider.
3321 stamp_coordinate: str, optional
3322 Coordinate about which to plot postage stamp plots. Defaults to
3323 ``"realization"``.
3324 single_plot: bool, optional
3325 If True, all postage stamp plots will be plotted in a single plot. If
3326 False, each postage stamp plot will be plotted separately. Is only valid
3327 if stamp_coordinate exists and has more than a single point.
3329 Returns
3330 -------
3331 iris.cube.Cube | iris.cube.CubeList
3332 The original Cube or CubeList (so further operations can be applied).
3333 Plotted data.
3335 Raises
3336 ------
3337 ValueError
3338 If the cube doesn't have the right dimensions.
3339 TypeError
3340 If the cube isn't a Cube or CubeList.
3341 """
3342 recipe_title = get_recipe_metadata().get("title", "Untitled")
3344 cubes = iter_maybe(cubes)
3346 # Internal plotting function.
3347 plotting_func = _plot_and_save_power_spectrum_series
3349 num_models = _get_num_models(cubes)
3351 _validate_cube_shape(cubes, num_models)
3353 # If several power spectra are plotted with time as sequence_coordinate for the
3354 # time slider option.
3355 for cube in cubes:
3356 try:
3357 cube.coord(sequence_coordinate)
3358 except iris.exceptions.CoordinateNotFoundError as err:
3359 raise ValueError(
3360 f"Cube must have a {sequence_coordinate} coordinate."
3361 ) from err
3363 # Make postage stamp plots if stamp_coordinate exists and has more than a
3364 # single point. If single_plot is True:
3365 # -- all postage stamp plots will be plotted in a single plot instead of
3366 # separate postage stamp plots.
3367 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3368 # produced per single model only
3369 if num_models == 1: 3369 ↛ 3382line 3369 didn't jump to line 3382 because the condition on line 3369 was always true
3370 if ( 3370 ↛ 3374line 3370 didn't jump to line 3374 because the condition on line 3370 was never true
3371 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3372 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3373 ):
3374 if single_plot:
3375 plotting_func = (
3376 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3377 )
3378 else:
3379 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3380 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3381 else:
3382 all_points = sorted(
3383 set(
3384 itertools.chain.from_iterable(
3385 cb.coord(sequence_coordinate).points for cb in cubes
3386 )
3387 )
3388 )
3389 all_slices = list(
3390 itertools.chain.from_iterable(
3391 cb.slices_over(sequence_coordinate) for cb in cubes
3392 )
3393 )
3394 # Matched slices (matched by seq coord point; it may happen that
3395 # evaluated models do not cover the same seq coord range, hence matching
3396 # necessary)
3397 cube_iterables = [
3398 iris.cube.CubeList(
3399 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3400 )
3401 for point in all_points
3402 ]
3404 plot_index = []
3405 nplot = np.size(cube.coord(sequence_coordinate).points)
3406 # Create a plot for each value of the sequence coordinate. Allowing for
3407 # multiple cubes in a CubeList to be plotted in the same plot for similar
3408 # sequence values. Passing a CubeList into the internal plotting function
3409 # for similar values of the sequence coordinate. cube_slice can be an
3410 # iris.cube.Cube or an iris.cube.CubeList.
3411 for cube_slice in cube_iterables:
3412 single_cube = cube_slice
3413 if isinstance(cube_slice, iris.cube.CubeList): 3413 ↛ 3414line 3413 didn't jump to line 3414 because the condition on line 3413 was never true
3414 single_cube = cube_slice[0]
3416 # Set stamp coordinate
3417 if stamp_coordinate == "realization": 3417 ↛ 3420line 3417 didn't jump to line 3420 because the condition on line 3417 was always true
3418 stamp_coordinate = check_stamp_coordinate(single_cube)
3419 # Set plot title and filenames based on sequence values
3420 seq_coord = single_cube.coord(sequence_coordinate)
3421 plot_title, plot_filename = _set_title_and_filename(
3422 seq_coord, nplot, recipe_title, filename
3423 )
3425 # Do the actual plotting.
3426 plotting_func(
3427 cube_slice,
3428 filename=plot_filename,
3429 stamp_coordinate=stamp_coordinate,
3430 title=plot_title,
3431 )
3432 plot_index.append(plot_filename)
3434 # Add list of plots to plot metadata.
3435 complete_plot_index = _append_to_plot_index(plot_index)
3437 # Make a page to display the plots.
3438 _make_plot_html_page(complete_plot_index)
3440 return cubes
3443def _DCT_ps(y_3d):
3444 """Calculate power spectra for regional domains.
3446 Parameters
3447 ----------
3448 y_3d: 3D array
3449 3 dimensional array to calculate spectrum for.
3450 (2D field data with 3rd dimension of time)
3452 Returns: ps_array
3453 Array of power spectra values calculated for input field (for each time)
3455 Method for regional domains:
3456 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3457 as described in Denis et al 2002 [Denis_etal_2002]_.
3459 References
3460 ----------
3461 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3462 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3463 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3464 Monthly Weather Review, Vol. 130, 1812-1828
3465 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3466 """
3467 Nt, Ny, Nx = y_3d.shape
3469 # Max coefficient
3470 Nmin = min(Nx - 1, Ny - 1)
3472 # Create alpha matrix (of wavenumbers)
3473 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3475 # Prepare output array
3476 ps_array = np.zeros((Nt, Nmin))
3478 # Loop over time to get spectrum for each time.
3479 for t in range(Nt):
3480 y_2d = y_3d[t]
3482 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3483 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3484 # cosine basis functions at different spatial frequencies.
3486 # normalise spectrum to allow comparison between models.
3487 fkk = fft.dctn(y_2d, norm="ortho")
3489 # Normalise fkk
3490 fkk = fkk / np.sqrt(Ny * Nx)
3492 # calculate variance of spectral coefficient
3493 sigma_2 = fkk**2 / Nx / Ny
3495 # Group ellipses of alphas into the same wavenumber k/Nmin
3496 for k in range(1, Nmin + 1):
3497 alpha = k / Nmin
3498 alpha_p1 = (k + 1) / Nmin
3500 # Sum up elements matching k
3501 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3502 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3504 return ps_array
3507def _create_alpha_matrix(Ny, Nx):
3508 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3510 Parameters
3511 ----------
3512 Ny, Nx:
3513 Dimensions of the 2D field for which the power spectra is calculated. Used to
3514 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3515 single-scale parameter.
3517 Returns: alpha_matrix
3518 normalisation of 2D wavenumber axes, transforming the spectral domain into
3519 an elliptic coordinate system.
3521 """
3522 # Create x_indices: each row is [1, 2, ..., Nx]
3523 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3525 # Create y_indices: each column is [1, 2, ..., Ny]
3526 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3528 # Compute alpha_matrix
3529 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3531 return alpha_matrix