Coverage for src / CSET / operators / plot.py: 83%
1085 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-13 07:38 +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):
286 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)
287 return cmap, levels, norm
289 # If no valid colormap has been defined, use defaults and return.
290 if not cmap:
291 logging.warning("No colorbar definition exists for %s.", cube.name())
292 cmap, levels, norm = mpl.colormaps["viridis"], None, None
293 return cmap, levels, norm
295 # Test if pressure-level specific settings are provided for cube.
296 if pressure_level:
297 try:
298 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]
299 except KeyError:
300 logging.debug(
301 "%s has no colorbar definition for pressure level %s.",
302 varname,
303 pressure_level,
304 )
306 # Check for availability of x-axis or y-axis user-specific overrides
307 # for setting level bounds for line plot types and return just levels.
308 # Line plots do not need a colormap, and just use the data range.
309 if axis:
310 if axis == "x":
311 try:
312 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]
313 except KeyError:
314 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
315 if axis == "y":
316 try:
317 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]
318 except KeyError:
319 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
320 # Check if user-specified auto-scaling for this variable
321 if vmin == "auto" or vmax == "auto":
322 levels = None
323 else:
324 levels = [vmin, vmax]
325 return None, levels, None
326 # Get and use the colorbar levels for this variable if spatial or histogram.
327 else:
328 try:
329 levels = var_colorbar["levels"]
330 # Use discrete bins when levels are specified, rather
331 # than a smooth range.
332 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
333 logging.debug("Using levels for %s colorbar.", varname)
334 logging.info("Using levels: %s", levels)
335 except KeyError:
336 # Get the range for this variable.
337 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
338 logging.debug("Using min and max for %s colorbar.", varname)
339 # Calculate levels from range.
340 if vmin == "auto" or vmax == "auto": 340 ↛ 341line 340 didn't jump to line 341 because the condition on line 340 was never true
341 levels = None
342 else:
343 levels = np.linspace(vmin, vmax, 101)
344 norm = None
346 # Overwrite cmap, levels and norm for specific variables that
347 # require custom colorbar_map as these can not be defined in the
348 # JSON file.
349 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)
350 cmap, levels, norm = _custom_colourmap_visibility_in_air(
351 cube, cmap, levels, norm
352 )
353 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)
354 return cmap, levels, norm
357def _setup_spatial_map(
358 cube: iris.cube.Cube,
359 figure,
360 cmap,
361 grid_size: tuple[int, int] | None = None,
362 subplot: int | None = None,
363):
364 """Define map projections, extent and add coastlines and borderlines for spatial plots.
366 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
367 is specified, and map extent defined based on the input data.
369 Parameters
370 ----------
371 cube: Cube
372 2 dimensional (lat and lon) Cube of the data to plot.
373 figure:
374 Matplotlib Figure object holding all plot elements.
375 cmap:
376 Matplotlib colormap.
377 grid_size: (int, int), optional
378 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
379 subplot: int, optional
380 Subplot index if multiple spatial subplots in figure.
382 Returns
383 -------
384 axes:
385 Matplotlib GeoAxes definition.
386 """
387 # Identify min/max plot bounds.
388 try:
389 lat_axis, lon_axis = get_cube_yxcoordname(cube)
390 x1 = np.min(cube.coord(lon_axis).points)
391 x2 = np.max(cube.coord(lon_axis).points)
392 y1 = np.min(cube.coord(lat_axis).points)
393 y2 = np.max(cube.coord(lat_axis).points)
395 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
396 if np.abs(x2 - x1) > 180.0:
397 x1 = x1 - 180.0
398 x2 = x2 - 180.0
399 logging.debug("Adjusting plot bounds to fit global extent.")
401 # Consider map projection orientation.
402 # Adapting orientation enables plotting across international dateline.
403 # Users can adapt the default central_longitude if alternative projections views.
404 if x2 > 180.0 or x1 < -180.0:
405 central_longitude = 180.0
406 else:
407 central_longitude = 0.0
409 # Define spatial map projection.
410 coord_system = cube.coord(lat_axis).coord_system
411 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
412 # Define rotated pole map projection for rotated pole inputs.
413 projection = ccrs.RotatedPole(
414 pole_longitude=coord_system.grid_north_pole_longitude,
415 pole_latitude=coord_system.grid_north_pole_latitude,
416 central_rotated_longitude=central_longitude,
417 )
418 crs = projection
419 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 419 ↛ 421line 419 didn't jump to line 421 because the condition on line 419 was never true
420 # Define Transverse Mercator projection for TM inputs.
421 projection = ccrs.TransverseMercator(
422 central_longitude=coord_system.longitude_of_central_meridian,
423 central_latitude=coord_system.latitude_of_projection_origin,
424 false_easting=coord_system.false_easting,
425 false_northing=coord_system.false_northing,
426 scale_factor=coord_system.scale_factor_at_central_meridian,
427 )
428 crs = projection
429 else:
430 # Define regular map projection for non-rotated pole inputs.
431 # Alternatives might include e.g. for global model outputs:
432 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
433 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
434 projection = ccrs.PlateCarree(central_longitude=central_longitude)
435 crs = ccrs.PlateCarree()
437 # Define axes for plot (or subplot) with required map projection.
438 if subplot is not None:
439 axes = figure.add_subplot(
440 grid_size[0], grid_size[1], subplot, projection=projection
441 )
442 else:
443 axes = figure.add_subplot(projection=projection)
445 # Add coastlines and borderlines if cube contains x and y map coordinates.
446 # Avoid adding lines for masked data or specific fixed ancillary spatial plots.
447 if iris.util.is_masked(cube.data) or any( 447 ↛ 450line 447 didn't jump to line 450 because the condition on line 447 was never true
448 name in cube.name() for name in ["land_", "orography", "altitude"]
449 ):
450 pass
451 else:
452 if cmap.name in ["viridis", "Greys"]:
453 coastcol = "magenta"
454 else:
455 coastcol = "black"
456 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
457 axes.coastlines(resolution="10m", color=coastcol)
458 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
460 # Add gridlines.
461 gl = axes.gridlines(
462 alpha=0.3,
463 draw_labels=True,
464 dms=False,
465 x_inline=False,
466 y_inline=False,
467 )
468 gl.top_labels = False
469 gl.right_labels = False
470 if subplot:
471 gl.bottom_labels = False
472 gl.left_labels = False
473 if subplot % grid_size[1] == 1:
474 gl.left_labels = True
475 if subplot > ((grid_size[0] - 1) * grid_size[1]): 475 ↛ 480line 475 didn't jump to line 480 because the condition on line 475 was always true
476 gl.bottom_labels = True
478 # If is lat/lon spatial map, fix extent to keep plot tight.
479 # Specifying crs within set_extent helps ensure only data region is shown.
480 if isinstance(coord_system, iris.coord_systems.GeogCS):
481 axes.set_extent([x1, x2, y1, y2], crs=crs)
483 except ValueError:
484 # Skip if not both x and y map coordinates.
485 axes = figure.gca()
486 pass
488 return axes
491def _get_plot_resolution() -> int:
492 """Get resolution of rasterised plots in pixels per inch."""
493 return get_recipe_metadata().get("plot_resolution", 100)
496def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
497 """Return title and filename based on start and end points or bounds."""
498 if use_bounds and seq_coord.has_bounds():
499 vals = seq_coord.bounds.flatten()
500 else:
501 vals = seq_coord.points
502 start = seq_coord.units.title(vals[0])
503 end = seq_coord.units.title(vals[-1])
505 if start == end:
506 sequence_title = f"\n [{start}]"
507 sequence_fname = f"_{filename_slugify(start)}"
508 else:
509 sequence_title = f"\n [{start} to {end}]"
510 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
512 # Do not include time if coord set to zero.
513 if (
514 seq_coord.units == "hours since 0001-01-01 00:00:00"
515 and vals[0] == 0
516 and vals[-1] == 0
517 ):
518 sequence_title = ""
519 sequence_fname = ""
521 return sequence_title, sequence_fname
524def _set_title_and_filename(
525 seq_coord: iris.coords.Coord,
526 nplot: int,
527 recipe_title: str,
528 filename: str,
529):
530 """Set plot title and filename based on cube coordinate.
532 Parameters
533 ----------
534 sequence_coordinate: iris.coords.Coord
535 Coordinate about which to make a plot sequence.
536 nplot: int
537 Number of output plots to generate - controls title/naming.
538 recipe_title: str
539 Default plot title, potentially to update.
540 filename: str
541 Input plot filename, potentially to update.
543 Returns
544 -------
545 plot_title: str
546 Output formatted plot title string, based on plotted data.
547 plot_filename: str
548 Output formatted plot filename string.
549 """
550 ndim = seq_coord.ndim
551 npoints = np.size(seq_coord.points)
552 sequence_title = ""
553 sequence_fname = ""
555 # Case 1: Multiple dimension sequence input - list number of aggregated cases
556 # (e.g. aggregation histogram plots)
557 if ndim > 1:
558 ncase = np.shape(seq_coord)[0]
559 sequence_title = f"\n [{ncase} cases]"
560 sequence_fname = f"_{ncase}cases"
562 # Case 2: Single dimension input
563 else:
564 # Single sequence point
565 if npoints == 1:
566 if nplot > 1:
567 # Default labels for sequence inputs
568 sequence_value = seq_coord.units.title(seq_coord.points[0])
569 sequence_title = f"\n [{sequence_value}]"
570 sequence_fname = f"_{filename_slugify(sequence_value)}"
571 else:
572 # Aggregated attribute available where input collapsed over aggregation
573 try:
574 ncase = seq_coord.attributes["number_reference_times"]
575 sequence_title = f"\n [{ncase} cases]"
576 sequence_fname = f"_{ncase}cases"
577 except KeyError:
578 sequence_title, sequence_fname = _get_start_end_strings(
579 seq_coord, use_bounds=seq_coord.has_bounds()
580 )
581 # Multiple sequence (e.g. time) points
582 else:
583 sequence_title, sequence_fname = _get_start_end_strings(
584 seq_coord, use_bounds=False
585 )
587 # Set plot title and filename
588 plot_title = f"{recipe_title}{sequence_title}"
590 # Set plot filename, defaulting to user input if provided.
591 if filename is None:
592 filename = slugify(recipe_title)
593 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
594 else:
595 if nplot > 1:
596 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
597 else:
598 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
600 return plot_title, plot_filename
603def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
604 """Control postage stamp plot output titles based on stamp coordinate."""
605 if stamp_coord.name() == "realization":
606 mtitle = "Member"
607 else:
608 mtitle = stamp_coord.name().capitalize()
610 if stamp_coord.name() == "time":
611 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
612 else:
613 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
615 return mtitle
618def _plot_and_save_spatial_plot(
619 cube: iris.cube.Cube,
620 filename: str,
621 title: str,
622 method: Literal["contourf", "pcolormesh"],
623 overlay_cube: iris.cube.Cube | None = None,
624 contour_cube: iris.cube.Cube | None = None,
625 **kwargs,
626):
627 """Plot and save a spatial plot.
629 Parameters
630 ----------
631 cube: Cube
632 2 dimensional (lat and lon) Cube of the data to plot.
633 filename: str
634 Filename of the plot to write.
635 title: str
636 Plot title.
637 method: "contourf" | "pcolormesh"
638 The plotting method to use.
639 overlay_cube: Cube, optional
640 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
641 contour_cube: Cube, optional
642 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
643 """
644 # Setup plot details, size, resolution, etc.
645 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
647 # Specify the color bar
648 cmap, levels, norm = _colorbar_map_levels(cube)
650 # If overplotting, set required colorbars
651 if overlay_cube:
652 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
653 if contour_cube:
654 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
656 # Setup plot map projection, extent and coastlines and borderlines.
657 axes = _setup_spatial_map(cube, fig, cmap)
659 # Plot the field.
660 if method == "contourf":
661 # Filled contour plot of the field.
662 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
663 elif method == "pcolormesh":
664 try:
665 vmin = min(levels)
666 vmax = max(levels)
667 except TypeError:
668 vmin, vmax = None, None
669 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
670 # if levels are defined.
671 if norm is not None:
672 vmin = None
673 vmax = None
674 logging.debug("Plotting using defined levels.")
675 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
676 else:
677 raise ValueError(f"Unknown plotting method: {method}")
679 # Overplot overlay field, if required
680 if overlay_cube:
681 try:
682 over_vmin = min(over_levels)
683 over_vmax = max(over_levels)
684 except TypeError:
685 over_vmin, over_vmax = None, None
686 if over_norm is not None: 686 ↛ 687line 686 didn't jump to line 687 because the condition on line 686 was never true
687 over_vmin = None
688 over_vmax = None
689 overlay = iplt.pcolormesh(
690 overlay_cube,
691 cmap=over_cmap,
692 norm=over_norm,
693 alpha=0.8,
694 vmin=over_vmin,
695 vmax=over_vmax,
696 )
697 # Overplot contour field, if required, with contour labelling.
698 if contour_cube:
699 contour = iplt.contour(
700 contour_cube,
701 colors="darkgray",
702 levels=cntr_levels,
703 norm=cntr_norm,
704 alpha=0.5,
705 linestyles="--",
706 linewidths=1,
707 )
708 plt.clabel(contour)
710 # Check to see if transect, and if so, adjust y axis.
711 if is_transect(cube):
712 if "pressure" in [coord.name() for coord in cube.coords()]:
713 axes.invert_yaxis()
714 axes.set_yscale("log")
715 axes.set_ylim(1100, 100)
716 # If both model_level_number and level_height exists, iplt can construct
717 # plot as a function of height above orography (NOT sea level).
718 elif {"model_level_number", "level_height"}.issubset( 718 ↛ 723line 718 didn't jump to line 723 because the condition on line 718 was always true
719 {coord.name() for coord in cube.coords()}
720 ):
721 axes.set_yscale("log")
723 axes.set_title(
724 f"{title}\n"
725 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
726 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
727 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
728 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
729 fontsize=16,
730 )
732 # Inset code
733 axins = inset_axes(
734 axes,
735 width="20%",
736 height="20%",
737 loc="upper right",
738 axes_class=GeoAxes,
739 axes_kwargs={"map_projection": ccrs.PlateCarree()},
740 )
742 # Slightly transparent to reduce plot blocking.
743 axins.patch.set_alpha(0.4)
745 axins.coastlines(resolution="50m")
746 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
748 SLat, SLon, ELat, ELon = (
749 float(coord) for coord in cube.attributes["transect_coords"].split("_")
750 )
752 # Draw line between them
753 axins.plot(
754 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
755 )
757 # Plot points (note: lon, lat order for Cartopy)
758 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
759 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
761 lon_min, lon_max = sorted([SLon, ELon])
762 lat_min, lat_max = sorted([SLat, ELat])
764 # Midpoints
765 lon_mid = (lon_min + lon_max) / 2
766 lat_mid = (lat_min + lat_max) / 2
768 # Maximum half-range
769 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
770 if half_range == 0: # points identical → provide small default 770 ↛ 774line 770 didn't jump to line 774 because the condition on line 770 was always true
771 half_range = 1
773 # Set square extent
774 axins.set_extent(
775 [
776 lon_mid - half_range,
777 lon_mid + half_range,
778 lat_mid - half_range,
779 lat_mid + half_range,
780 ],
781 crs=ccrs.PlateCarree(),
782 )
784 # Ensure square aspect
785 axins.set_aspect("equal")
787 else:
788 # Add title.
789 axes.set_title(title, fontsize=16)
791 # Adjust padding if spatial plot or transect
792 if is_transect(cube):
793 yinfopad = -0.1
794 ycbarpad = 0.1
795 else:
796 yinfopad = 0.01
797 ycbarpad = 0.042
799 # Add watermark with min/max/mean. Currently not user togglable.
800 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
801 axes.annotate(
802 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
803 xy=(0.025, yinfopad),
804 xycoords="axes fraction",
805 xytext=(-5, 5),
806 textcoords="offset points",
807 ha="left",
808 va="bottom",
809 size=11,
810 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
811 )
813 # Add secondary colour bar for overlay_cube field if required.
814 if overlay_cube:
815 cbarB = fig.colorbar(
816 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
817 )
818 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
819 # add ticks and tick_labels for every levels if less than 20 levels exist
820 if over_levels is not None and len(over_levels) < 20: 820 ↛ 821line 820 didn't jump to line 821 because the condition on line 820 was never true
821 cbarB.set_ticks(over_levels)
822 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
823 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
824 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
825 logging.debug("Set secondary colorbar ticks and labels.")
827 # Add main colour bar.
828 cbar = fig.colorbar(
829 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
830 )
832 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
833 # add ticks and tick_labels for every levels if less than 20 levels exist
834 if levels is not None and len(levels) < 20:
835 cbar.set_ticks(levels)
836 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
837 if "rainfall" or "snowfall" or "visibility" in cube.name(): 837 ↛ 839line 837 didn't jump to line 839 because the condition on line 837 was always true
838 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
839 logging.debug("Set colorbar ticks and labels.")
841 # Save plot.
842 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
843 logging.info("Saved spatial plot to %s", filename)
844 plt.close(fig)
847def _plot_and_save_postage_stamp_spatial_plot(
848 cube: iris.cube.Cube,
849 filename: str,
850 stamp_coordinate: str,
851 title: str,
852 method: Literal["contourf", "pcolormesh"],
853 overlay_cube: iris.cube.Cube | None = None,
854 contour_cube: iris.cube.Cube | None = None,
855 **kwargs,
856):
857 """Plot postage stamp spatial plots from an ensemble.
859 Parameters
860 ----------
861 cube: Cube
862 Iris cube of data to be plotted. It must have the stamp coordinate.
863 filename: str
864 Filename of the plot to write.
865 stamp_coordinate: str
866 Coordinate that becomes different plots.
867 method: "contourf" | "pcolormesh"
868 The plotting method to use.
869 overlay_cube: Cube, optional
870 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
871 contour_cube: Cube, optional
872 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
874 Raises
875 ------
876 ValueError
877 If the cube doesn't have the right dimensions.
878 """
879 # Use the smallest square grid that will fit the members.
880 nmember = len(cube.coord(stamp_coordinate).points)
881 grid_rows = int(math.sqrt(nmember))
882 grid_size = math.ceil(nmember / grid_rows)
884 fig = plt.figure(
885 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
886 )
888 # Specify the color bar
889 cmap, levels, norm = _colorbar_map_levels(cube)
890 # If overplotting, set required colorbars
891 if overlay_cube: 891 ↛ 892line 891 didn't jump to line 892 because the condition on line 891 was never true
892 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
893 if contour_cube: 893 ↛ 894line 893 didn't jump to line 894 because the condition on line 893 was never true
894 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
896 # Make a subplot for each member.
897 for member, subplot in zip(
898 cube.slices_over(stamp_coordinate),
899 range(1, grid_size * grid_rows + 1),
900 strict=False,
901 ):
902 # Setup subplot map projection, extent and coastlines and borderlines.
903 axes = _setup_spatial_map(
904 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
905 )
906 if method == "contourf":
907 # Filled contour plot of the field.
908 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
909 elif method == "pcolormesh":
910 if levels is not None:
911 vmin = min(levels)
912 vmax = max(levels)
913 else:
914 raise TypeError("Unknown vmin and vmax range.")
915 vmin, vmax = None, None
916 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
917 # if levels are defined.
918 if norm is not None: 918 ↛ 919line 918 didn't jump to line 919 because the condition on line 918 was never true
919 vmin = None
920 vmax = None
921 # pcolormesh plot of the field.
922 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
923 else:
924 raise ValueError(f"Unknown plotting method: {method}")
926 # Overplot overlay field, if required
927 if overlay_cube: 927 ↛ 928line 927 didn't jump to line 928 because the condition on line 927 was never true
928 try:
929 over_vmin = min(over_levels)
930 over_vmax = max(over_levels)
931 except TypeError:
932 over_vmin, over_vmax = None, None
933 if over_norm is not None:
934 over_vmin = None
935 over_vmax = None
936 iplt.pcolormesh(
937 overlay_cube[member.coord(stamp_coordinate).points[0]],
938 cmap=over_cmap,
939 norm=over_norm,
940 alpha=0.6,
941 vmin=over_vmin,
942 vmax=over_vmax,
943 )
944 # Overplot contour field, if required
945 if contour_cube: 945 ↛ 946line 945 didn't jump to line 946 because the condition on line 945 was never true
946 iplt.contour(
947 contour_cube[member.coord(stamp_coordinate).points[0]],
948 colors="darkgray",
949 levels=cntr_levels,
950 norm=cntr_norm,
951 alpha=0.6,
952 linestyles="--",
953 linewidths=1,
954 )
955 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
956 axes.set_title(f"{mtitle}")
958 # Put the shared colorbar in its own axes.
959 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
960 colorbar = fig.colorbar(
961 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
962 )
963 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
965 # Overall figure title.
966 fig.suptitle(title, fontsize=16)
968 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
969 logging.info("Saved contour postage stamp plot to %s", filename)
970 plt.close(fig)
973def _plot_and_save_line_series(
974 cubes: iris.cube.CubeList,
975 coords: list[iris.coords.Coord],
976 ensemble_coord: str,
977 filename: str,
978 title: str,
979 **kwargs,
980):
981 """Plot and save a 1D line series.
983 Parameters
984 ----------
985 cubes: Cube or CubeList
986 Cube or CubeList containing the cubes to plot on the y-axis.
987 coords: list[Coord]
988 Coordinates to plot on the x-axis, one per cube.
989 ensemble_coord: str
990 Ensemble coordinate in the cube.
991 filename: str
992 Filename of the plot to write.
993 title: str
994 Plot title.
995 """
996 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
998 model_colors_map = _get_model_colors_map(cubes)
1000 # Store min/max ranges.
1001 y_levels = []
1003 # Check match-up across sequence coords gives consistent sizes
1004 _validate_cubes_coords(cubes, coords)
1006 for cube, coord in zip(cubes, coords, strict=True):
1007 label = None
1008 color = "black"
1009 if model_colors_map:
1010 label = cube.attributes.get("model_name")
1011 color = model_colors_map.get(label)
1012 for cube_slice in cube.slices_over(ensemble_coord):
1013 # Label with (control) if part of an ensemble or not otherwise.
1014 if cube_slice.coord(ensemble_coord).points == [0]:
1015 iplt.plot(
1016 coord,
1017 cube_slice,
1018 color=color,
1019 marker="o",
1020 ls="-",
1021 lw=3,
1022 label=f"{label} (control)"
1023 if len(cube.coord(ensemble_coord).points) > 1
1024 else label,
1025 )
1026 # Label with (perturbed) if part of an ensemble and not the control.
1027 else:
1028 iplt.plot(
1029 coord,
1030 cube_slice,
1031 color=color,
1032 ls="-",
1033 lw=1.5,
1034 alpha=0.75,
1035 label=f"{label} (member)",
1036 )
1038 # Calculate the global min/max if multiple cubes are given.
1039 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1040 if levels is not None: 1040 ↛ 1041line 1040 didn't jump to line 1041 because the condition on line 1040 was never true
1041 y_levels.append(min(levels))
1042 y_levels.append(max(levels))
1044 # Get the current axes.
1045 ax = plt.gca()
1047 # Add some labels and tweak the style.
1048 # check if cubes[0] works for single cube if not CubeList
1049 if coords[0].name() == "time":
1050 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1051 else:
1052 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1053 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1054 ax.set_title(title, fontsize=16)
1056 ax.ticklabel_format(axis="y", useOffset=False)
1057 ax.tick_params(axis="x", labelrotation=15)
1058 ax.tick_params(axis="both", labelsize=12)
1060 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1061 if y_levels: 1061 ↛ 1062line 1061 didn't jump to line 1062 because the condition on line 1061 was never true
1062 ax.set_ylim(min(y_levels), max(y_levels))
1063 # Add zero line.
1064 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1065 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1066 logging.debug(
1067 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1068 )
1069 else:
1070 ax.autoscale()
1072 # Add gridlines
1073 ax.grid(linestyle="--", color="grey", linewidth=1)
1074 # Ientify unique labels for legend
1075 handles = list(
1076 {
1077 label: handle
1078 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1079 }.values()
1080 )
1081 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1083 # Save plot.
1084 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1085 logging.info("Saved line plot to %s", filename)
1086 plt.close(fig)
1089def _plot_and_save_vertical_line_series(
1090 cubes: iris.cube.CubeList,
1091 coords: list[iris.coords.Coord],
1092 ensemble_coord: str,
1093 filename: str,
1094 series_coordinate: str,
1095 title: str,
1096 vmin: float,
1097 vmax: float,
1098 **kwargs,
1099):
1100 """Plot and save a 1D line series in vertical.
1102 Parameters
1103 ----------
1104 cubes: CubeList
1105 1 dimensional Cube or CubeList of the data to plot on x-axis.
1106 coord: list[Coord]
1107 Coordinates to plot on the y-axis, one per cube.
1108 ensemble_coord: str
1109 Ensemble coordinate in the cube.
1110 filename: str
1111 Filename of the plot to write.
1112 series_coordinate: str
1113 Coordinate to use as vertical axis.
1114 title: str
1115 Plot title.
1116 vmin: float
1117 Minimum value for the x-axis.
1118 vmax: float
1119 Maximum value for the x-axis.
1120 """
1121 # plot the vertical pressure axis using log scale
1122 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1124 model_colors_map = _get_model_colors_map(cubes)
1126 # Check match-up across sequence coords gives consistent sizes
1127 _validate_cubes_coords(cubes, coords)
1129 for cube, coord in zip(cubes, coords, strict=True):
1130 label = None
1131 color = "black"
1132 if model_colors_map: 1132 ↛ 1133line 1132 didn't jump to line 1133 because the condition on line 1132 was never true
1133 label = cube.attributes.get("model_name")
1134 color = model_colors_map.get(label)
1136 for cube_slice in cube.slices_over(ensemble_coord):
1137 # If ensemble data given plot control member with (control)
1138 # unless single forecast.
1139 if cube_slice.coord(ensemble_coord).points == [0]:
1140 iplt.plot(
1141 cube_slice,
1142 coord,
1143 color=color,
1144 marker="o",
1145 ls="-",
1146 lw=3,
1147 label=f"{label} (control)"
1148 if len(cube.coord(ensemble_coord).points) > 1
1149 else label,
1150 )
1151 # If ensemble data given plot perturbed members with (perturbed).
1152 else:
1153 iplt.plot(
1154 cube_slice,
1155 coord,
1156 color=color,
1157 ls="-",
1158 lw=1.5,
1159 alpha=0.75,
1160 label=f"{label} (member)",
1161 )
1163 # Get the current axis
1164 ax = plt.gca()
1166 # Special handling for pressure level data.
1167 if series_coordinate == "pressure": 1167 ↛ 1189line 1167 didn't jump to line 1189 because the condition on line 1167 was always true
1168 # Invert y-axis and set to log scale.
1169 ax.invert_yaxis()
1170 ax.set_yscale("log")
1172 # Define y-ticks and labels for pressure log axis.
1173 y_tick_labels = [
1174 "1000",
1175 "850",
1176 "700",
1177 "500",
1178 "300",
1179 "200",
1180 "100",
1181 ]
1182 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1184 # Set y-axis limits and ticks.
1185 ax.set_ylim(1100, 100)
1187 # Test if series_coordinate is model level data. The UM data uses
1188 # model_level_number and lfric uses full_levels as coordinate.
1189 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1190 # Define y-ticks and labels for vertical axis.
1191 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1192 y_tick_labels = [str(int(i)) for i in y_ticks]
1193 ax.set_ylim(min(y_ticks), max(y_ticks))
1195 ax.set_yticks(y_ticks)
1196 ax.set_yticklabels(y_tick_labels)
1198 # Set x-axis limits.
1199 ax.set_xlim(vmin, vmax)
1200 # Mark y=0 if present in plot.
1201 if vmin < 0.0 and vmax > 0.0: 1201 ↛ 1202line 1201 didn't jump to line 1202 because the condition on line 1201 was never true
1202 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1204 # Add some labels and tweak the style.
1205 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1206 ax.set_xlabel(
1207 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1208 )
1209 ax.set_title(title, fontsize=16)
1210 ax.ticklabel_format(axis="x")
1211 ax.tick_params(axis="y")
1212 ax.tick_params(axis="both", labelsize=12)
1214 # Add gridlines
1215 ax.grid(linestyle="--", color="grey", linewidth=1)
1216 # Ientify unique labels for legend
1217 handles = list(
1218 {
1219 label: handle
1220 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1221 }.values()
1222 )
1223 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1225 # Save plot.
1226 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1227 logging.info("Saved line plot to %s", filename)
1228 plt.close(fig)
1231def _plot_and_save_scatter_plot(
1232 cube_x: iris.cube.Cube | iris.cube.CubeList,
1233 cube_y: iris.cube.Cube | iris.cube.CubeList,
1234 filename: str,
1235 title: str,
1236 one_to_one: bool,
1237 model_names: list[str] = None,
1238 **kwargs,
1239):
1240 """Plot and save a 2D scatter plot.
1242 Parameters
1243 ----------
1244 cube_x: Cube | CubeList
1245 1 dimensional Cube or CubeList of the data to plot on x-axis.
1246 cube_y: Cube | CubeList
1247 1 dimensional Cube or CubeList of the data to plot on y-axis.
1248 filename: str
1249 Filename of the plot to write.
1250 title: str
1251 Plot title.
1252 one_to_one: bool
1253 Whether a 1:1 line is plotted.
1254 """
1255 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1256 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1257 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1258 # over the pairs simultaneously.
1260 # Ensure cube_x and cube_y are iterable
1261 cube_x_iterable = iter_maybe(cube_x)
1262 cube_y_iterable = iter_maybe(cube_y)
1264 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1265 iplt.scatter(cube_x_iter, cube_y_iter)
1266 if one_to_one is True:
1267 plt.plot(
1268 [
1269 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1270 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1271 ],
1272 [
1273 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1274 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1275 ],
1276 "k",
1277 linestyle="--",
1278 )
1279 ax = plt.gca()
1281 # Add some labels and tweak the style.
1282 if model_names is None:
1283 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1284 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1285 else:
1286 # Add the model names, these should be order of base (x) and other (y).
1287 ax.set_xlabel(
1288 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1289 )
1290 ax.set_ylabel(
1291 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1292 )
1293 ax.set_title(title, fontsize=16)
1294 ax.ticklabel_format(axis="y", useOffset=False)
1295 ax.tick_params(axis="x", labelrotation=15)
1296 ax.tick_params(axis="both", labelsize=12)
1297 ax.autoscale()
1299 # Save plot.
1300 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1301 logging.info("Saved scatter plot to %s", filename)
1302 plt.close(fig)
1305def _plot_and_save_vector_plot(
1306 cube_u: iris.cube.Cube,
1307 cube_v: iris.cube.Cube,
1308 filename: str,
1309 title: str,
1310 method: Literal["contourf", "pcolormesh"],
1311 **kwargs,
1312):
1313 """Plot and save a 2D vector plot.
1315 Parameters
1316 ----------
1317 cube_u: Cube
1318 2 dimensional Cube of u component of the data.
1319 cube_v: Cube
1320 2 dimensional Cube of v component of the data.
1321 filename: str
1322 Filename of the plot to write.
1323 title: str
1324 Plot title.
1325 """
1326 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1328 # Create a cube containing the magnitude of the vector field.
1329 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1330 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1332 # Specify the color bar
1333 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1335 # Setup plot map projection, extent and coastlines and borderlines.
1336 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1338 if method == "contourf":
1339 # Filled contour plot of the field.
1340 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1341 elif method == "pcolormesh":
1342 try:
1343 vmin = min(levels)
1344 vmax = max(levels)
1345 except TypeError:
1346 vmin, vmax = None, None
1347 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1348 # if levels are defined.
1349 if norm is not None:
1350 vmin = None
1351 vmax = None
1352 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1353 else:
1354 raise ValueError(f"Unknown plotting method: {method}")
1356 # Check to see if transect, and if so, adjust y axis.
1357 if is_transect(cube_vec_mag):
1358 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1359 axes.invert_yaxis()
1360 axes.set_yscale("log")
1361 axes.set_ylim(1100, 100)
1362 # If both model_level_number and level_height exists, iplt can construct
1363 # plot as a function of height above orography (NOT sea level).
1364 elif {"model_level_number", "level_height"}.issubset(
1365 {coord.name() for coord in cube_vec_mag.coords()}
1366 ):
1367 axes.set_yscale("log")
1369 axes.set_title(
1370 f"{title}\n"
1371 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1372 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1373 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1374 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1375 fontsize=16,
1376 )
1378 else:
1379 # Add title.
1380 axes.set_title(title, fontsize=16)
1382 # Add watermark with min/max/mean. Currently not user togglable.
1383 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1384 axes.annotate(
1385 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}",
1386 xy=(0.05, -0.05),
1387 xycoords="axes fraction",
1388 xytext=(-5, 5),
1389 textcoords="offset points",
1390 ha="right",
1391 va="bottom",
1392 size=11,
1393 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1394 )
1396 # Add colour bar.
1397 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1398 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1399 # add ticks and tick_labels for every levels if less than 20 levels exist
1400 if levels is not None and len(levels) < 20:
1401 cbar.set_ticks(levels)
1402 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1404 # 30 barbs along the longest axis of the plot, or a barb per point for data
1405 # with less than 30 points.
1406 step = max(max(cube_u.shape) // 30, 1)
1407 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1409 # Save plot.
1410 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1411 logging.info("Saved vector plot to %s", filename)
1412 plt.close(fig)
1415def _plot_and_save_histogram_series(
1416 cubes: iris.cube.Cube | iris.cube.CubeList,
1417 filename: str,
1418 title: str,
1419 vmin: float,
1420 vmax: float,
1421 **kwargs,
1422):
1423 """Plot and save a histogram series.
1425 Parameters
1426 ----------
1427 cubes: Cube or CubeList
1428 2 dimensional Cube or CubeList of the data to plot as histogram.
1429 filename: str
1430 Filename of the plot to write.
1431 title: str
1432 Plot title.
1433 vmin: float
1434 minimum for colorbar
1435 vmax: float
1436 maximum for colorbar
1437 """
1438 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1439 ax = plt.gca()
1441 model_colors_map = _get_model_colors_map(cubes)
1443 # Set default that histograms will produce probability density function
1444 # at each bin (integral over range sums to 1).
1445 density = True
1447 for cube in iter_maybe(cubes):
1448 # Easier to check title (where var name originates)
1449 # than seeing if long names exist etc.
1450 # Exception case, where distribution better fits log scales/bins.
1451 if "surface_microphysical" in title:
1452 if "amount" in title: 1452 ↛ 1454line 1452 didn't jump to line 1454 because the condition on line 1452 was never true
1453 # Compute histogram following Klingaman et al. (2017): ASoP
1454 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1455 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1456 density = False
1457 else:
1458 bins = 10.0 ** (
1459 np.arange(-10, 27, 1) / 10.0
1460 ) # Suggestion from RMED toolbox.
1461 bins = np.insert(bins, 0, 0)
1462 ax.set_yscale("log")
1463 vmin = bins[1]
1464 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1465 ax.set_xscale("log")
1466 elif "lightning" in title:
1467 bins = [0, 1, 2, 3, 4, 5]
1468 else:
1469 bins = np.linspace(vmin, vmax, 51)
1470 logging.debug(
1471 "Plotting histogram with %s bins %s - %s.",
1472 np.size(bins),
1473 np.min(bins),
1474 np.max(bins),
1475 )
1477 # Reshape cube data into a single array to allow for a single histogram.
1478 # Otherwise we plot xdim histograms stacked.
1479 cube_data_1d = (cube.data).flatten()
1481 label = None
1482 color = "black"
1483 if model_colors_map: 1483 ↛ 1484line 1483 didn't jump to line 1484 because the condition on line 1483 was never true
1484 label = cube.attributes.get("model_name")
1485 color = model_colors_map[label]
1486 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1488 # Compute area under curve.
1489 if "surface_microphysical" in title and "amount" in title: 1489 ↛ 1490line 1489 didn't jump to line 1490 because the condition on line 1489 was never true
1490 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1491 x = x * bin_mean / x.sum()
1492 x = x[1:]
1493 y = y[1:]
1495 ax.plot(
1496 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1497 )
1499 # Add some labels and tweak the style.
1500 ax.set_title(title, fontsize=16)
1501 ax.set_xlabel(
1502 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1503 )
1504 ax.set_ylabel("Normalised probability density", fontsize=14)
1505 if "surface_microphysical" in title and "amount" in title: 1505 ↛ 1506line 1505 didn't jump to line 1506 because the condition on line 1505 was never true
1506 ax.set_ylabel(
1507 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1508 )
1509 try:
1510 ax.set_xlim(vmin, vmax)
1511 except ValueError:
1512 pass
1513 ax.tick_params(axis="both", labelsize=12)
1515 # Overlay grid-lines onto histogram plot.
1516 ax.grid(linestyle="--", color="grey", linewidth=1)
1517 if model_colors_map: 1517 ↛ 1518line 1517 didn't jump to line 1518 because the condition on line 1517 was never true
1518 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1520 # Save plot.
1521 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1522 logging.info("Saved histogram plot to %s", filename)
1523 plt.close(fig)
1526def _plot_and_save_postage_stamp_histogram_series(
1527 cube: iris.cube.Cube,
1528 filename: str,
1529 title: str,
1530 stamp_coordinate: str,
1531 vmin: float,
1532 vmax: float,
1533 **kwargs,
1534):
1535 """Plot and save postage (ensemble members) stamps for a histogram series.
1537 Parameters
1538 ----------
1539 cube: Cube
1540 2 dimensional Cube of the data to plot as histogram.
1541 filename: str
1542 Filename of the plot to write.
1543 title: str
1544 Plot title.
1545 stamp_coordinate: str
1546 Coordinate that becomes different plots.
1547 vmin: float
1548 minimum for pdf x-axis
1549 vmax: float
1550 maximum for pdf x-axis
1551 """
1552 # Use the smallest square grid that will fit the members.
1553 nmember = len(cube.coord(stamp_coordinate).points)
1554 grid_rows = int(math.sqrt(nmember))
1555 grid_size = math.ceil(nmember / grid_rows)
1557 fig = plt.figure(
1558 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1559 )
1560 # Make a subplot for each member.
1561 for member, subplot in zip(
1562 cube.slices_over(stamp_coordinate),
1563 range(1, grid_size * grid_rows + 1),
1564 strict=False,
1565 ):
1566 # Implicit interface is much easier here, due to needing to have the
1567 # cartopy GeoAxes generated.
1568 plt.subplot(grid_rows, grid_size, subplot)
1569 # Reshape cube data into a single array to allow for a single histogram.
1570 # Otherwise we plot xdim histograms stacked.
1571 member_data_1d = (member.data).flatten()
1572 plt.hist(member_data_1d, density=True, stacked=True)
1573 axes = plt.gca()
1574 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1575 axes.set_title(f"{mtitle}")
1576 axes.set_xlim(vmin, vmax)
1578 # Overall figure title.
1579 fig.suptitle(title, fontsize=16)
1581 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1582 logging.info("Saved histogram postage stamp plot to %s", filename)
1583 plt.close(fig)
1586def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1587 cube: iris.cube.Cube,
1588 filename: str,
1589 title: str,
1590 stamp_coordinate: str,
1591 vmin: float,
1592 vmax: float,
1593 **kwargs,
1594):
1595 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1596 ax.set_title(title, fontsize=16)
1597 ax.set_xlim(vmin, vmax)
1598 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1599 ax.set_ylabel("normalised probability density", fontsize=14)
1600 # Loop over all slices along the stamp_coordinate
1601 for member in cube.slices_over(stamp_coordinate):
1602 # Flatten the member data to 1D
1603 member_data_1d = member.data.flatten()
1604 # Plot the histogram using plt.hist
1605 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1606 plt.hist(
1607 member_data_1d,
1608 density=True,
1609 stacked=True,
1610 label=f"{mtitle}",
1611 )
1613 # Add a legend
1614 ax.legend(fontsize=16)
1616 # Save the figure to a file
1617 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1618 logging.info("Saved histogram postage stamp plot to %s", filename)
1620 # Close the figure
1621 plt.close(fig)
1624def _plot_and_save_scattermap_plot(
1625 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1626):
1627 """Plot and save a geographical scatter plot.
1629 Parameters
1630 ----------
1631 cube: Cube
1632 1 dimensional Cube of the data points with auxiliary latitude and
1633 longitude coordinates,
1634 filename: str
1635 Filename of the plot to write.
1636 title: str
1637 Plot title.
1638 projection: str
1639 Mapping projection to be used by cartopy.
1640 """
1641 # Setup plot details, size, resolution, etc.
1642 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1643 if projection is not None:
1644 # Apart from the default, the only projection we currently support is
1645 # a stereographic projection over the North Pole.
1646 if projection == "NP_Stereo":
1647 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1648 else:
1649 raise ValueError(f"Unknown projection: {projection}")
1650 else:
1651 axes = plt.axes(projection=ccrs.PlateCarree())
1653 # Scatter plot of the field. The marker size is chosen to give
1654 # symbols that decrease in size as the number of observations
1655 # increases, although the fraction of the figure covered by
1656 # symbols increases roughly as N^(1/2), disregarding overlaps,
1657 # and has been selected for the default figure size of (10, 10).
1658 # Should this be changed, the marker size should be adjusted in
1659 # proportion to the area of the figure.
1660 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1661 klon = None
1662 klat = None
1663 for kc in range(len(cube.aux_coords)):
1664 if cube.aux_coords[kc].standard_name == "latitude":
1665 klat = kc
1666 elif cube.aux_coords[kc].standard_name == "longitude":
1667 klon = kc
1668 scatter_map = iplt.scatter(
1669 cube.aux_coords[klon],
1670 cube.aux_coords[klat],
1671 c=cube.data[:],
1672 s=mrk_size,
1673 cmap="jet",
1674 edgecolors="k",
1675 )
1677 # Add coastlines and borderlines.
1678 try:
1679 axes.coastlines(resolution="10m")
1680 axes.add_feature(cfeature.BORDERS)
1681 except AttributeError:
1682 pass
1684 # Add title.
1685 axes.set_title(title, fontsize=16)
1687 # Add colour bar.
1688 cbar = fig.colorbar(scatter_map)
1689 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1691 # Save plot.
1692 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1693 logging.info("Saved geographical scatter plot to %s", filename)
1694 plt.close(fig)
1697def _plot_and_save_power_spectrum_series(
1698 cubes: iris.cube.Cube | iris.cube.CubeList,
1699 filename: str,
1700 title: str,
1701 **kwargs,
1702):
1703 """Plot and save a power spectrum series.
1705 Parameters
1706 ----------
1707 cubes: Cube or CubeList
1708 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1709 filename: str
1710 Filename of the plot to write.
1711 title: str
1712 Plot title.
1713 """
1714 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1715 ax = plt.gca()
1717 model_colors_map = _get_model_colors_map(cubes)
1719 for cube in iter_maybe(cubes):
1720 # Calculate power spectrum
1722 # Extract time coordinate and convert to datetime
1723 time_coord = cube.coord("time")
1724 time_points = time_coord.units.num2date(time_coord.points)
1726 # Choose one time point (e.g., the first one)
1727 target_time = time_points[0]
1729 # Bind target_time inside the lambda using a default argument
1730 time_constraint = iris.Constraint(
1731 time=lambda cell, target_time=target_time: cell.point == target_time
1732 )
1734 cube = cube.extract(time_constraint)
1736 if cube.ndim == 2:
1737 cube_3d = cube.data[np.newaxis, :, :]
1738 logging.debug("Adding in new axis for a 2 dimensional cube.")
1739 elif cube.ndim == 3: 1739 ↛ 1740line 1739 didn't jump to line 1740 because the condition on line 1739 was never true
1740 cube_3d = cube.data
1741 else:
1742 raise ValueError("Cube dimensions unsuitable for power spectra code")
1743 raise ValueError(
1744 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1745 )
1747 # Calculate spectra
1748 ps_array = _DCT_ps(cube_3d)
1750 ps_cube = iris.cube.Cube(
1751 ps_array,
1752 long_name="power_spectra",
1753 )
1755 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1757 # Create a frequency/wavelength array for coordinate
1758 ps_len = ps_cube.data.shape[1]
1759 freqs = np.arange(1, ps_len + 1)
1760 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1762 # Convert datetime to numeric time using original units
1763 numeric_time = time_coord.units.date2num(time_points)
1764 # Create a new DimCoord with numeric time
1765 new_time_coord = iris.coords.DimCoord(
1766 numeric_time, standard_name="time", units=time_coord.units
1767 )
1769 # Add time and frequency coordinate to spectra cube.
1770 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1771 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1773 # Extract data from the cube
1774 frequency = ps_cube.coord("frequency").points
1775 power_spectrum = ps_cube.data
1777 label = None
1778 color = "black"
1779 if model_colors_map: 1779 ↛ 1780line 1779 didn't jump to line 1780 because the condition on line 1779 was never true
1780 label = ps_cube.attributes.get("model_name")
1781 color = model_colors_map[label]
1782 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1784 # Add some labels and tweak the style.
1785 ax.set_title(title, fontsize=16)
1786 ax.set_xlabel("Wavenumber", fontsize=14)
1787 ax.set_ylabel("Power", fontsize=14)
1788 ax.tick_params(axis="both", labelsize=12)
1790 # Set log-log scale
1791 ax.set_xscale("log")
1792 ax.set_yscale("log")
1794 # Overlay grid-lines onto power spectrum plot.
1795 ax.grid(linestyle="--", color="grey", linewidth=1)
1796 if model_colors_map: 1796 ↛ 1797line 1796 didn't jump to line 1797 because the condition on line 1796 was never true
1797 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1799 # Save plot.
1800 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1801 logging.info("Saved power spectrum plot to %s", filename)
1802 plt.close(fig)
1805def _plot_and_save_postage_stamp_power_spectrum_series(
1806 cube: iris.cube.Cube,
1807 filename: str,
1808 title: str,
1809 stamp_coordinate: str,
1810 **kwargs,
1811):
1812 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1814 Parameters
1815 ----------
1816 cube: Cube
1817 2 dimensional Cube of the data to plot as power spectrum.
1818 filename: str
1819 Filename of the plot to write.
1820 title: str
1821 Plot title.
1822 stamp_coordinate: str
1823 Coordinate that becomes different plots.
1824 """
1825 # Use the smallest square grid that will fit the members.
1826 nmember = len(cube.coord(stamp_coordinate).points)
1827 grid_rows = int(math.sqrt(nmember))
1828 grid_size = math.ceil(nmember / grid_rows)
1830 fig = plt.figure(
1831 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1832 )
1834 # Make a subplot for each member.
1835 for member, subplot in zip(
1836 cube.slices_over(stamp_coordinate),
1837 range(1, grid_size * grid_rows + 1),
1838 strict=False,
1839 ):
1840 # Implicit interface is much easier here, due to needing to have the
1841 # cartopy GeoAxes generated.
1842 plt.subplot(grid_rows, grid_size, subplot)
1844 frequency = member.coord("frequency").points
1846 axes = plt.gca()
1847 axes.plot(frequency, member.data)
1848 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1850 # Overall figure title.
1851 fig.suptitle(title, fontsize=16)
1853 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1854 logging.info("Saved power spectra postage stamp plot to %s", filename)
1855 plt.close(fig)
1858def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1859 cube: iris.cube.Cube,
1860 filename: str,
1861 title: str,
1862 stamp_coordinate: str,
1863 **kwargs,
1864):
1865 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1866 ax.set_title(title, fontsize=16)
1867 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1868 ax.set_ylabel("Power", fontsize=14)
1869 # Loop over all slices along the stamp_coordinate
1870 for member in cube.slices_over(stamp_coordinate):
1871 frequency = member.coord("frequency").points
1872 ax.plot(
1873 frequency,
1874 member.data,
1875 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1876 )
1878 # Add a legend
1879 ax.legend(fontsize=16)
1881 # Save the figure to a file
1882 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1883 logging.info("Saved power spectra plot to %s", filename)
1885 # Close the figure
1886 plt.close(fig)
1889def _spatial_plot(
1890 method: Literal["contourf", "pcolormesh"],
1891 cube: iris.cube.Cube,
1892 filename: str | None,
1893 sequence_coordinate: str,
1894 stamp_coordinate: str,
1895 overlay_cube: iris.cube.Cube | None = None,
1896 contour_cube: iris.cube.Cube | None = None,
1897 **kwargs,
1898):
1899 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1901 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1902 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1903 is present then postage stamp plots will be produced.
1905 If an overlay_cube and/or contour_cube are specified, multiple variables can
1906 be overplotted on the same figure.
1908 Parameters
1909 ----------
1910 method: "contourf" | "pcolormesh"
1911 The plotting method to use.
1912 cube: Cube
1913 Iris cube of the data to plot. It should have two spatial dimensions,
1914 such as lat and lon, and may also have a another two dimension to be
1915 plotted sequentially and/or as postage stamp plots.
1916 filename: str | None
1917 Name of the plot to write, used as a prefix for plot sequences. If None
1918 uses the recipe name.
1919 sequence_coordinate: str
1920 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1921 This coordinate must exist in the cube.
1922 stamp_coordinate: str
1923 Coordinate about which to plot postage stamp plots. Defaults to
1924 ``"realization"``.
1925 overlay_cube: Cube | None, optional
1926 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1927 contour_cube: Cube | None, optional
1928 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1930 Raises
1931 ------
1932 ValueError
1933 If the cube doesn't have the right dimensions.
1934 TypeError
1935 If the cube isn't a single cube.
1936 """
1937 recipe_title = get_recipe_metadata().get("title", "Untitled")
1939 # Ensure we've got a single cube.
1940 cube = _check_single_cube(cube)
1942 # Check if there is a valid stamp coordinate in cube dimensions.
1943 if stamp_coordinate == "realization": 1943 ↛ 1948line 1943 didn't jump to line 1948 because the condition on line 1943 was always true
1944 stamp_coordinate = check_stamp_coordinate(cube)
1946 # Make postage stamp plots if stamp_coordinate exists and has more than a
1947 # single point.
1948 plotting_func = _plot_and_save_spatial_plot
1949 try:
1950 if cube.coord(stamp_coordinate).shape[0] > 1:
1951 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1952 except iris.exceptions.CoordinateNotFoundError:
1953 pass
1955 # Produce a geographical scatter plot if the data have a
1956 # dimension called observation or model_obs_error
1957 if any( 1957 ↛ 1961line 1957 didn't jump to line 1961 because the condition on line 1957 was never true
1958 crd.var_name == "station" or crd.var_name == "model_obs_error"
1959 for crd in cube.coords()
1960 ):
1961 plotting_func = _plot_and_save_scattermap_plot
1963 # Must have a sequence coordinate.
1964 try:
1965 cube.coord(sequence_coordinate)
1966 except iris.exceptions.CoordinateNotFoundError as err:
1967 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1969 # Create a plot for each value of the sequence coordinate.
1970 plot_index = []
1971 nplot = np.size(cube.coord(sequence_coordinate).points)
1973 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1974 # Set plot titles and filename
1975 seq_coord = cube_slice.coord(sequence_coordinate)
1976 plot_title, plot_filename = _set_title_and_filename(
1977 seq_coord, nplot, recipe_title, filename
1978 )
1980 # Extract sequence slice for overlay_cube and contour_cube if required.
1981 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1982 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1984 # Do the actual plotting.
1985 plotting_func(
1986 cube_slice,
1987 filename=plot_filename,
1988 stamp_coordinate=stamp_coordinate,
1989 title=plot_title,
1990 method=method,
1991 overlay_cube=overlay_slice,
1992 contour_cube=contour_slice,
1993 **kwargs,
1994 )
1995 plot_index.append(plot_filename)
1997 # Add list of plots to plot metadata.
1998 complete_plot_index = _append_to_plot_index(plot_index)
2000 # Make a page to display the plots.
2001 _make_plot_html_page(complete_plot_index)
2004def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2005 """Get colourmap for mask.
2007 If "mask_for_" appears anywhere in the name of a cube this function will be called
2008 regardless of the name of the variable to ensure a consistent plot.
2010 Parameters
2011 ----------
2012 cube: Cube
2013 Cube of variable for which the colorbar information is desired.
2014 axis: "x", "y", optional
2015 Select the levels for just this axis of a line plot. The min and max
2016 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2017 setting a universal range is not desirable (e.g. temperature), users
2018 can set ymin/ymax values to "auto" in the colorbar definitions file.
2019 Where no additional xmin/xmax or ymin/ymax values are provided, the
2020 axis bounds default to use the vmin/vmax values provided.
2022 Returns
2023 -------
2024 cmap:
2025 Matplotlib colormap.
2026 levels:
2027 List of levels to use for plotting. For continuous plots the min and max
2028 should be taken as the range.
2029 norm:
2030 BoundaryNorm information.
2031 """
2032 if "difference" not in cube.long_name:
2033 if axis:
2034 levels = [0, 1]
2035 # Complete settings based on levels.
2036 return None, levels, None
2037 else:
2038 # Define the levels and colors.
2039 levels = [0, 1, 2]
2040 colors = ["white", "dodgerblue"]
2041 # Create a custom color map.
2042 cmap = mcolors.ListedColormap(colors)
2043 # Normalize the levels.
2044 norm = mcolors.BoundaryNorm(levels, cmap.N)
2045 logging.debug("Colourmap for %s.", cube.long_name)
2046 return cmap, levels, norm
2047 else:
2048 if axis:
2049 levels = [-1, 1]
2050 return None, levels, None
2051 else:
2052 # Search for if mask difference, set to +/- 0.5 as values plotted <
2053 # not <=.
2054 levels = [-2, -0.5, 0.5, 2]
2055 colors = ["goldenrod", "white", "teal"]
2056 cmap = mcolors.ListedColormap(colors)
2057 norm = mcolors.BoundaryNorm(levels, cmap.N)
2058 logging.debug("Colourmap for %s.", cube.long_name)
2059 return cmap, levels, norm
2062def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2063 """Get a custom colorbar for a cube in the Beaufort Scale.
2065 Specific variable ranges can be separately set in user-supplied definition
2066 for x- or y-axis limits, or indicate where automated range preferred.
2068 Parameters
2069 ----------
2070 cube: Cube
2071 Cube of variable with Beaufort Scale in name.
2072 axis: "x", "y", optional
2073 Select the levels for just this axis of a line plot. The min and max
2074 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2075 setting a universal range is not desirable (e.g. temperature), users
2076 can set ymin/ymax values to "auto" in the colorbar definitions file.
2077 Where no additional xmin/xmax or ymin/ymax values are provided, the
2078 axis bounds default to use the vmin/vmax values provided.
2080 Returns
2081 -------
2082 cmap:
2083 Matplotlib colormap.
2084 levels:
2085 List of levels to use for plotting. For continuous plots the min and max
2086 should be taken as the range.
2087 norm:
2088 BoundaryNorm information.
2089 """
2090 if "difference" not in cube.long_name:
2091 if axis:
2092 levels = [0, 12]
2093 return None, levels, None
2094 else:
2095 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
2096 colors = [
2097 "black",
2098 (0, 0, 0.6),
2099 "blue",
2100 "cyan",
2101 "green",
2102 "yellow",
2103 (1, 0.5, 0),
2104 "red",
2105 "pink",
2106 "magenta",
2107 "purple",
2108 "maroon",
2109 "white",
2110 ]
2111 cmap = mcolors.ListedColormap(colors)
2112 norm = mcolors.BoundaryNorm(levels, cmap.N)
2113 logging.info("change colormap for Beaufort Scale colorbar.")
2114 return cmap, levels, norm
2115 else:
2116 if axis:
2117 levels = [-4, 4]
2118 return None, levels, None
2119 else:
2120 levels = [
2121 -3.5,
2122 -2.5,
2123 -1.5,
2124 -0.5,
2125 0.5,
2126 1.5,
2127 2.5,
2128 3.5,
2129 ]
2130 cmap = plt.get_cmap("bwr", 8)
2131 norm = mcolors.BoundaryNorm(levels, cmap.N)
2132 return cmap, levels, norm
2135def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
2136 """Return altered colourmap for temperature with change in units to Celsius.
2138 If "Celsius" appears anywhere in the name of a cube this function will be called.
2140 Parameters
2141 ----------
2142 cube: Cube
2143 Cube of variable for which the colorbar information is desired.
2144 cmap: Matplotlib colormap.
2145 levels: List
2146 List of levels to use for plotting. For continuous plots the min and max
2147 should be taken as the range.
2148 norm: BoundaryNorm.
2150 Returns
2151 -------
2152 cmap: Matplotlib colormap.
2153 levels: List
2154 List of levels to use for plotting. For continuous plots the min and max
2155 should be taken as the range.
2156 norm: BoundaryNorm.
2157 """
2158 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2159 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
2160 levels = np.array(levels)
2161 levels -= 273
2162 levels = levels.tolist()
2163 else:
2164 # Do nothing keep the existing colourbar attributes
2165 levels = levels
2166 cmap = cmap
2167 norm = norm
2168 return cmap, levels, norm
2171def _custom_colormap_probability(
2172 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
2173):
2174 """Get a custom colorbar for a probability cube.
2176 Specific variable ranges can be separately set in user-supplied definition
2177 for x- or y-axis limits, or indicate where automated range preferred.
2179 Parameters
2180 ----------
2181 cube: Cube
2182 Cube of variable with probability in name.
2183 axis: "x", "y", optional
2184 Select the levels for just this axis of a line plot. The min and max
2185 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2186 setting a universal range is not desirable (e.g. temperature), users
2187 can set ymin/ymax values to "auto" in the colorbar definitions file.
2188 Where no additional xmin/xmax or ymin/ymax values are provided, the
2189 axis bounds default to use the vmin/vmax values provided.
2191 Returns
2192 -------
2193 cmap:
2194 Matplotlib colormap.
2195 levels:
2196 List of levels to use for plotting. For continuous plots the min and max
2197 should be taken as the range.
2198 norm:
2199 BoundaryNorm information.
2200 """
2201 if axis:
2202 levels = [0, 1]
2203 return None, levels, None
2204 else:
2205 cmap = mcolors.ListedColormap(
2206 [
2207 "#FFFFFF",
2208 "#636363",
2209 "#e1dada",
2210 "#B5CAFF",
2211 "#8FB3FF",
2212 "#7F97FF",
2213 "#ABCF63",
2214 "#E8F59E",
2215 "#FFFA14",
2216 "#FFD121",
2217 "#FFA30A",
2218 ]
2219 )
2220 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]
2221 norm = mcolors.BoundaryNorm(levels, cmap.N)
2222 return cmap, levels, norm
2225def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
2226 """Return a custom colourmap for the current recipe."""
2227 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2228 if (
2229 any("surface_microphysical" in name for name in varnames)
2230 and "difference" not in cube.long_name
2231 and "mask" not in cube.long_name
2232 ):
2233 # Define the levels and colors
2234 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
2235 colors = [
2236 "w",
2237 (0, 0, 0.6),
2238 "b",
2239 "c",
2240 "g",
2241 "y",
2242 (1, 0.5, 0),
2243 "r",
2244 "pink",
2245 "m",
2246 "purple",
2247 "maroon",
2248 "gray",
2249 ]
2250 # Create a custom colormap
2251 cmap = mcolors.ListedColormap(colors)
2252 # Normalize the levels
2253 norm = mcolors.BoundaryNorm(levels, cmap.N)
2254 logging.info("change colormap for surface_microphysical variable colorbar.")
2255 else:
2256 # do nothing and keep existing colorbar attributes
2257 cmap = cmap
2258 levels = levels
2259 norm = norm
2260 return cmap, levels, norm
2263def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
2264 """Return custom colourmap for aviation colour state.
2266 If "aviation_colour_state" appears anywhere in the name of a cube
2267 this function will be called.
2269 Parameters
2270 ----------
2271 cube: Cube
2272 Cube of variable for which the colorbar information is desired.
2274 Returns
2275 -------
2276 cmap: Matplotlib colormap.
2277 levels: List
2278 List of levels to use for plotting. For continuous plots the min and max
2279 should be taken as the range.
2280 norm: BoundaryNorm.
2281 """
2282 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
2283 colors = [
2284 "#87ceeb",
2285 "#ffffff",
2286 "#8ced69",
2287 "#ffff00",
2288 "#ffd700",
2289 "#ffa500",
2290 "#fe3620",
2291 ]
2292 # Create a custom colormap
2293 cmap = mcolors.ListedColormap(colors)
2294 # Normalise the levels
2295 norm = mcolors.BoundaryNorm(levels, cmap.N)
2296 return cmap, levels, norm
2299def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
2300 """Return a custom colourmap for the current recipe."""
2301 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2302 if (
2303 any("visibility_in_air" in name for name in varnames)
2304 and "difference" not in cube.long_name
2305 and "mask" not in cube.long_name
2306 ):
2307 # Define the levels and colors (in km)
2308 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]
2309 norm = mcolors.BoundaryNorm(levels, cmap.N)
2310 colours = [
2311 "#8f00d6",
2312 "#d10000",
2313 "#ff9700",
2314 "#ffff00",
2315 "#00007f",
2316 "#6c9ccd",
2317 "#aae8ff",
2318 "#37a648",
2319 "#8edc64",
2320 "#c5ffc5",
2321 "#dcdcdc",
2322 "#ffffff",
2323 ]
2324 # Create a custom colormap
2325 cmap = mcolors.ListedColormap(colours)
2326 # Normalize the levels
2327 norm = mcolors.BoundaryNorm(levels, cmap.N)
2328 logging.info("change colormap for visibility_in_air variable colorbar.")
2329 else:
2330 # do nothing and keep existing colorbar attributes
2331 cmap = cmap
2332 levels = levels
2333 norm = norm
2334 return cmap, levels, norm
2337def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
2338 """Return number of models based on cube attributes."""
2339 model_names = list(
2340 filter(
2341 lambda x: x is not None,
2342 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
2343 )
2344 )
2345 if not model_names:
2346 logging.debug("Missing model names. Will assume single model.")
2347 return 1
2348 else:
2349 return len(model_names)
2352def _validate_cube_shape(
2353 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2354) -> None:
2355 """Check all cubes have a model name."""
2356 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2356 ↛ 2357line 2356 didn't jump to line 2357 because the condition on line 2356 was never true
2357 raise ValueError(
2358 f"The number of model names ({num_models}) should equal the number "
2359 f"of cubes ({len(cube)})."
2360 )
2363def _validate_cubes_coords(
2364 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2365) -> None:
2366 """Check same number of cubes as sequence coordinate for zip functions."""
2367 if len(cubes) != len(coords): 2367 ↛ 2368line 2367 didn't jump to line 2368 because the condition on line 2367 was never true
2368 raise ValueError(
2369 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2370 f"of sequence coordinates ({len(coords)})."
2371 f"Check that number of time entries in input data are consistent if "
2372 f"performing time-averaging steps prior to plotting outputs."
2373 )
2376####################
2377# Public functions #
2378####################
2381def spatial_contour_plot(
2382 cube: iris.cube.Cube,
2383 filename: str = None,
2384 sequence_coordinate: str = "time",
2385 stamp_coordinate: str = "realization",
2386 **kwargs,
2387) -> iris.cube.Cube:
2388 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2390 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2391 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2392 is present then postage stamp plots will be produced.
2394 Parameters
2395 ----------
2396 cube: Cube
2397 Iris cube of the data to plot. It should have two spatial dimensions,
2398 such as lat and lon, and may also have a another two dimension to be
2399 plotted sequentially and/or as postage stamp plots.
2400 filename: str, optional
2401 Name of the plot to write, used as a prefix for plot sequences. Defaults
2402 to the recipe name.
2403 sequence_coordinate: str, optional
2404 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2405 This coordinate must exist in the cube.
2406 stamp_coordinate: str, optional
2407 Coordinate about which to plot postage stamp plots. Defaults to
2408 ``"realization"``.
2410 Returns
2411 -------
2412 Cube
2413 The original cube (so further operations can be applied).
2415 Raises
2416 ------
2417 ValueError
2418 If the cube doesn't have the right dimensions.
2419 TypeError
2420 If the cube isn't a single cube.
2421 """
2422 _spatial_plot(
2423 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2424 )
2425 return cube
2428def spatial_pcolormesh_plot(
2429 cube: iris.cube.Cube,
2430 filename: str = None,
2431 sequence_coordinate: str = "time",
2432 stamp_coordinate: str = "realization",
2433 **kwargs,
2434) -> iris.cube.Cube:
2435 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2437 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2438 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2439 is present then postage stamp plots will be produced.
2441 This function is significantly faster than ``spatial_contour_plot``,
2442 especially at high resolutions, and should be preferred unless contiguous
2443 contour areas are important.
2445 Parameters
2446 ----------
2447 cube: Cube
2448 Iris cube of the data to plot. It should have two spatial dimensions,
2449 such as lat and lon, and may also have a another two dimension to be
2450 plotted sequentially and/or as postage stamp plots.
2451 filename: str, optional
2452 Name of the plot to write, used as a prefix for plot sequences. Defaults
2453 to the recipe name.
2454 sequence_coordinate: str, optional
2455 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2456 This coordinate must exist in the cube.
2457 stamp_coordinate: str, optional
2458 Coordinate about which to plot postage stamp plots. Defaults to
2459 ``"realization"``.
2461 Returns
2462 -------
2463 Cube
2464 The original cube (so further operations can be applied).
2466 Raises
2467 ------
2468 ValueError
2469 If the cube doesn't have the right dimensions.
2470 TypeError
2471 If the cube isn't a single cube.
2472 """
2473 _spatial_plot(
2474 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2475 )
2476 return cube
2479def spatial_multi_pcolormesh_plot(
2480 cube: iris.cube.Cube,
2481 overlay_cube: iris.cube.Cube,
2482 contour_cube: iris.cube.Cube,
2483 filename: str = None,
2484 sequence_coordinate: str = "time",
2485 stamp_coordinate: str = "realization",
2486 **kwargs,
2487) -> iris.cube.Cube:
2488 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2490 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2491 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2492 is present then postage stamp plots will be produced.
2494 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2496 If specified, contours of a contour_cube can be overplotted on top of those.
2498 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2500 This function is significantly faster than ``spatial_contour_plot``,
2501 especially at high resolutions, and should be preferred unless contiguous
2502 contour areas are important.
2504 Parameters
2505 ----------
2506 cube: Cube
2507 Iris cube of the data to plot. It should have two spatial dimensions,
2508 such as lat and lon, and may also have a another two dimension to be
2509 plotted sequentially and/or as postage stamp plots.
2510 overlay_cube: Cube
2511 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2512 such as lat and lon, and may also have a another two dimension to be
2513 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.
2514 contour_cube: Cube
2515 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,
2516 such as lat and lon, and may also have a another two dimension to be
2517 plotted sequentially and/or as postage stamp plots.
2518 filename: str, optional
2519 Name of the plot to write, used as a prefix for plot sequences. Defaults
2520 to the recipe name.
2521 sequence_coordinate: str, optional
2522 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2523 This coordinate must exist in the cube.
2524 stamp_coordinate: str, optional
2525 Coordinate about which to plot postage stamp plots. Defaults to
2526 ``"realization"``.
2528 Returns
2529 -------
2530 Cube
2531 The original cube (so further operations can be applied).
2533 Raises
2534 ------
2535 ValueError
2536 If the cube doesn't have the right dimensions.
2537 TypeError
2538 If the cube isn't a single cube.
2539 """
2540 _spatial_plot(
2541 "pcolormesh",
2542 cube,
2543 filename,
2544 sequence_coordinate,
2545 stamp_coordinate,
2546 overlay_cube=overlay_cube,
2547 contour_cube=contour_cube,
2548 )
2549 return cube, overlay_cube, contour_cube
2552# TODO: Expand function to handle ensemble data.
2553# line_coordinate: str, optional
2554# Coordinate about which to plot multiple lines. Defaults to
2555# ``"realization"``.
2556def plot_line_series(
2557 cube: iris.cube.Cube | iris.cube.CubeList,
2558 filename: str = None,
2559 series_coordinate: str = "time",
2560 # line_coordinate: str = "realization",
2561 **kwargs,
2562) -> iris.cube.Cube | iris.cube.CubeList:
2563 """Plot a line plot for the specified coordinate.
2565 The Cube or CubeList must be 1D.
2567 Parameters
2568 ----------
2569 iris.cube | iris.cube.CubeList
2570 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2571 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2572 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2573 filename: str, optional
2574 Name of the plot to write, used as a prefix for plot sequences. Defaults
2575 to the recipe name.
2576 series_coordinate: str, optional
2577 Coordinate about which to make a series. Defaults to ``"time"``. This
2578 coordinate must exist in the cube.
2580 Returns
2581 -------
2582 iris.cube.Cube | iris.cube.CubeList
2583 The original Cube or CubeList (so further operations can be applied).
2584 plotted data.
2586 Raises
2587 ------
2588 ValueError
2589 If the cubes don't have the right dimensions.
2590 TypeError
2591 If the cube isn't a Cube or CubeList.
2592 """
2593 # Ensure we have a name for the plot file.
2594 recipe_title = get_recipe_metadata().get("title", "Untitled")
2596 num_models = _get_num_models(cube)
2598 _validate_cube_shape(cube, num_models)
2600 # Iterate over all cubes and extract coordinate to plot.
2601 cubes = iter_maybe(cube)
2602 coords = []
2603 for cube in cubes:
2604 try:
2605 coords.append(cube.coord(series_coordinate))
2606 except iris.exceptions.CoordinateNotFoundError as err:
2607 raise ValueError(
2608 f"Cube must have a {series_coordinate} coordinate."
2609 ) from err
2610 if cube.ndim > 2 or not cube.coords("realization"):
2611 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2613 # Format the title and filename using plotted series coordinate
2614 nplot = 1
2615 seq_coord = coords[0]
2616 plot_title, plot_filename = _set_title_and_filename(
2617 seq_coord, nplot, recipe_title, filename
2618 )
2620 # Do the actual plotting.
2621 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, plot_title)
2623 # Add list of plots to plot metadata.
2624 plot_index = _append_to_plot_index([plot_filename])
2626 # Make a page to display the plots.
2627 _make_plot_html_page(plot_index)
2629 return cube
2632def plot_vertical_line_series(
2633 cubes: iris.cube.Cube | iris.cube.CubeList,
2634 filename: str = None,
2635 series_coordinate: str = "model_level_number",
2636 sequence_coordinate: str = "time",
2637 # line_coordinate: str = "realization",
2638 **kwargs,
2639) -> iris.cube.Cube | iris.cube.CubeList:
2640 """Plot a line plot against a type of vertical coordinate.
2642 The Cube or CubeList must be 1D.
2644 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2645 then a sequence of plots will be produced.
2647 Parameters
2648 ----------
2649 iris.cube | iris.cube.CubeList
2650 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2651 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2652 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2653 filename: str, optional
2654 Name of the plot to write, used as a prefix for plot sequences. Defaults
2655 to the recipe name.
2656 series_coordinate: str, optional
2657 Coordinate to plot on the y-axis. Can be ``pressure`` or
2658 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2659 for LFRic. Defaults to ``model_level_number``.
2660 This coordinate must exist in the cube.
2661 sequence_coordinate: str, optional
2662 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2663 This coordinate must exist in the cube.
2665 Returns
2666 -------
2667 iris.cube.Cube | iris.cube.CubeList
2668 The original Cube or CubeList (so further operations can be applied).
2669 Plotted data.
2671 Raises
2672 ------
2673 ValueError
2674 If the cubes doesn't have the right dimensions.
2675 TypeError
2676 If the cube isn't a Cube or CubeList.
2677 """
2678 # Ensure we have a name for the plot file.
2679 recipe_title = get_recipe_metadata().get("title", "Untitled")
2681 cubes = iter_maybe(cubes)
2682 # Initialise empty list to hold all data from all cubes in a CubeList
2683 all_data = []
2685 # Store min/max ranges for x range.
2686 x_levels = []
2688 num_models = _get_num_models(cubes)
2690 _validate_cube_shape(cubes, num_models)
2692 # Iterate over all cubes in cube or CubeList and plot.
2693 coords = []
2694 for cube in cubes:
2695 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2696 try:
2697 coords.append(cube.coord(series_coordinate))
2698 except iris.exceptions.CoordinateNotFoundError as err:
2699 raise ValueError(
2700 f"Cube must have a {series_coordinate} coordinate."
2701 ) from err
2703 try:
2704 if cube.ndim > 1 or not cube.coords("realization"): 2704 ↛ 2712line 2704 didn't jump to line 2712 because the condition on line 2704 was always true
2705 cube.coord(sequence_coordinate)
2706 except iris.exceptions.CoordinateNotFoundError as err:
2707 raise ValueError(
2708 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2709 ) from err
2711 # Get minimum and maximum from levels information.
2712 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2713 if levels is not None: 2713 ↛ 2717line 2713 didn't jump to line 2717 because the condition on line 2713 was always true
2714 x_levels.append(min(levels))
2715 x_levels.append(max(levels))
2716 else:
2717 all_data.append(cube.data)
2719 if len(x_levels) == 0: 2719 ↛ 2721line 2719 didn't jump to line 2721 because the condition on line 2719 was never true
2720 # Combine all data into a single NumPy array
2721 combined_data = np.concatenate(all_data)
2723 # Set the lower and upper limit for the x-axis to ensure all plots have
2724 # same range. This needs to read the whole cube over the range of the
2725 # sequence and if applicable postage stamp coordinate.
2726 vmin = np.floor(combined_data.min())
2727 vmax = np.ceil(combined_data.max())
2728 else:
2729 vmin = min(x_levels)
2730 vmax = max(x_levels)
2732 # Matching the slices (matching by seq coord point; it may happen that
2733 # evaluated models do not cover the same seq coord range, hence matching
2734 # necessary)
2735 def filter_cube_iterables(cube_iterables) -> bool:
2736 return len(cube_iterables) == len(coords)
2738 cube_iterables = filter(
2739 filter_cube_iterables,
2740 (
2741 iris.cube.CubeList(
2742 s
2743 for s in itertools.chain.from_iterable(
2744 cb.slices_over(sequence_coordinate) for cb in cubes
2745 )
2746 if s.coord(sequence_coordinate).points[0] == point
2747 )
2748 for point in sorted(
2749 set(
2750 itertools.chain.from_iterable(
2751 cb.coord(sequence_coordinate).points for cb in cubes
2752 )
2753 )
2754 )
2755 ),
2756 )
2758 # Create a plot for each value of the sequence coordinate.
2759 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2760 # similar sequence values. Passing a CubeList into the internal plotting function
2761 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2762 # or an iris.cube.CubeList.
2763 plot_index = []
2764 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2765 for cubes_slice in cube_iterables:
2766 # Format the coordinate value in a unit appropriate way.
2767 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2768 plot_title, plot_filename = _set_title_and_filename(
2769 seq_coord, nplot, recipe_title, filename
2770 )
2772 # Do the actual plotting.
2773 _plot_and_save_vertical_line_series(
2774 cubes_slice,
2775 coords,
2776 "realization",
2777 plot_filename,
2778 series_coordinate,
2779 title=plot_title,
2780 vmin=vmin,
2781 vmax=vmax,
2782 )
2783 plot_index.append(plot_filename)
2785 # Add list of plots to plot metadata.
2786 complete_plot_index = _append_to_plot_index(plot_index)
2788 # Make a page to display the plots.
2789 _make_plot_html_page(complete_plot_index)
2791 return cubes
2794def qq_plot(
2795 cubes: iris.cube.CubeList,
2796 coordinates: list[str],
2797 percentiles: list[float],
2798 model_names: list[str],
2799 filename: str = None,
2800 one_to_one: bool = True,
2801 **kwargs,
2802) -> iris.cube.CubeList:
2803 """Plot a Quantile-Quantile plot between two models for common time points.
2805 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2806 collapsed within the operator over all specified coordinates such as
2807 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2808 ensemble members to ensure a 1D cube (array).
2810 Parameters
2811 ----------
2812 cubes: iris.cube.CubeList
2813 Two cubes of the same variable with different models.
2814 coordinate: list[str]
2815 The list of coordinates to collapse over. This list should be
2816 every coordinate within the cube to result in a 1D cube around
2817 the percentile coordinate.
2818 percent: list[float]
2819 A list of percentiles to appear in the plot.
2820 model_names: list[str]
2821 A list of model names to appear on the axis of the plot.
2822 filename: str, optional
2823 Filename of the plot to write.
2824 one_to_one: bool, optional
2825 If True a 1:1 line is plotted; if False it is not. Default is True.
2827 Raises
2828 ------
2829 ValueError
2830 When the cubes are not compatible.
2832 Notes
2833 -----
2834 The quantile-quantile plot is a variant on the scatter plot representing
2835 two datasets by their quantiles (percentiles) for common time points.
2836 This plot does not use a theoretical distribution to compare against, but
2837 compares percentiles of two datasets. This plot does
2838 not use all raw data points, but plots the selected percentiles (quantiles) of
2839 each variable instead for the two datasets, thereby normalising the data for a
2840 direct comparison between the selected percentiles of the two dataset distributions.
2842 Quantile-quantile plots are valuable for comparing against
2843 observations and other models. Identical percentiles between the variables
2844 will lie on the one-to-one line implying the values correspond well to each
2845 other. Where there is a deviation from the one-to-one line a range of
2846 possibilities exist depending on how and where the data is shifted (e.g.,
2847 Wilks 2011 [Wilks2011]_).
2849 For distributions above the one-to-one line the distribution is left-skewed;
2850 below is right-skewed. A distinct break implies a bimodal distribution, and
2851 closer values/values further apart at the tails imply poor representation of
2852 the extremes.
2854 References
2855 ----------
2856 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2857 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2858 """
2859 # Check cubes using same functionality as the difference operator.
2860 if len(cubes) != 2:
2861 raise ValueError("cubes should contain exactly 2 cubes.")
2862 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2863 other: Cube = cubes.extract_cube(
2864 iris.Constraint(
2865 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2866 )
2867 )
2869 # Get spatial coord names.
2870 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2871 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2873 # Ensure cubes to compare are on common differencing grid.
2874 # This is triggered if either
2875 # i) latitude and longitude shapes are not the same. Note grid points
2876 # are not compared directly as these can differ through rounding
2877 # errors.
2878 # ii) or variables are known to often sit on different grid staggering
2879 # in different models (e.g. cell center vs cell edge), as is the case
2880 # for UM and LFRic comparisons.
2881 # In future greater choice of regridding method might be applied depending
2882 # on variable type. Linear regridding can in general be appropriate for smooth
2883 # variables. Care should be taken with interpretation of differences
2884 # given this dependency on regridding.
2885 if (
2886 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2887 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2888 ) or (
2889 base.long_name
2890 in [
2891 "eastward_wind_at_10m",
2892 "northward_wind_at_10m",
2893 "northward_wind_at_cell_centres",
2894 "eastward_wind_at_cell_centres",
2895 "zonal_wind_at_pressure_levels",
2896 "meridional_wind_at_pressure_levels",
2897 "potential_vorticity_at_pressure_levels",
2898 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2899 ]
2900 ):
2901 logging.debug(
2902 "Linear regridding base cube to other grid to compute differences"
2903 )
2904 base = regrid_onto_cube(base, other, method="Linear")
2906 # Extract just common time points.
2907 base, other = _extract_common_time_points(base, other)
2909 # Equalise attributes so we can merge.
2910 fully_equalise_attributes([base, other])
2911 logging.debug("Base: %s\nOther: %s", base, other)
2913 # Collapse cubes.
2914 base = collapse(
2915 base,
2916 coordinate=coordinates,
2917 method="PERCENTILE",
2918 additional_percent=percentiles,
2919 )
2920 other = collapse(
2921 other,
2922 coordinate=coordinates,
2923 method="PERCENTILE",
2924 additional_percent=percentiles,
2925 )
2927 # Ensure we have a name for the plot file.
2928 recipe_title = get_recipe_metadata().get("title", "Untitled")
2929 title = f"{recipe_title}"
2931 if filename is None:
2932 filename = slugify(recipe_title)
2934 # Add file extension.
2935 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2937 # Do the actual plotting on a scatter plot
2938 _plot_and_save_scatter_plot(
2939 base, other, plot_filename, title, one_to_one, model_names
2940 )
2942 # Add list of plots to plot metadata.
2943 plot_index = _append_to_plot_index([plot_filename])
2945 # Make a page to display the plots.
2946 _make_plot_html_page(plot_index)
2948 return iris.cube.CubeList([base, other])
2951def scatter_plot(
2952 cube_x: iris.cube.Cube | iris.cube.CubeList,
2953 cube_y: iris.cube.Cube | iris.cube.CubeList,
2954 filename: str = None,
2955 one_to_one: bool = True,
2956 **kwargs,
2957) -> iris.cube.CubeList:
2958 """Plot a scatter plot between two variables.
2960 Both cubes must be 1D.
2962 Parameters
2963 ----------
2964 cube_x: Cube | CubeList
2965 1 dimensional Cube of the data to plot on y-axis.
2966 cube_y: Cube | CubeList
2967 1 dimensional Cube of the data to plot on x-axis.
2968 filename: str, optional
2969 Filename of the plot to write.
2970 one_to_one: bool, optional
2971 If True a 1:1 line is plotted; if False it is not. Default is True.
2973 Returns
2974 -------
2975 cubes: CubeList
2976 CubeList of the original x and y cubes for further processing.
2978 Raises
2979 ------
2980 ValueError
2981 If the cube doesn't have the right dimensions and cubes not the same
2982 size.
2983 TypeError
2984 If the cube isn't a single cube.
2986 Notes
2987 -----
2988 Scatter plots are used for determining if there is a relationship between
2989 two variables. Positive relations have a slope going from bottom left to top
2990 right; Negative relations have a slope going from top left to bottom right.
2991 """
2992 # Iterate over all cubes in cube or CubeList and plot.
2993 for cube_iter in iter_maybe(cube_x):
2994 # Check cubes are correct shape.
2995 cube_iter = _check_single_cube(cube_iter)
2996 if cube_iter.ndim > 1:
2997 raise ValueError("cube_x must be 1D.")
2999 # Iterate over all cubes in cube or CubeList and plot.
3000 for cube_iter in iter_maybe(cube_y):
3001 # Check cubes are correct shape.
3002 cube_iter = _check_single_cube(cube_iter)
3003 if cube_iter.ndim > 1:
3004 raise ValueError("cube_y must be 1D.")
3006 # Ensure we have a name for the plot file.
3007 recipe_title = get_recipe_metadata().get("title", "Untitled")
3008 title = f"{recipe_title}"
3010 if filename is None:
3011 filename = slugify(recipe_title)
3013 # Add file extension.
3014 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3016 # Do the actual plotting.
3017 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3019 # Add list of plots to plot metadata.
3020 plot_index = _append_to_plot_index([plot_filename])
3022 # Make a page to display the plots.
3023 _make_plot_html_page(plot_index)
3025 return iris.cube.CubeList([cube_x, cube_y])
3028def vector_plot(
3029 cube_u: iris.cube.Cube,
3030 cube_v: iris.cube.Cube,
3031 filename: str = None,
3032 sequence_coordinate: str = "time",
3033 **kwargs,
3034) -> iris.cube.CubeList:
3035 """Plot a vector plot based on the input u and v components."""
3036 recipe_title = get_recipe_metadata().get("title", "Untitled")
3038 # Cubes must have a matching sequence coordinate.
3039 try:
3040 # Check that the u and v cubes have the same sequence coordinate.
3041 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3041 ↛ anywhereline 3041 didn't jump anywhere: it always raised an exception.
3042 raise ValueError("Coordinates do not match.")
3043 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3044 raise ValueError(
3045 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3046 ) from err
3048 # Create a plot for each value of the sequence coordinate.
3049 plot_index = []
3050 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3051 for cube_u_slice, cube_v_slice in zip(
3052 cube_u.slices_over(sequence_coordinate),
3053 cube_v.slices_over(sequence_coordinate),
3054 strict=True,
3055 ):
3056 # Format the coordinate value in a unit appropriate way.
3057 seq_coord = cube_u_slice.coord(sequence_coordinate)
3058 plot_title, plot_filename = _set_title_and_filename(
3059 seq_coord, nplot, recipe_title, filename
3060 )
3062 # Do the actual plotting.
3063 _plot_and_save_vector_plot(
3064 cube_u_slice,
3065 cube_v_slice,
3066 filename=plot_filename,
3067 title=plot_title,
3068 method="contourf",
3069 )
3070 plot_index.append(plot_filename)
3072 # Add list of plots to plot metadata.
3073 complete_plot_index = _append_to_plot_index(plot_index)
3075 # Make a page to display the plots.
3076 _make_plot_html_page(complete_plot_index)
3078 return iris.cube.CubeList([cube_u, cube_v])
3081def plot_histogram_series(
3082 cubes: iris.cube.Cube | iris.cube.CubeList,
3083 filename: str = None,
3084 sequence_coordinate: str = "time",
3085 stamp_coordinate: str = "realization",
3086 single_plot: bool = False,
3087 **kwargs,
3088) -> iris.cube.Cube | iris.cube.CubeList:
3089 """Plot a histogram plot for each vertical level provided.
3091 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3092 is present then a sequence of plots will be produced using the time slider
3093 functionality to scroll through histograms against time. If a
3094 stamp_coordinate is present then postage stamp plots will be produced. If
3095 stamp_coordinate and single_plot is True, all postage stamp plots will be
3096 plotted in a single plot instead of separate postage stamp plots.
3098 Parameters
3099 ----------
3100 cubes: Cube | iris.cube.CubeList
3101 Iris cube or CubeList of the data to plot. It should have a single dimension other
3102 than the stamp coordinate.
3103 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3104 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3105 filename: str, optional
3106 Name of the plot to write, used as a prefix for plot sequences. Defaults
3107 to the recipe name.
3108 sequence_coordinate: str, optional
3109 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3110 This coordinate must exist in the cube and will be used for the time
3111 slider.
3112 stamp_coordinate: str, optional
3113 Coordinate about which to plot postage stamp plots. Defaults to
3114 ``"realization"``.
3115 single_plot: bool, optional
3116 If True, all postage stamp plots will be plotted in a single plot. If
3117 False, each postage stamp plot will be plotted separately. Is only valid
3118 if stamp_coordinate exists and has more than a single point.
3120 Returns
3121 -------
3122 iris.cube.Cube | iris.cube.CubeList
3123 The original Cube or CubeList (so further operations can be applied).
3124 Plotted data.
3126 Raises
3127 ------
3128 ValueError
3129 If the cube doesn't have the right dimensions.
3130 TypeError
3131 If the cube isn't a Cube or CubeList.
3132 """
3133 recipe_title = get_recipe_metadata().get("title", "Untitled")
3135 cubes = iter_maybe(cubes)
3137 # Internal plotting function.
3138 plotting_func = _plot_and_save_histogram_series
3140 num_models = _get_num_models(cubes)
3142 _validate_cube_shape(cubes, num_models)
3144 # If several histograms are plotted with time as sequence_coordinate for the
3145 # time slider option.
3146 for cube in cubes:
3147 try:
3148 cube.coord(sequence_coordinate)
3149 except iris.exceptions.CoordinateNotFoundError as err:
3150 raise ValueError(
3151 f"Cube must have a {sequence_coordinate} coordinate."
3152 ) from err
3154 # Get minimum and maximum from levels information.
3155 levels = None
3156 for cube in cubes: 3156 ↛ 3172line 3156 didn't jump to line 3172 because the loop on line 3156 didn't complete
3157 # First check if user-specified "auto" range variable.
3158 # This maintains the value of levels as None, so proceed.
3159 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3160 if levels is None:
3161 break
3162 # If levels is changed, recheck to use the vmin,vmax or
3163 # levels-based ranges for histogram plots.
3164 _, levels, _ = _colorbar_map_levels(cube)
3165 logging.debug("levels: %s", levels)
3166 if levels is not None: 3166 ↛ 3156line 3166 didn't jump to line 3156 because the condition on line 3166 was always true
3167 vmin = min(levels)
3168 vmax = max(levels)
3169 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
3170 break
3172 if levels is None:
3173 vmin = min(cb.data.min() for cb in cubes)
3174 vmax = max(cb.data.max() for cb in cubes)
3176 # Make postage stamp plots if stamp_coordinate exists and has more than a
3177 # single point. If single_plot is True:
3178 # -- all postage stamp plots will be plotted in a single plot instead of
3179 # separate postage stamp plots.
3180 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3181 # produced per single model only
3182 if num_models == 1: 3182 ↛ 3195line 3182 didn't jump to line 3195 because the condition on line 3182 was always true
3183 if ( 3183 ↛ 3187line 3183 didn't jump to line 3187 because the condition on line 3183 was never true
3184 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3185 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3186 ):
3187 if single_plot:
3188 plotting_func = (
3189 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3190 )
3191 else:
3192 plotting_func = _plot_and_save_postage_stamp_histogram_series
3193 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3194 else:
3195 all_points = sorted(
3196 set(
3197 itertools.chain.from_iterable(
3198 cb.coord(sequence_coordinate).points for cb in cubes
3199 )
3200 )
3201 )
3202 all_slices = list(
3203 itertools.chain.from_iterable(
3204 cb.slices_over(sequence_coordinate) for cb in cubes
3205 )
3206 )
3207 # Matched slices (matched by seq coord point; it may happen that
3208 # evaluated models do not cover the same seq coord range, hence matching
3209 # necessary)
3210 cube_iterables = [
3211 iris.cube.CubeList(
3212 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3213 )
3214 for point in all_points
3215 ]
3217 plot_index = []
3218 nplot = np.size(cube.coord(sequence_coordinate).points)
3219 # Create a plot for each value of the sequence coordinate. Allowing for
3220 # multiple cubes in a CubeList to be plotted in the same plot for similar
3221 # sequence values. Passing a CubeList into the internal plotting function
3222 # for similar values of the sequence coordinate. cube_slice can be an
3223 # iris.cube.Cube or an iris.cube.CubeList.
3224 for cube_slice in cube_iterables:
3225 single_cube = cube_slice
3226 if isinstance(cube_slice, iris.cube.CubeList): 3226 ↛ 3227line 3226 didn't jump to line 3227 because the condition on line 3226 was never true
3227 single_cube = cube_slice[0]
3229 # Ensure valid stamp coordinate in cube dimensions
3230 if stamp_coordinate == "realization": 3230 ↛ 3233line 3230 didn't jump to line 3233 because the condition on line 3230 was always true
3231 stamp_coordinate = check_stamp_coordinate(single_cube)
3232 # Set plot titles and filename, based on sequence coordinate
3233 seq_coord = single_cube.coord(sequence_coordinate)
3234 # Use time coordinate in title and filename if single histogram output.
3235 if sequence_coordinate == "realization" and nplot == 1: 3235 ↛ 3236line 3235 didn't jump to line 3236 because the condition on line 3235 was never true
3236 seq_coord = single_cube.coord("time")
3237 plot_title, plot_filename = _set_title_and_filename(
3238 seq_coord, nplot, recipe_title, filename
3239 )
3241 # Do the actual plotting.
3242 plotting_func(
3243 cube_slice,
3244 filename=plot_filename,
3245 stamp_coordinate=stamp_coordinate,
3246 title=plot_title,
3247 vmin=vmin,
3248 vmax=vmax,
3249 )
3250 plot_index.append(plot_filename)
3252 # Add list of plots to plot metadata.
3253 complete_plot_index = _append_to_plot_index(plot_index)
3255 # Make a page to display the plots.
3256 _make_plot_html_page(complete_plot_index)
3258 return cubes
3261def plot_power_spectrum_series(
3262 cubes: iris.cube.Cube | iris.cube.CubeList,
3263 filename: str = None,
3264 sequence_coordinate: str = "time",
3265 stamp_coordinate: str = "realization",
3266 single_plot: bool = False,
3267 **kwargs,
3268) -> iris.cube.Cube | iris.cube.CubeList:
3269 """Plot a power spectrum plot for each vertical level provided.
3271 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
3272 is present then a sequence of plots will be produced using the time slider
3273 functionality to scroll through power spectrum against time. If a
3274 stamp_coordinate is present then postage stamp plots will be produced. If
3275 stamp_coordinate and single_plot is True, all postage stamp plots will be
3276 plotted in a single plot instead of separate postage stamp plots.
3278 Parameters
3279 ----------
3280 cubes: Cube | iris.cube.CubeList
3281 Iris cube or CubeList of the data to plot. It should have a single dimension other
3282 than the stamp coordinate.
3283 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3284 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3285 filename: str, optional
3286 Name of the plot to write, used as a prefix for plot sequences. Defaults
3287 to the recipe name.
3288 sequence_coordinate: str, optional
3289 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3290 This coordinate must exist in the cube and will be used for the time
3291 slider.
3292 stamp_coordinate: str, optional
3293 Coordinate about which to plot postage stamp plots. Defaults to
3294 ``"realization"``.
3295 single_plot: bool, optional
3296 If True, all postage stamp plots will be plotted in a single plot. If
3297 False, each postage stamp plot will be plotted separately. Is only valid
3298 if stamp_coordinate exists and has more than a single point.
3300 Returns
3301 -------
3302 iris.cube.Cube | iris.cube.CubeList
3303 The original Cube or CubeList (so further operations can be applied).
3304 Plotted data.
3306 Raises
3307 ------
3308 ValueError
3309 If the cube doesn't have the right dimensions.
3310 TypeError
3311 If the cube isn't a Cube or CubeList.
3312 """
3313 recipe_title = get_recipe_metadata().get("title", "Untitled")
3315 cubes = iter_maybe(cubes)
3317 # Internal plotting function.
3318 plotting_func = _plot_and_save_power_spectrum_series
3320 num_models = _get_num_models(cubes)
3322 _validate_cube_shape(cubes, num_models)
3324 # If several power spectra are plotted with time as sequence_coordinate for the
3325 # time slider option.
3326 for cube in cubes:
3327 try:
3328 cube.coord(sequence_coordinate)
3329 except iris.exceptions.CoordinateNotFoundError as err:
3330 raise ValueError(
3331 f"Cube must have a {sequence_coordinate} coordinate."
3332 ) from err
3334 # Make postage stamp plots if stamp_coordinate exists and has more than a
3335 # single point. If single_plot is True:
3336 # -- all postage stamp plots will be plotted in a single plot instead of
3337 # separate postage stamp plots.
3338 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3339 # produced per single model only
3340 if num_models == 1: 3340 ↛ 3353line 3340 didn't jump to line 3353 because the condition on line 3340 was always true
3341 if ( 3341 ↛ 3345line 3341 didn't jump to line 3345 because the condition on line 3341 was never true
3342 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3343 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3344 ):
3345 if single_plot:
3346 plotting_func = (
3347 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3348 )
3349 else:
3350 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3351 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3352 else:
3353 all_points = sorted(
3354 set(
3355 itertools.chain.from_iterable(
3356 cb.coord(sequence_coordinate).points for cb in cubes
3357 )
3358 )
3359 )
3360 all_slices = list(
3361 itertools.chain.from_iterable(
3362 cb.slices_over(sequence_coordinate) for cb in cubes
3363 )
3364 )
3365 # Matched slices (matched by seq coord point; it may happen that
3366 # evaluated models do not cover the same seq coord range, hence matching
3367 # necessary)
3368 cube_iterables = [
3369 iris.cube.CubeList(
3370 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3371 )
3372 for point in all_points
3373 ]
3375 plot_index = []
3376 nplot = np.size(cube.coord(sequence_coordinate).points)
3377 # Create a plot for each value of the sequence coordinate. Allowing for
3378 # multiple cubes in a CubeList to be plotted in the same plot for similar
3379 # sequence values. Passing a CubeList into the internal plotting function
3380 # for similar values of the sequence coordinate. cube_slice can be an
3381 # iris.cube.Cube or an iris.cube.CubeList.
3382 for cube_slice in cube_iterables:
3383 single_cube = cube_slice
3384 if isinstance(cube_slice, iris.cube.CubeList): 3384 ↛ 3385line 3384 didn't jump to line 3385 because the condition on line 3384 was never true
3385 single_cube = cube_slice[0]
3387 # Set stamp coordinate
3388 if stamp_coordinate == "realization": 3388 ↛ 3391line 3388 didn't jump to line 3391 because the condition on line 3388 was always true
3389 stamp_coordinate = check_stamp_coordinate(single_cube)
3390 # Set plot title and filenames based on sequence values
3391 seq_coord = single_cube.coord(sequence_coordinate)
3392 plot_title, plot_filename = _set_title_and_filename(
3393 seq_coord, nplot, recipe_title, filename
3394 )
3396 # Do the actual plotting.
3397 plotting_func(
3398 cube_slice,
3399 filename=plot_filename,
3400 stamp_coordinate=stamp_coordinate,
3401 title=plot_title,
3402 )
3403 plot_index.append(plot_filename)
3405 # Add list of plots to plot metadata.
3406 complete_plot_index = _append_to_plot_index(plot_index)
3408 # Make a page to display the plots.
3409 _make_plot_html_page(complete_plot_index)
3411 return cubes
3414def _DCT_ps(y_3d):
3415 """Calculate power spectra for regional domains.
3417 Parameters
3418 ----------
3419 y_3d: 3D array
3420 3 dimensional array to calculate spectrum for.
3421 (2D field data with 3rd dimension of time)
3423 Returns: ps_array
3424 Array of power spectra values calculated for input field (for each time)
3426 Method for regional domains:
3427 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3428 as described in Denis et al 2002 [Denis_etal_2002]_.
3430 References
3431 ----------
3432 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3433 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3434 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3435 Monthly Weather Review, Vol. 130, 1812-1828
3436 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3437 """
3438 Nt, Ny, Nx = y_3d.shape
3440 # Max coefficient
3441 Nmin = min(Nx - 1, Ny - 1)
3443 # Create alpha matrix (of wavenumbers)
3444 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3446 # Prepare output array
3447 ps_array = np.zeros((Nt, Nmin))
3449 # Loop over time to get spectrum for each time.
3450 for t in range(Nt):
3451 y_2d = y_3d[t]
3453 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3454 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3455 # cosine basis functions at different spatial frequencies.
3457 # normalise spectrum to allow comparison between models.
3458 fkk = fft.dctn(y_2d, norm="ortho")
3460 # Normalise fkk
3461 fkk = fkk / np.sqrt(Ny * Nx)
3463 # calculate variance of spectral coefficient
3464 sigma_2 = fkk**2 / Nx / Ny
3466 # Group ellipses of alphas into the same wavenumber k/Nmin
3467 for k in range(1, Nmin + 1):
3468 alpha = k / Nmin
3469 alpha_p1 = (k + 1) / Nmin
3471 # Sum up elements matching k
3472 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3473 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3475 return ps_array
3478def _create_alpha_matrix(Ny, Nx):
3479 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3481 Parameters
3482 ----------
3483 Ny, Nx:
3484 Dimensions of the 2D field for which the power spectra is calculated. Used to
3485 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3486 single-scale parameter.
3488 Returns: alpha_matrix
3489 normalisation of 2D wavenumber axes, transforming the spectral domain into
3490 an elliptic coordinate system.
3492 """
3493 # Create x_indices: each row is [1, 2, ..., Nx]
3494 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3496 # Create y_indices: each column is [1, 2, ..., Ny]
3497 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3499 # Compute alpha_matrix
3500 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3502 return alpha_matrix