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