Coverage for src / CSET / operators / plot.py: 86%
968 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-11-28 12:47 +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
25import sys
26from typing import Literal
28import cartopy.crs as ccrs
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 iris.cube import Cube
40from markdown_it import MarkdownIt
42from CSET._common import (
43 combine_dicts,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._utils import (
50 fully_equalise_attributes,
51 get_cube_yxcoordname,
52 is_transect,
53)
54from CSET.operators.collapse import collapse
55from CSET.operators.misc import _extract_common_time_points
56from CSET.operators.regrid import regrid_onto_cube
58# Use a non-interactive plotting backend.
59mpl.use("agg")
61DEFAULT_DISCRETE_COLORS = mpl.colormaps["tab10"].colors + mpl.colormaps["Accent"].colors
63############################
64# Private helper functions #
65############################
68def _append_to_plot_index(plot_index: list) -> list:
69 """Add plots into the plot index, returning the complete plot index."""
70 with open("meta.json", "r+t", encoding="UTF-8") as fp:
71 fcntl.flock(fp, fcntl.LOCK_EX)
72 fp.seek(0)
73 meta = json.load(fp)
74 complete_plot_index = meta.get("plots", [])
75 complete_plot_index = complete_plot_index + plot_index
76 meta["plots"] = complete_plot_index
77 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
78 os.getenv("DO_CASE_AGGREGATION")
79 ):
80 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
81 fp.seek(0)
82 fp.truncate()
83 json.dump(meta, fp, indent=2)
84 return complete_plot_index
87def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:
88 """Ensure a single cube is given.
90 If a CubeList of length one is given that the contained cube is returned,
91 otherwise an error is raised.
93 Parameters
94 ----------
95 cube: Cube | CubeList
96 The cube to check.
98 Returns
99 -------
100 cube: Cube
101 The checked cube.
103 Raises
104 ------
105 TypeError
106 If the input cube is not a Cube or CubeList of a single Cube.
107 """
108 if isinstance(cube, iris.cube.Cube):
109 return cube
110 if isinstance(cube, iris.cube.CubeList):
111 if len(cube) == 1:
112 return cube[0]
113 raise TypeError("Must have a single cube", cube)
116def _py312_importlib_resources_files_shim():
117 """Importlib behaviour changed in 3.12 to avoid circular dependencies.
119 This shim is needed until python 3.12 is our oldest supported version, after
120 which it can just be replaced by directly using importlib.resources.files.
121 """
122 if sys.version_info.minor >= 12:
123 files = importlib.resources.files()
124 else:
125 import CSET.operators
127 files = importlib.resources.files(CSET.operators)
128 return files
131def _make_plot_html_page(plots: list):
132 """Create a HTML page to display a plot image."""
133 # Debug check that plots actually contains some strings.
134 assert isinstance(plots[0], str)
136 # Load HTML template file.
137 operator_files = _py312_importlib_resources_files_shim()
138 template_file = operator_files.joinpath("_plot_page_template.html")
140 # Get some metadata.
141 meta = get_recipe_metadata()
142 title = meta.get("title", "Untitled")
143 description = MarkdownIt().render(meta.get("description", "*No description.*"))
145 # Prepare template variables.
146 variables = {
147 "title": title,
148 "description": description,
149 "initial_plot": plots[0],
150 "plots": plots,
151 "title_slug": slugify(title),
152 }
154 # Render template.
155 html = render_file(template_file, **variables)
157 # Save completed HTML.
158 with open("index.html", "wt", encoding="UTF-8") as fp:
159 fp.write(html)
162@functools.cache
163def _load_colorbar_map(user_colorbar_file: str = None) -> dict:
164 """Load the colorbar definitions from a file.
166 This is a separate function to make it cacheable.
167 """
168 colorbar_file = _py312_importlib_resources_files_shim().joinpath(
169 "_colorbar_definition.json"
170 )
171 with open(colorbar_file, "rt", encoding="UTF-8") as fp:
172 colorbar = json.load(fp)
174 logging.debug("User colour bar file: %s", user_colorbar_file)
175 override_colorbar = {}
176 if user_colorbar_file:
177 try:
178 with open(user_colorbar_file, "rt", encoding="UTF-8") as fp:
179 override_colorbar = json.load(fp)
180 except FileNotFoundError:
181 logging.warning("Colorbar file does not exist. Using default values.")
183 # Overwrite values with the user supplied colorbar definition.
184 colorbar = combine_dicts(colorbar, override_colorbar)
185 return colorbar
188def _get_model_colors_map(cubes: iris.cube.CubeList | iris.cube.Cube) -> dict:
189 """Get an appropriate colors for model lines in line plots.
191 For each model in the list of cubes colors either from user provided
192 color definition file (so-called style file) or from default colors are mapped
193 to model_name attribute.
195 Parameters
196 ----------
197 cubes: CubeList or Cube
198 Cubes with model_name attribute
200 Returns
201 -------
202 model_colors_map:
203 Dictionary mapping model_name attribute to colors
204 """
205 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
206 colorbar = _load_colorbar_map(user_colorbar_file)
207 model_names = sorted(
208 filter(
209 lambda x: x is not None,
210 (cube.attributes.get("model_name", None) for cube in iter_maybe(cubes)),
211 )
212 )
213 if not model_names:
214 return {}
215 use_user_colors = all(mname in colorbar.keys() for mname in model_names)
216 if use_user_colors: 216 ↛ 217line 216 didn't jump to line 217 because the condition on line 216 was never true
217 return {mname: colorbar[mname] for mname in model_names}
219 color_list = itertools.cycle(DEFAULT_DISCRETE_COLORS)
220 return {mname: color for mname, color in zip(model_names, color_list, strict=False)}
223def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
224 """Get an appropriate colorbar for the given cube.
226 For the given variable the appropriate colorbar is looked up from a
227 combination of the built-in CSET colorbar definitions, and any user supplied
228 definitions. As well as varying on variables, these definitions may also
229 exist for specific pressure levels to account for variables with
230 significantly different ranges at different heights. The colorbars also exist
231 for masks and mask differences for considering variable presence diagnostics.
232 Specific variable ranges can be separately set in user-supplied definition
233 for x- or y-axis limits, or indicate where automated range preferred.
235 Parameters
236 ----------
237 cube: Cube
238 Cube of variable for which the colorbar information is desired.
239 axis: "x", "y", optional
240 Select the levels for just this axis of a line plot. The min and max
241 can be set by xmin/xmax or ymin/ymax respectively. For variables where
242 setting a universal range is not desirable (e.g. temperature), users
243 can set ymin/ymax values to "auto" in the colorbar definitions file.
244 Where no additional xmin/xmax or ymin/ymax values are provided, the
245 axis bounds default to use the vmin/vmax values provided.
247 Returns
248 -------
249 cmap:
250 Matplotlib colormap.
251 levels:
252 List of levels to use for plotting. For continuous plots the min and max
253 should be taken as the range.
254 norm:
255 BoundaryNorm information.
256 """
257 # Grab the colorbar file from the recipe global metadata.
258 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
259 colorbar = _load_colorbar_map(user_colorbar_file)
260 cmap = None
262 try:
263 # We assume that pressure is a scalar coordinate here.
264 pressure_level_raw = cube.coord("pressure").points[0]
265 # Ensure pressure_level is a string, as it is used as a JSON key.
266 pressure_level = str(int(pressure_level_raw))
267 except iris.exceptions.CoordinateNotFoundError:
268 pressure_level = None
270 # First try long name, then standard name, then var name. This order is used
271 # as long name is the one we correct between models, so it most likely to be
272 # consistent.
273 varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name]))
274 for varname in varnames:
275 # Get the colormap for this variable.
276 try:
277 var_colorbar = colorbar[varname]
278 cmap = plt.get_cmap(colorbar[varname]["cmap"], 51)
279 varname_key = varname
280 break
281 except KeyError:
282 logging.debug("Cube name %s has no colorbar definition.", varname)
284 # Get colormap if it is a mask.
285 if any("mask_for_" in name for name in varnames):
286 cmap, levels, norm = _custom_colormap_mask(cube, axis=axis)
287 return cmap, levels, norm
288 # If winds on Beaufort Scale use custom colorbar and levels
289 if any("Beaufort_Scale" in name for name in varnames):
290 cmap, levels, norm = _custom_beaufort_scale(cube, axis=axis)
291 return cmap, levels, norm
292 # If probability is plotted use custom colorbar and levels
293 if any("probability_of_" in name for name in varnames):
294 cmap, levels, norm = _custom_colormap_probability(cube, axis=axis)
295 return cmap, levels, norm
296 # If aviation colour state use custom colorbar and levels
297 if any("aviation_colour_state" in name for name in varnames): 297 ↛ 298line 297 didn't jump to line 298 because the condition on line 297 was never true
298 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)
299 return cmap, levels, norm
301 # If no valid colormap has been defined, use defaults and return.
302 if not cmap:
303 logging.warning("No colorbar definition exists for %s.", cube.name())
304 cmap, levels, norm = mpl.colormaps["viridis"], None, None
305 return cmap, levels, norm
307 # Test if pressure-level specific settings are provided for cube.
308 if pressure_level:
309 try:
310 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]
311 except KeyError:
312 logging.debug(
313 "%s has no colorbar definition for pressure level %s.",
314 varname,
315 pressure_level,
316 )
318 # Check for availability of x-axis or y-axis user-specific overrides
319 # for setting level bounds for line plot types and return just levels.
320 # Line plots do not need a colormap, and just use the data range.
321 if axis:
322 if axis == "x":
323 try:
324 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]
325 except KeyError:
326 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
327 if axis == "y":
328 try:
329 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]
330 except KeyError:
331 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
332 # Check if user-specified auto-scaling for this variable
333 if vmin == "auto" or vmax == "auto":
334 levels = None
335 else:
336 levels = [vmin, vmax]
337 return None, levels, None
338 # Get and use the colorbar levels for this variable if spatial or histogram.
339 else:
340 try:
341 levels = var_colorbar["levels"]
342 # Use discrete bins when levels are specified, rather
343 # than a smooth range.
344 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
345 logging.debug("Using levels for %s colorbar.", varname)
346 logging.info("Using levels: %s", levels)
347 except KeyError:
348 # Get the range for this variable.
349 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
350 logging.debug("Using min and max for %s colorbar.", varname)
351 # Calculate levels from range.
352 levels = np.linspace(vmin, vmax, 101)
353 norm = None
355 # Overwrite cmap, levels and norm for specific variables that
356 # require custom colorbar_map as these can not be defined in the
357 # JSON file.
358 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)
359 cmap, levels, norm = _custom_colourmap_visibility_in_air(
360 cube, cmap, levels, norm
361 )
362 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)
363 return cmap, levels, norm
366def _setup_spatial_map(
367 cube: iris.cube.Cube,
368 figure,
369 cmap,
370 grid_size: int | None = None,
371 subplot: int | None = None,
372):
373 """Define map projections, extent and add coastlines for spatial plots.
375 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
376 is specified, and map extent defined based on the input data.
378 Parameters
379 ----------
380 cube: Cube
381 2 dimensional (lat and lon) Cube of the data to plot.
382 figure:
383 Matplotlib Figure object holding all plot elements.
384 cmap:
385 Matplotlib colormap.
386 grid_size: int, optional
387 Size of grid for subplots if multiple spatial subplots in figure.
388 subplot: int, optional
389 Subplot index if multiple spatial subplots in figure.
391 Returns
392 -------
393 axes:
394 Matplotlib GeoAxes definition.
395 """
396 # Identify min/max plot bounds.
397 try:
398 lat_axis, lon_axis = get_cube_yxcoordname(cube)
399 x1 = np.min(cube.coord(lon_axis).points)
400 x2 = np.max(cube.coord(lon_axis).points)
401 y1 = np.min(cube.coord(lat_axis).points)
402 y2 = np.max(cube.coord(lat_axis).points)
404 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
405 if np.abs(x2 - x1) > 180.0:
406 x1 = x1 - 180.0
407 x2 = x2 - 180.0
408 logging.debug("Adjusting plot bounds to fit global extent.")
410 # Consider map projection orientation.
411 # Adapting orientation enables plotting across international dateline.
412 # Users can adapt the default central_longitude if alternative projections views.
413 if x2 > 180.0:
414 central_longitude = 180.0
415 else:
416 central_longitude = 0.0
418 # Define spatial map projection.
419 coord_system = cube.coord(lat_axis).coord_system
420 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
421 # Define rotated pole map projection for rotated pole inputs.
422 projection = ccrs.RotatedPole(
423 pole_longitude=coord_system.grid_north_pole_longitude,
424 pole_latitude=coord_system.grid_north_pole_latitude,
425 central_rotated_longitude=0.0,
426 )
427 crs = projection
428 else:
429 # Define regular map projection for non-rotated pole inputs.
430 # Alternatives might include e.g. for global model outputs:
431 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
432 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
433 projection = ccrs.PlateCarree(central_longitude=central_longitude)
434 crs = ccrs.PlateCarree()
436 # Define axes for plot (or subplot) with required map projection.
437 if subplot is not None:
438 axes = figure.add_subplot(
439 grid_size, grid_size, subplot, projection=projection
440 )
441 else:
442 axes = figure.add_subplot(projection=projection)
444 # Add coastlines if cube contains x and y map coordinates.
445 if cmap.name in ["viridis", "Greys"]:
446 coastcol = "magenta"
447 else:
448 coastcol = "black"
449 logging.debug("Plotting coastlines in colour %s.", coastcol)
450 axes.coastlines(resolution="10m", color=coastcol)
452 # If is lat/lon spatial map, fix extent to keep plot tight.
453 # Specifying crs within set_extent helps ensure only data region is shown.
454 if isinstance(coord_system, iris.coord_systems.GeogCS):
455 axes.set_extent([x1, x2, y1, y2], crs=crs)
457 except ValueError:
458 # Skip if not both x and y map coordinates.
459 axes = figure.gca()
460 pass
462 return axes
465def _get_plot_resolution() -> int:
466 """Get resolution of rasterised plots in pixels per inch."""
467 return get_recipe_metadata().get("plot_resolution", 100)
470def _plot_and_save_spatial_plot(
471 cube: iris.cube.Cube,
472 filename: str,
473 title: str,
474 method: Literal["contourf", "pcolormesh"],
475 **kwargs,
476):
477 """Plot and save a spatial plot.
479 Parameters
480 ----------
481 cube: Cube
482 2 dimensional (lat and lon) Cube of the data to plot.
483 filename: str
484 Filename of the plot to write.
485 title: str
486 Plot title.
487 method: "contourf" | "pcolormesh"
488 The plotting method to use.
489 """
490 # Setup plot details, size, resolution, etc.
491 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
493 # Specify the color bar
494 cmap, levels, norm = _colorbar_map_levels(cube)
496 # Setup plot map projection, extent and coastlines.
497 axes = _setup_spatial_map(cube, fig, cmap)
499 # Plot the field.
500 if method == "contourf":
501 # Filled contour plot of the field.
502 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
503 elif method == "pcolormesh":
504 try:
505 vmin = min(levels)
506 vmax = max(levels)
507 except TypeError:
508 vmin, vmax = None, None
509 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
510 # if levels are defined.
511 if norm is not None:
512 vmin = None
513 vmax = None
514 logging.debug("Plotting using defined levels.")
515 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
516 else:
517 raise ValueError(f"Unknown plotting method: {method}")
519 # Check to see if transect, and if so, adjust y axis.
520 if is_transect(cube):
521 if "pressure" in [coord.name() for coord in cube.coords()]:
522 axes.invert_yaxis()
523 axes.set_yscale("log")
524 axes.set_ylim(1100, 100)
525 # If both model_level_number and level_height exists, iplt can construct
526 # plot as a function of height above orography (NOT sea level).
527 elif {"model_level_number", "level_height"}.issubset( 527 ↛ 532line 527 didn't jump to line 532 because the condition on line 527 was always true
528 {coord.name() for coord in cube.coords()}
529 ):
530 axes.set_yscale("log")
532 axes.set_title(
533 f"{title}\n"
534 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
535 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
536 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
537 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
538 fontsize=16,
539 )
541 else:
542 # Add title.
543 axes.set_title(title, fontsize=16)
545 # Add watermark with min/max/mean. Currently not user togglable.
546 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
547 axes.annotate(
548 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
549 xy=(1, -0.05),
550 xycoords="axes fraction",
551 xytext=(-5, 5),
552 textcoords="offset points",
553 ha="right",
554 va="bottom",
555 size=11,
556 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
557 )
559 # Add colour bar.
560 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
561 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
562 # add ticks and tick_labels for every levels if less than 20 levels exist
563 if levels is not None and len(levels) < 20:
564 cbar.set_ticks(levels)
565 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
566 if "visibility" in cube.name(): 566 ↛ 567line 566 didn't jump to line 567 because the condition on line 566 was never true
567 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
568 logging.debug("Set colorbar ticks and labels.")
570 # Save plot.
571 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
572 logging.info("Saved spatial plot to %s", filename)
573 plt.close(fig)
576def _plot_and_save_postage_stamp_spatial_plot(
577 cube: iris.cube.Cube,
578 filename: str,
579 stamp_coordinate: str,
580 title: str,
581 method: Literal["contourf", "pcolormesh"],
582 **kwargs,
583):
584 """Plot postage stamp spatial plots from an ensemble.
586 Parameters
587 ----------
588 cube: Cube
589 Iris cube of data to be plotted. It must have the stamp coordinate.
590 filename: str
591 Filename of the plot to write.
592 stamp_coordinate: str
593 Coordinate that becomes different plots.
594 method: "contourf" | "pcolormesh"
595 The plotting method to use.
597 Raises
598 ------
599 ValueError
600 If the cube doesn't have the right dimensions.
601 """
602 # Use the smallest square grid that will fit the members.
603 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
605 fig = plt.figure(figsize=(10, 10))
607 # Specify the color bar
608 cmap, levels, norm = _colorbar_map_levels(cube)
610 # Make a subplot for each member.
611 for member, subplot in zip(
612 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
613 ):
614 # Setup subplot map projection, extent and coastlines.
615 axes = _setup_spatial_map(
616 member, fig, cmap, grid_size=grid_size, subplot=subplot
617 )
618 if method == "contourf":
619 # Filled contour plot of the field.
620 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
621 elif method == "pcolormesh":
622 if levels is not None:
623 vmin = min(levels)
624 vmax = max(levels)
625 else:
626 raise TypeError("Unknown vmin and vmax range.")
627 vmin, vmax = None, None
628 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
629 # if levels are defined.
630 if norm is not None: 630 ↛ 631line 630 didn't jump to line 631 because the condition on line 630 was never true
631 vmin = None
632 vmax = None
633 # pcolormesh plot of the field.
634 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
635 else:
636 raise ValueError(f"Unknown plotting method: {method}")
637 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
638 axes.set_axis_off()
640 # Put the shared colorbar in its own axes.
641 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])
642 colorbar = fig.colorbar(
643 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
644 )
645 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
647 # Overall figure title.
648 fig.suptitle(title, fontsize=16)
650 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
651 logging.info("Saved contour postage stamp plot to %s", filename)
652 plt.close(fig)
655def _plot_and_save_line_series(
656 cubes: iris.cube.CubeList,
657 coords: list[iris.coords.Coord],
658 ensemble_coord: str,
659 filename: str,
660 title: str,
661 **kwargs,
662):
663 """Plot and save a 1D line series.
665 Parameters
666 ----------
667 cubes: Cube or CubeList
668 Cube or CubeList containing the cubes to plot on the y-axis.
669 coords: list[Coord]
670 Coordinates to plot on the x-axis, one per cube.
671 ensemble_coord: str
672 Ensemble coordinate in the cube.
673 filename: str
674 Filename of the plot to write.
675 title: str
676 Plot title.
677 """
678 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
680 model_colors_map = _get_model_colors_map(cubes)
682 # Store min/max ranges.
683 y_levels = []
685 # Check match-up across sequence coords gives consistent sizes
686 _validate_cubes_coords(cubes, coords)
688 for cube, coord in zip(cubes, coords, strict=True):
689 label = None
690 color = "black"
691 if model_colors_map:
692 label = cube.attributes.get("model_name")
693 color = model_colors_map.get(label)
694 for cube_slice in cube.slices_over(ensemble_coord):
695 # Label with (control) if part of an ensemble or not otherwise.
696 if cube_slice.coord(ensemble_coord).points == [0]:
697 iplt.plot(
698 coord,
699 cube_slice,
700 color=color,
701 marker="o",
702 ls="-",
703 lw=3,
704 label=f"{label} (control)"
705 if len(cube.coord(ensemble_coord).points) > 1
706 else label,
707 )
708 # Label with (perturbed) if part of an ensemble and not the control.
709 else:
710 iplt.plot(
711 coord,
712 cube_slice,
713 color=color,
714 ls="-",
715 lw=1.5,
716 alpha=0.75,
717 label=f"{label} (member)",
718 )
720 # Calculate the global min/max if multiple cubes are given.
721 _, levels, _ = _colorbar_map_levels(cube, axis="y")
722 if levels is not None: 722 ↛ 723line 722 didn't jump to line 723 because the condition on line 722 was never true
723 y_levels.append(min(levels))
724 y_levels.append(max(levels))
726 # Get the current axes.
727 ax = plt.gca()
729 # Add some labels and tweak the style.
730 # check if cubes[0] works for single cube if not CubeList
731 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
732 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
733 ax.set_title(title, fontsize=16)
735 ax.ticklabel_format(axis="y", useOffset=False)
736 ax.tick_params(axis="x", labelrotation=15)
737 ax.tick_params(axis="both", labelsize=12)
739 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
740 if y_levels: 740 ↛ 741line 740 didn't jump to line 741 because the condition on line 740 was never true
741 ax.set_ylim(min(y_levels), max(y_levels))
742 # Add zero line.
743 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
744 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
745 logging.debug(
746 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
747 )
748 else:
749 ax.autoscale()
751 # Add gridlines
752 ax.grid(linestyle="--", color="grey", linewidth=1)
753 # Ientify unique labels for legend
754 handles = list(
755 {
756 label: handle
757 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
758 }.values()
759 )
760 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
762 # Save plot.
763 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
764 logging.info("Saved line plot to %s", filename)
765 plt.close(fig)
768def _plot_and_save_vertical_line_series(
769 cubes: iris.cube.CubeList,
770 coords: list[iris.coords.Coord],
771 ensemble_coord: str,
772 filename: str,
773 series_coordinate: str,
774 title: str,
775 vmin: float,
776 vmax: float,
777 **kwargs,
778):
779 """Plot and save a 1D line series in vertical.
781 Parameters
782 ----------
783 cubes: CubeList
784 1 dimensional Cube or CubeList of the data to plot on x-axis.
785 coord: list[Coord]
786 Coordinates to plot on the y-axis, one per cube.
787 ensemble_coord: str
788 Ensemble coordinate in the cube.
789 filename: str
790 Filename of the plot to write.
791 series_coordinate: str
792 Coordinate to use as vertical axis.
793 title: str
794 Plot title.
795 vmin: float
796 Minimum value for the x-axis.
797 vmax: float
798 Maximum value for the x-axis.
799 """
800 # plot the vertical pressure axis using log scale
801 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
803 model_colors_map = _get_model_colors_map(cubes)
805 # Check match-up across sequence coords gives consistent sizes
806 _validate_cubes_coords(cubes, coords)
808 for cube, coord in zip(cubes, coords, strict=True):
809 label = None
810 color = "black"
811 if model_colors_map: 811 ↛ 812line 811 didn't jump to line 812 because the condition on line 811 was never true
812 label = cube.attributes.get("model_name")
813 color = model_colors_map.get(label)
815 for cube_slice in cube.slices_over(ensemble_coord):
816 # If ensemble data given plot control member with (control)
817 # unless single forecast.
818 if cube_slice.coord(ensemble_coord).points == [0]:
819 iplt.plot(
820 cube_slice,
821 coord,
822 color=color,
823 marker="o",
824 ls="-",
825 lw=3,
826 label=f"{label} (control)"
827 if len(cube.coord(ensemble_coord).points) > 1
828 else label,
829 )
830 # If ensemble data given plot perturbed members with (perturbed).
831 else:
832 iplt.plot(
833 cube_slice,
834 coord,
835 color=color,
836 ls="-",
837 lw=1.5,
838 alpha=0.75,
839 label=f"{label} (member)",
840 )
842 # Get the current axis
843 ax = plt.gca()
845 # Special handling for pressure level data.
846 if series_coordinate == "pressure": 846 ↛ 868line 846 didn't jump to line 868 because the condition on line 846 was always true
847 # Invert y-axis and set to log scale.
848 ax.invert_yaxis()
849 ax.set_yscale("log")
851 # Define y-ticks and labels for pressure log axis.
852 y_tick_labels = [
853 "1000",
854 "850",
855 "700",
856 "500",
857 "300",
858 "200",
859 "100",
860 ]
861 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
863 # Set y-axis limits and ticks.
864 ax.set_ylim(1100, 100)
866 # Test if series_coordinate is model level data. The UM data uses
867 # model_level_number and lfric uses full_levels as coordinate.
868 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
869 # Define y-ticks and labels for vertical axis.
870 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
871 y_tick_labels = [str(int(i)) for i in y_ticks]
872 ax.set_ylim(min(y_ticks), max(y_ticks))
874 ax.set_yticks(y_ticks)
875 ax.set_yticklabels(y_tick_labels)
877 # Set x-axis limits.
878 ax.set_xlim(vmin, vmax)
879 # Mark y=0 if present in plot.
880 if vmin < 0.0 and vmax > 0.0: 880 ↛ 881line 880 didn't jump to line 881 because the condition on line 880 was never true
881 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
883 # Add some labels and tweak the style.
884 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
885 ax.set_xlabel(
886 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
887 )
888 ax.set_title(title, fontsize=16)
889 ax.ticklabel_format(axis="x")
890 ax.tick_params(axis="y")
891 ax.tick_params(axis="both", labelsize=12)
893 # Add gridlines
894 ax.grid(linestyle="--", color="grey", linewidth=1)
895 # Ientify unique labels for legend
896 handles = list(
897 {
898 label: handle
899 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
900 }.values()
901 )
902 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
904 # Save plot.
905 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
906 logging.info("Saved line plot to %s", filename)
907 plt.close(fig)
910def _plot_and_save_scatter_plot(
911 cube_x: iris.cube.Cube | iris.cube.CubeList,
912 cube_y: iris.cube.Cube | iris.cube.CubeList,
913 filename: str,
914 title: str,
915 one_to_one: bool,
916 model_names: list[str] = None,
917 **kwargs,
918):
919 """Plot and save a 2D scatter plot.
921 Parameters
922 ----------
923 cube_x: Cube | CubeList
924 1 dimensional Cube or CubeList of the data to plot on x-axis.
925 cube_y: Cube | CubeList
926 1 dimensional Cube or CubeList of the data to plot on y-axis.
927 filename: str
928 Filename of the plot to write.
929 title: str
930 Plot title.
931 one_to_one: bool
932 Whether a 1:1 line is plotted.
933 """
934 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
935 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
936 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
937 # over the pairs simultaneously.
939 # Ensure cube_x and cube_y are iterable
940 cube_x_iterable = iter_maybe(cube_x)
941 cube_y_iterable = iter_maybe(cube_y)
943 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
944 iplt.scatter(cube_x_iter, cube_y_iter)
945 if one_to_one is True:
946 plt.plot(
947 [
948 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
949 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
950 ],
951 [
952 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
953 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
954 ],
955 "k",
956 linestyle="--",
957 )
958 ax = plt.gca()
960 # Add some labels and tweak the style.
961 if model_names is None:
962 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
963 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
964 else:
965 # Add the model names, these should be order of base (x) and other (y).
966 ax.set_xlabel(
967 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
968 )
969 ax.set_ylabel(
970 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
971 )
972 ax.set_title(title, fontsize=16)
973 ax.ticklabel_format(axis="y", useOffset=False)
974 ax.tick_params(axis="x", labelrotation=15)
975 ax.tick_params(axis="both", labelsize=12)
976 ax.autoscale()
978 # Save plot.
979 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
980 logging.info("Saved scatter plot to %s", filename)
981 plt.close(fig)
984def _plot_and_save_vector_plot(
985 cube_u: iris.cube.Cube,
986 cube_v: iris.cube.Cube,
987 filename: str,
988 title: str,
989 method: Literal["contourf", "pcolormesh"],
990 **kwargs,
991):
992 """Plot and save a 2D vector plot.
994 Parameters
995 ----------
996 cube_u: Cube
997 2 dimensional Cube of u component of the data.
998 cube_v: Cube
999 2 dimensional Cube of v component of the data.
1000 filename: str
1001 Filename of the plot to write.
1002 title: str
1003 Plot title.
1004 """
1005 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1007 # Create a cube containing the magnitude of the vector field.
1008 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1009 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1011 # Specify the color bar
1012 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1014 # Setup plot map projection, extent and coastlines.
1015 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1017 if method == "contourf": 1017 ↛ 1020line 1017 didn't jump to line 1020 because the condition on line 1017 was always true
1018 # Filled contour plot of the field.
1019 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1020 elif method == "pcolormesh":
1021 try:
1022 vmin = min(levels)
1023 vmax = max(levels)
1024 except TypeError:
1025 vmin, vmax = None, None
1026 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1027 # if levels are defined.
1028 if norm is not None:
1029 vmin = None
1030 vmax = None
1031 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1032 else:
1033 raise ValueError(f"Unknown plotting method: {method}")
1035 # Check to see if transect, and if so, adjust y axis.
1036 if is_transect(cube_vec_mag): 1036 ↛ 1037line 1036 didn't jump to line 1037 because the condition on line 1036 was never true
1037 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1038 axes.invert_yaxis()
1039 axes.set_yscale("log")
1040 axes.set_ylim(1100, 100)
1041 # If both model_level_number and level_height exists, iplt can construct
1042 # plot as a function of height above orography (NOT sea level).
1043 elif {"model_level_number", "level_height"}.issubset(
1044 {coord.name() for coord in cube_vec_mag.coords()}
1045 ):
1046 axes.set_yscale("log")
1048 axes.set_title(
1049 f"{title}\n"
1050 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1051 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1052 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1053 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1054 fontsize=16,
1055 )
1057 else:
1058 # Add title.
1059 axes.set_title(title, fontsize=16)
1061 # Add watermark with min/max/mean. Currently not user togglable.
1062 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1063 axes.annotate(
1064 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}",
1065 xy=(1, -0.05),
1066 xycoords="axes fraction",
1067 xytext=(-5, 5),
1068 textcoords="offset points",
1069 ha="right",
1070 va="bottom",
1071 size=11,
1072 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1073 )
1075 # Add colour bar.
1076 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1077 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1078 # add ticks and tick_labels for every levels if less than 20 levels exist
1079 if levels is not None and len(levels) < 20: 1079 ↛ 1080line 1079 didn't jump to line 1080 because the condition on line 1079 was never true
1080 cbar.set_ticks(levels)
1081 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1083 # 30 barbs along the longest axis of the plot, or a barb per point for data
1084 # with less than 30 points.
1085 step = max(max(cube_u.shape) // 30, 1)
1086 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1088 # Save plot.
1089 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1090 logging.info("Saved vector plot to %s", filename)
1091 plt.close(fig)
1094def _plot_and_save_histogram_series(
1095 cubes: iris.cube.Cube | iris.cube.CubeList,
1096 filename: str,
1097 title: str,
1098 vmin: float,
1099 vmax: float,
1100 **kwargs,
1101):
1102 """Plot and save a histogram series.
1104 Parameters
1105 ----------
1106 cubes: Cube or CubeList
1107 2 dimensional Cube or CubeList of the data to plot as histogram.
1108 filename: str
1109 Filename of the plot to write.
1110 title: str
1111 Plot title.
1112 vmin: float
1113 minimum for colorbar
1114 vmax: float
1115 maximum for colorbar
1116 """
1117 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1118 ax = plt.gca()
1120 model_colors_map = _get_model_colors_map(cubes)
1122 # Set default that histograms will produce probability density function
1123 # at each bin (integral over range sums to 1).
1124 density = True
1126 for cube in iter_maybe(cubes):
1127 # Easier to check title (where var name originates)
1128 # than seeing if long names exist etc.
1129 # Exception case, where distribution better fits log scales/bins.
1130 if "surface_microphysical" in title:
1131 if "amount" in title: 1131 ↛ 1133line 1131 didn't jump to line 1133 because the condition on line 1131 was never true
1132 # Compute histogram following Klingaman et al. (2017): ASoP
1133 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1134 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1135 density = False
1136 else:
1137 bins = 10.0 ** (
1138 np.arange(-10, 27, 1) / 10.0
1139 ) # Suggestion from RMED toolbox.
1140 bins = np.insert(bins, 0, 0)
1141 ax.set_yscale("log")
1142 vmin = bins[1]
1143 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1144 ax.set_xscale("log")
1145 elif "lightning" in title:
1146 bins = [0, 1, 2, 3, 4, 5]
1147 else:
1148 bins = np.linspace(vmin, vmax, 51)
1149 logging.debug(
1150 "Plotting histogram with %s bins %s - %s.",
1151 np.size(bins),
1152 np.min(bins),
1153 np.max(bins),
1154 )
1156 # Reshape cube data into a single array to allow for a single histogram.
1157 # Otherwise we plot xdim histograms stacked.
1158 cube_data_1d = (cube.data).flatten()
1160 label = None
1161 color = "black"
1162 if model_colors_map: 1162 ↛ 1163line 1162 didn't jump to line 1163 because the condition on line 1162 was never true
1163 label = cube.attributes.get("model_name")
1164 color = model_colors_map[label]
1165 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1167 # Compute area under curve.
1168 if "surface_microphysical" in title and "amount" in title: 1168 ↛ 1169line 1168 didn't jump to line 1169 because the condition on line 1168 was never true
1169 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1170 x = x * bin_mean / x.sum()
1171 x = x[1:]
1172 y = y[1:]
1174 ax.plot(
1175 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1176 )
1178 # Add some labels and tweak the style.
1179 ax.set_title(title, fontsize=16)
1180 ax.set_xlabel(
1181 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1182 )
1183 ax.set_ylabel("Normalised probability density", fontsize=14)
1184 if "surface_microphysical" in title and "amount" in title: 1184 ↛ 1185line 1184 didn't jump to line 1185 because the condition on line 1184 was never true
1185 ax.set_ylabel(
1186 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1187 )
1188 ax.set_xlim(vmin, vmax)
1189 ax.tick_params(axis="both", labelsize=12)
1191 # Overlay grid-lines onto histogram plot.
1192 ax.grid(linestyle="--", color="grey", linewidth=1)
1193 if model_colors_map: 1193 ↛ 1194line 1193 didn't jump to line 1194 because the condition on line 1193 was never true
1194 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1196 # Save plot.
1197 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1198 logging.info("Saved line plot to %s", filename)
1199 plt.close(fig)
1202def _plot_and_save_postage_stamp_histogram_series(
1203 cube: iris.cube.Cube,
1204 filename: str,
1205 title: str,
1206 stamp_coordinate: str,
1207 vmin: float,
1208 vmax: float,
1209 **kwargs,
1210):
1211 """Plot and save postage (ensemble members) stamps for a histogram series.
1213 Parameters
1214 ----------
1215 cube: Cube
1216 2 dimensional Cube of the data to plot as histogram.
1217 filename: str
1218 Filename of the plot to write.
1219 title: str
1220 Plot title.
1221 stamp_coordinate: str
1222 Coordinate that becomes different plots.
1223 vmin: float
1224 minimum for pdf x-axis
1225 vmax: float
1226 maximum for pdf x-axis
1227 """
1228 # Use the smallest square grid that will fit the members.
1229 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1231 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1232 # Make a subplot for each member.
1233 for member, subplot in zip(
1234 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1235 ):
1236 # Implicit interface is much easier here, due to needing to have the
1237 # cartopy GeoAxes generated.
1238 plt.subplot(grid_size, grid_size, subplot)
1239 # Reshape cube data into a single array to allow for a single histogram.
1240 # Otherwise we plot xdim histograms stacked.
1241 member_data_1d = (member.data).flatten()
1242 plt.hist(member_data_1d, density=True, stacked=True)
1243 ax = plt.gca()
1244 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1245 ax.set_xlim(vmin, vmax)
1247 # Overall figure title.
1248 fig.suptitle(title, fontsize=16)
1250 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1251 logging.info("Saved histogram postage stamp plot to %s", filename)
1252 plt.close(fig)
1255def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1256 cube: iris.cube.Cube,
1257 filename: str,
1258 title: str,
1259 stamp_coordinate: str,
1260 vmin: float,
1261 vmax: float,
1262 **kwargs,
1263):
1264 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1265 ax.set_title(title, fontsize=16)
1266 ax.set_xlim(vmin, vmax)
1267 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1268 ax.set_ylabel("normalised probability density", fontsize=14)
1269 # Loop over all slices along the stamp_coordinate
1270 for member in cube.slices_over(stamp_coordinate):
1271 # Flatten the member data to 1D
1272 member_data_1d = member.data.flatten()
1273 # Plot the histogram using plt.hist
1274 plt.hist(
1275 member_data_1d,
1276 density=True,
1277 stacked=True,
1278 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1279 )
1281 # Add a legend
1282 ax.legend(fontsize=16)
1284 # Save the figure to a file
1285 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1287 # Close the figure
1288 plt.close(fig)
1291def _plot_and_save_scattermap_plot(
1292 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1293):
1294 """Plot and save a geographical scatter plot.
1296 Parameters
1297 ----------
1298 cube: Cube
1299 1 dimensional Cube of the data points with auxiliary latitude and
1300 longitude coordinates,
1301 filename: str
1302 Filename of the plot to write.
1303 title: str
1304 Plot title.
1305 projection: str
1306 Mapping projection to be used by cartopy.
1307 """
1308 # Setup plot details, size, resolution, etc.
1309 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1310 if projection is not None:
1311 # Apart from the default, the only projection we currently support is
1312 # a stereographic projection over the North Pole.
1313 if projection == "NP_Stereo":
1314 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1315 else:
1316 raise ValueError(f"Unknown projection: {projection}")
1317 else:
1318 axes = plt.axes(projection=ccrs.PlateCarree())
1320 # Scatter plot of the field. The marker size is chosen to give
1321 # symbols that decrease in size as the number of observations
1322 # increases, although the fraction of the figure covered by
1323 # symbols increases roughly as N^(1/2), disregarding overlaps,
1324 # and has been selected for the default figure size of (10, 10).
1325 # Should this be changed, the marker size should be adjusted in
1326 # proportion to the area of the figure.
1327 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1328 klon = None
1329 klat = None
1330 for kc in range(len(cube.aux_coords)):
1331 if cube.aux_coords[kc].standard_name == "latitude":
1332 klat = kc
1333 elif cube.aux_coords[kc].standard_name == "longitude":
1334 klon = kc
1335 scatter_map = iplt.scatter(
1336 cube.aux_coords[klon],
1337 cube.aux_coords[klat],
1338 c=cube.data[:],
1339 s=mrk_size,
1340 cmap="jet",
1341 edgecolors="k",
1342 )
1344 # Add coastlines.
1345 try:
1346 axes.coastlines(resolution="10m")
1347 except AttributeError:
1348 pass
1350 # Add title.
1351 axes.set_title(title, fontsize=16)
1353 # Add colour bar.
1354 cbar = fig.colorbar(scatter_map)
1355 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1357 # Save plot.
1358 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1359 logging.info("Saved geographical scatter plot to %s", filename)
1360 plt.close(fig)
1363def _plot_and_save_power_spectrum_series(
1364 cubes: iris.cube.Cube | iris.cube.CubeList,
1365 filename: str,
1366 title: str,
1367 **kwargs,
1368):
1369 """Plot and save a power spectrum series.
1371 Parameters
1372 ----------
1373 cubes: Cube or CubeList
1374 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1375 filename: str
1376 Filename of the plot to write.
1377 title: str
1378 Plot title.
1379 """
1380 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1381 ax = plt.gca()
1383 model_colors_map = _get_model_colors_map(cubes)
1385 for cube in iter_maybe(cubes):
1386 # Calculate power spectrum
1388 # Extract time coordinate and convert to datetime
1389 time_coord = cube.coord("time")
1390 time_points = time_coord.units.num2date(time_coord.points)
1392 # Choose one time point (e.g., the first one)
1393 target_time = time_points[0]
1395 # Bind target_time inside the lambda using a default argument
1396 time_constraint = iris.Constraint(
1397 time=lambda cell, target_time=target_time: cell.point == target_time
1398 )
1400 cube = cube.extract(time_constraint)
1402 if cube.ndim == 2:
1403 cube_3d = cube.data[np.newaxis, :, :]
1404 logging.debug("Adding in new axis for a 2 dimensional cube.")
1405 elif cube.ndim == 3: 1405 ↛ 1406line 1405 didn't jump to line 1406 because the condition on line 1405 was never true
1406 cube_3d = cube.data
1407 else:
1408 raise ValueError("Cube dimensions unsuitable for power spectra code")
1409 raise ValueError(
1410 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1411 )
1413 # Calculate spectra
1414 ps_array = _DCT_ps(cube_3d)
1416 ps_cube = iris.cube.Cube(
1417 ps_array,
1418 long_name="power_spectra",
1419 )
1421 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1423 # Create a frequency/wavelength array for coordinate
1424 ps_len = ps_cube.data.shape[1]
1425 freqs = np.arange(1, ps_len + 1)
1426 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1428 # Convert datetime to numeric time using original units
1429 numeric_time = time_coord.units.date2num(time_points)
1430 # Create a new DimCoord with numeric time
1431 new_time_coord = iris.coords.DimCoord(
1432 numeric_time, standard_name="time", units=time_coord.units
1433 )
1435 # Add time and frequency coordinate to spectra cube.
1436 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1437 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1439 # Extract data from the cube
1440 frequency = ps_cube.coord("frequency").points
1441 power_spectrum = ps_cube.data
1443 label = None
1444 color = "black"
1445 if model_colors_map: 1445 ↛ 1446line 1445 didn't jump to line 1446 because the condition on line 1445 was never true
1446 label = ps_cube.attributes.get("model_name")
1447 color = model_colors_map[label]
1448 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1450 # Add some labels and tweak the style.
1451 ax.set_title(title, fontsize=16)
1452 ax.set_xlabel("Wavenumber", fontsize=14)
1453 ax.set_ylabel("Power", fontsize=14)
1454 ax.tick_params(axis="both", labelsize=12)
1456 # Set log-log scale
1457 ax.set_xscale("log")
1458 ax.set_yscale("log")
1460 # Overlay grid-lines onto power spectrum plot.
1461 ax.grid(linestyle="--", color="grey", linewidth=1)
1462 if model_colors_map: 1462 ↛ 1463line 1462 didn't jump to line 1463 because the condition on line 1462 was never true
1463 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1465 # Save plot.
1466 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1467 logging.info("Saved line plot to %s", filename)
1468 plt.close(fig)
1471def _plot_and_save_postage_stamp_power_spectrum_series(
1472 cube: iris.cube.Cube,
1473 filename: str,
1474 title: str,
1475 stamp_coordinate: str,
1476 **kwargs,
1477):
1478 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1480 Parameters
1481 ----------
1482 cube: Cube
1483 2 dimensional Cube of the data to plot as power spectrum.
1484 filename: str
1485 Filename of the plot to write.
1486 title: str
1487 Plot title.
1488 stamp_coordinate: str
1489 Coordinate that becomes different plots.
1490 """
1491 # Use the smallest square grid that will fit the members.
1492 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1494 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1495 # Make a subplot for each member.
1496 for member, subplot in zip(
1497 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1498 ):
1499 # Implicit interface is much easier here, due to needing to have the
1500 # cartopy GeoAxes generated.
1501 plt.subplot(grid_size, grid_size, subplot)
1503 frequency = member.coord("frequency").points
1504 power_spectrum = member.data
1506 ax = plt.gca()
1507 ax.plot(frequency, power_spectrum[0])
1508 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1510 # Overall figure title.
1511 fig.suptitle(title, fontsize=16)
1513 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1514 logging.info("Saved power spectra postage stamp plot to %s", filename)
1515 plt.close(fig)
1518def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1519 cube: iris.cube.Cube,
1520 filename: str,
1521 title: str,
1522 stamp_coordinate: str,
1523 **kwargs,
1524):
1525 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1526 ax.set_title(title, fontsize=16)
1527 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1528 ax.set_ylabel("Power", fontsize=14)
1529 # Loop over all slices along the stamp_coordinate
1530 for member in cube.slices_over(stamp_coordinate):
1531 frequency = member.coord("frequency").points
1532 power_spectrum = member.data
1533 ax.plot(
1534 frequency,
1535 power_spectrum[0],
1536 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1537 )
1539 # Add a legend
1540 ax.legend(fontsize=16)
1542 # Save the figure to a file
1543 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1545 # Close the figure
1546 plt.close(fig)
1549def _spatial_plot(
1550 method: Literal["contourf", "pcolormesh"],
1551 cube: iris.cube.Cube,
1552 filename: str | None,
1553 sequence_coordinate: str,
1554 stamp_coordinate: str,
1555 **kwargs,
1556):
1557 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1559 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1560 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1561 is present then postage stamp plots will be produced.
1563 Parameters
1564 ----------
1565 method: "contourf" | "pcolormesh"
1566 The plotting method to use.
1567 cube: Cube
1568 Iris cube of the data to plot. It should have two spatial dimensions,
1569 such as lat and lon, and may also have a another two dimension to be
1570 plotted sequentially and/or as postage stamp plots.
1571 filename: str | None
1572 Name of the plot to write, used as a prefix for plot sequences. If None
1573 uses the recipe name.
1574 sequence_coordinate: str
1575 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1576 This coordinate must exist in the cube.
1577 stamp_coordinate: str
1578 Coordinate about which to plot postage stamp plots. Defaults to
1579 ``"realization"``.
1581 Raises
1582 ------
1583 ValueError
1584 If the cube doesn't have the right dimensions.
1585 TypeError
1586 If the cube isn't a single cube.
1587 """
1588 recipe_title = get_recipe_metadata().get("title", "Untitled")
1590 # Ensure we have a name for the plot file.
1591 if filename is None:
1592 filename = slugify(recipe_title)
1594 # Ensure we've got a single cube.
1595 cube = _check_single_cube(cube)
1597 # Make postage stamp plots if stamp_coordinate exists and has more than a
1598 # single point.
1599 plotting_func = _plot_and_save_spatial_plot
1600 try:
1601 if cube.coord(stamp_coordinate).shape[0] > 1:
1602 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1603 except iris.exceptions.CoordinateNotFoundError:
1604 pass
1606 # Produce a geographical scatter plot if the data have a
1607 # dimension called observation or model_obs_error
1608 if any( 1608 ↛ 1612line 1608 didn't jump to line 1612 because the condition on line 1608 was never true
1609 crd.var_name == "station" or crd.var_name == "model_obs_error"
1610 for crd in cube.coords()
1611 ):
1612 plotting_func = _plot_and_save_scattermap_plot
1614 # Must have a sequence coordinate.
1615 try:
1616 cube.coord(sequence_coordinate)
1617 except iris.exceptions.CoordinateNotFoundError as err:
1618 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1620 # Create a plot for each value of the sequence coordinate.
1621 plot_index = []
1622 nplot = np.size(cube.coord(sequence_coordinate).points)
1623 for cube_slice in cube.slices_over(sequence_coordinate):
1624 # Use sequence value so multiple sequences can merge.
1625 sequence_value = cube_slice.coord(sequence_coordinate).points[0]
1626 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
1627 coord = cube_slice.coord(sequence_coordinate)
1628 # Format the coordinate value in a unit appropriate way.
1629 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
1630 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
1631 if nplot == 1 and coord.has_bounds:
1632 if np.size(coord.bounds) > 1:
1633 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
1634 # Do the actual plotting.
1635 plotting_func(
1636 cube_slice,
1637 filename=plot_filename,
1638 stamp_coordinate=stamp_coordinate,
1639 title=title,
1640 method=method,
1641 **kwargs,
1642 )
1643 plot_index.append(plot_filename)
1645 # Add list of plots to plot metadata.
1646 complete_plot_index = _append_to_plot_index(plot_index)
1648 # Make a page to display the plots.
1649 _make_plot_html_page(complete_plot_index)
1652def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1653 """Get colourmap for mask.
1655 If "mask_for_" appears anywhere in the name of a cube this function will be called
1656 regardless of the name of the variable to ensure a consistent plot.
1658 Parameters
1659 ----------
1660 cube: Cube
1661 Cube of variable for which the colorbar information is desired.
1662 axis: "x", "y", optional
1663 Select the levels for just this axis of a line plot. The min and max
1664 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1665 setting a universal range is not desirable (e.g. temperature), users
1666 can set ymin/ymax values to "auto" in the colorbar definitions file.
1667 Where no additional xmin/xmax or ymin/ymax values are provided, the
1668 axis bounds default to use the vmin/vmax values provided.
1670 Returns
1671 -------
1672 cmap:
1673 Matplotlib colormap.
1674 levels:
1675 List of levels to use for plotting. For continuous plots the min and max
1676 should be taken as the range.
1677 norm:
1678 BoundaryNorm information.
1679 """
1680 if "difference" not in cube.long_name:
1681 if axis:
1682 levels = [0, 1]
1683 # Complete settings based on levels.
1684 return None, levels, None
1685 else:
1686 # Define the levels and colors.
1687 levels = [0, 1, 2]
1688 colors = ["white", "dodgerblue"]
1689 # Create a custom color map.
1690 cmap = mcolors.ListedColormap(colors)
1691 # Normalize the levels.
1692 norm = mcolors.BoundaryNorm(levels, cmap.N)
1693 logging.debug("Colourmap for %s.", cube.long_name)
1694 return cmap, levels, norm
1695 else:
1696 if axis:
1697 levels = [-1, 1]
1698 return None, levels, None
1699 else:
1700 # Search for if mask difference, set to +/- 0.5 as values plotted <
1701 # not <=.
1702 levels = [-2, -0.5, 0.5, 2]
1703 colors = ["goldenrod", "white", "teal"]
1704 cmap = mcolors.ListedColormap(colors)
1705 norm = mcolors.BoundaryNorm(levels, cmap.N)
1706 logging.debug("Colourmap for %s.", cube.long_name)
1707 return cmap, levels, norm
1710def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1711 """Get a custom colorbar for a cube in the Beaufort Scale.
1713 Specific variable ranges can be separately set in user-supplied definition
1714 for x- or y-axis limits, or indicate where automated range preferred.
1716 Parameters
1717 ----------
1718 cube: Cube
1719 Cube of variable with Beaufort Scale in name.
1720 axis: "x", "y", optional
1721 Select the levels for just this axis of a line plot. The min and max
1722 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1723 setting a universal range is not desirable (e.g. temperature), users
1724 can set ymin/ymax values to "auto" in the colorbar definitions file.
1725 Where no additional xmin/xmax or ymin/ymax values are provided, the
1726 axis bounds default to use the vmin/vmax values provided.
1728 Returns
1729 -------
1730 cmap:
1731 Matplotlib colormap.
1732 levels:
1733 List of levels to use for plotting. For continuous plots the min and max
1734 should be taken as the range.
1735 norm:
1736 BoundaryNorm information.
1737 """
1738 if "difference" not in cube.long_name:
1739 if axis:
1740 levels = [0, 12]
1741 return None, levels, None
1742 else:
1743 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
1744 colors = [
1745 "black",
1746 (0, 0, 0.6),
1747 "blue",
1748 "cyan",
1749 "green",
1750 "yellow",
1751 (1, 0.5, 0),
1752 "red",
1753 "pink",
1754 "magenta",
1755 "purple",
1756 "maroon",
1757 "white",
1758 ]
1759 cmap = mcolors.ListedColormap(colors)
1760 norm = mcolors.BoundaryNorm(levels, cmap.N)
1761 logging.info("change colormap for Beaufort Scale colorbar.")
1762 return cmap, levels, norm
1763 else:
1764 if axis:
1765 levels = [-4, 4]
1766 return None, levels, None
1767 else:
1768 levels = [
1769 -3.5,
1770 -2.5,
1771 -1.5,
1772 -0.5,
1773 0.5,
1774 1.5,
1775 2.5,
1776 3.5,
1777 ]
1778 cmap = plt.get_cmap("bwr", 8)
1779 norm = mcolors.BoundaryNorm(levels, cmap.N)
1780 return cmap, levels, norm
1783def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
1784 """Return altered colourmap for temperature with change in units to Celsius.
1786 If "Celsius" appears anywhere in the name of a cube this function will be called.
1788 Parameters
1789 ----------
1790 cube: Cube
1791 Cube of variable for which the colorbar information is desired.
1792 cmap: Matplotlib colormap.
1793 levels: List
1794 List of levels to use for plotting. For continuous plots the min and max
1795 should be taken as the range.
1796 norm: BoundaryNorm.
1798 Returns
1799 -------
1800 cmap: Matplotlib colormap.
1801 levels: List
1802 List of levels to use for plotting. For continuous plots the min and max
1803 should be taken as the range.
1804 norm: BoundaryNorm.
1805 """
1806 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1807 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
1808 levels = np.array(levels)
1809 levels -= 273
1810 levels = levels.tolist()
1811 else:
1812 # Do nothing keep the existing colourbar attributes
1813 levels = levels
1814 cmap = cmap
1815 norm = norm
1816 return cmap, levels, norm
1819def _custom_colormap_probability(
1820 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
1821):
1822 """Get a custom colorbar for a probability cube.
1824 Specific variable ranges can be separately set in user-supplied definition
1825 for x- or y-axis limits, or indicate where automated range preferred.
1827 Parameters
1828 ----------
1829 cube: Cube
1830 Cube of variable with probability in name.
1831 axis: "x", "y", optional
1832 Select the levels for just this axis of a line plot. The min and max
1833 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1834 setting a universal range is not desirable (e.g. temperature), users
1835 can set ymin/ymax values to "auto" in the colorbar definitions file.
1836 Where no additional xmin/xmax or ymin/ymax values are provided, the
1837 axis bounds default to use the vmin/vmax values provided.
1839 Returns
1840 -------
1841 cmap:
1842 Matplotlib colormap.
1843 levels:
1844 List of levels to use for plotting. For continuous plots the min and max
1845 should be taken as the range.
1846 norm:
1847 BoundaryNorm information.
1848 """
1849 if axis:
1850 levels = [0, 1]
1851 return None, levels, None
1852 else:
1853 cmap = mcolors.ListedColormap(
1854 [
1855 "#FFFFFF",
1856 "#636363",
1857 "#e1dada",
1858 "#B5CAFF",
1859 "#8FB3FF",
1860 "#7F97FF",
1861 "#ABCF63",
1862 "#E8F59E",
1863 "#FFFA14",
1864 "#FFD121",
1865 "#FFA30A",
1866 ]
1867 )
1868 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]
1869 norm = mcolors.BoundaryNorm(levels, cmap.N)
1870 return cmap, levels, norm
1873def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
1874 """Return a custom colourmap for the current recipe."""
1875 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1876 if (
1877 any("surface_microphysical" in name for name in varnames)
1878 and "difference" not in cube.long_name
1879 and "mask" not in cube.long_name
1880 ):
1881 # Define the levels and colors
1882 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
1883 colors = [
1884 "w",
1885 (0, 0, 0.6),
1886 "b",
1887 "c",
1888 "g",
1889 "y",
1890 (1, 0.5, 0),
1891 "r",
1892 "pink",
1893 "m",
1894 "purple",
1895 "maroon",
1896 "gray",
1897 ]
1898 # Create a custom colormap
1899 cmap = mcolors.ListedColormap(colors)
1900 # Normalize the levels
1901 norm = mcolors.BoundaryNorm(levels, cmap.N)
1902 logging.info("change colormap for surface_microphysical variable colorbar.")
1903 else:
1904 # do nothing and keep existing colorbar attributes
1905 cmap = cmap
1906 levels = levels
1907 norm = norm
1908 return cmap, levels, norm
1911def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
1912 """Return custom colourmap for aviation colour state.
1914 If "aviation_colour_state" appears anywhere in the name of a cube
1915 this function will be called.
1917 Parameters
1918 ----------
1919 cube: Cube
1920 Cube of variable for which the colorbar information is desired.
1922 Returns
1923 -------
1924 cmap: Matplotlib colormap.
1925 levels: List
1926 List of levels to use for plotting. For continuous plots the min and max
1927 should be taken as the range.
1928 norm: BoundaryNorm.
1929 """
1930 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
1931 colors = [
1932 "#87ceeb",
1933 "#ffffff",
1934 "#8ced69",
1935 "#ffff00",
1936 "#ffd700",
1937 "#ffa500",
1938 "#fe3620",
1939 ]
1940 # Create a custom colormap
1941 cmap = mcolors.ListedColormap(colors)
1942 # Normalise the levels
1943 norm = mcolors.BoundaryNorm(levels, cmap.N)
1944 return cmap, levels, norm
1947def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
1948 """Return a custom colourmap for the current recipe."""
1949 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1950 if (
1951 any("visibility_in_air" in name for name in varnames)
1952 and "difference" not in cube.long_name
1953 and "mask" not in cube.long_name
1954 ):
1955 # Define the levels and colors (in km)
1956 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]
1957 norm = mcolors.BoundaryNorm(levels, cmap.N)
1958 colours = [
1959 "#8f00d6",
1960 "#d10000",
1961 "#ff9700",
1962 "#ffff00",
1963 "#00007f",
1964 "#6c9ccd",
1965 "#aae8ff",
1966 "#37a648",
1967 "#8edc64",
1968 "#c5ffc5",
1969 "#dcdcdc",
1970 "#ffffff",
1971 ]
1972 # Create a custom colormap
1973 cmap = mcolors.ListedColormap(colours)
1974 # Normalize the levels
1975 norm = mcolors.BoundaryNorm(levels, cmap.N)
1976 logging.info("change colormap for visibility_in_air variable colorbar.")
1977 else:
1978 # do nothing and keep existing colorbar attributes
1979 cmap = cmap
1980 levels = levels
1981 norm = norm
1982 return cmap, levels, norm
1985def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
1986 """Return number of models based on cube attributes."""
1987 model_names = list(
1988 filter(
1989 lambda x: x is not None,
1990 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
1991 )
1992 )
1993 if not model_names:
1994 logging.debug("Missing model names. Will assume single model.")
1995 return 1
1996 else:
1997 return len(model_names)
2000def _validate_cube_shape(
2001 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2002) -> None:
2003 """Check all cubes have a model name."""
2004 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2004 ↛ 2005line 2004 didn't jump to line 2005 because the condition on line 2004 was never true
2005 raise ValueError(
2006 f"The number of model names ({num_models}) should equal the number "
2007 f"of cubes ({len(cube)})."
2008 )
2011def _validate_cubes_coords(
2012 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2013) -> None:
2014 """Check same number of cubes as sequence coordinate for zip functions."""
2015 if len(cubes) != len(coords): 2015 ↛ 2016line 2015 didn't jump to line 2016 because the condition on line 2015 was never true
2016 raise ValueError(
2017 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2018 f"of sequence coordinates ({len(coords)})."
2019 f"Check that number of time entries in input data are consistent if "
2020 f"performing time-averaging steps prior to plotting outputs."
2021 )
2024def _calculate_CFAD(
2025 cube: iris.cube.Cube, vertical_coordinate: str, bin_edges: list[float]
2026) -> iris.cube.Cube:
2027 """Calculate a Contour Frequency by Altitude Diagram (CFAD).
2029 Parameters
2030 ----------
2031 cube: iris.cube.Cube
2032 A cube of the data to be turned into a CFAD. It should be a minimum
2033 of two dimensions with one being a user specified vertical coordinate.
2034 vertical_coordinate: str
2035 The vertical coordinate of the cube for the CFAD to be calculated over.
2036 bin_edges: list[float]
2037 The bin edges for the histogram. The bins need to be specified to
2038 ensure consistency across the CFAD, otherwise it cannot be interpreted.
2040 Notes
2041 -----
2042 Contour Frequency by Altitude Diagrams (CFADs) were first designed by
2043 Yuter and Houze (1995)[YuterandHouze95]. They are calculated by binning the
2044 data by altitude and then by variable bins (e.g. temperature). The variable
2045 bins are then normalised by each altitude. This essenitally creates a
2046 normalised frequency distribution for each altitude. These are then stacked
2047 and combined in a single plot.
2049 References
2050 ----------
2051 .. [YuterandHouze95] Yuter S.E., and Houze, R.A. (1995) "Three-Dimensional
2052 Kinematic and Microphysical Evolution of Florida Cumulonimbus. Part II:
2053 Frequency Distributions of Vertical Velocity, Reflectivity, and
2054 Differential Reflectivity" Monthly Weather Review, vol. 123, 1941-1963,
2055 doi: 10.1175/1520-0493(1995)123<1941:TDKAME>2.0.CO;2
2056 """
2057 # Setup empty array for containing the CFAD data.
2058 CFAD_values = np.zeros(
2059 (len(cube.coord(vertical_coordinate).points), len(bin_edges) - 1)
2060 )
2062 # Calculate the CFAD as a histogram summing to one for each level.
2063 for i, level_cube in enumerate(cube.slices_over(vertical_coordinate)):
2064 # Note setting density to True does not produce the correct
2065 # normalization for a CFAD, where each row must sum to one.
2066 CFAD_values[i, :] = (
2067 np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bin_edges)[
2068 0
2069 ]
2070 / level_cube.data.size
2071 )
2072 # Calculate central points for bins.
2073 bins = (np.array(bin_edges[:-1]) + np.array(bin_edges[1:])) / 2.0
2074 bin_bounds = np.array((bin_edges[:-1], bin_edges[1:])).T
2075 # Now construct the coordinates for the cube.
2076 vert_coord = cube.coord(vertical_coordinate)
2077 bin_coord = iris.coords.DimCoord(
2078 bins, bounds=bin_bounds, standard_name=cube.standard_name, units=cube.units
2079 )
2080 # Now construct the cube that is to be output.
2081 CFAD = iris.cube.Cube(
2082 CFAD_values,
2083 dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)],
2084 long_name=f"{cube.name()}_cfad",
2085 units="1",
2086 )
2087 CFAD.rename(f"{cube.name()}_cfad")
2088 return CFAD
2091####################
2092# Public functions #
2093####################
2096def spatial_contour_plot(
2097 cube: iris.cube.Cube,
2098 filename: str = None,
2099 sequence_coordinate: str = "time",
2100 stamp_coordinate: str = "realization",
2101 **kwargs,
2102) -> iris.cube.Cube:
2103 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2105 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2106 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2107 is present then postage stamp plots will be produced.
2109 Parameters
2110 ----------
2111 cube: Cube
2112 Iris cube of the data to plot. It should have two spatial dimensions,
2113 such as lat and lon, and may also have a another two dimension to be
2114 plotted sequentially and/or as postage stamp plots.
2115 filename: str, optional
2116 Name of the plot to write, used as a prefix for plot sequences. Defaults
2117 to the recipe name.
2118 sequence_coordinate: str, optional
2119 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2120 This coordinate must exist in the cube.
2121 stamp_coordinate: str, optional
2122 Coordinate about which to plot postage stamp plots. Defaults to
2123 ``"realization"``.
2125 Returns
2126 -------
2127 Cube
2128 The original cube (so further operations can be applied).
2130 Raises
2131 ------
2132 ValueError
2133 If the cube doesn't have the right dimensions.
2134 TypeError
2135 If the cube isn't a single cube.
2136 """
2137 _spatial_plot(
2138 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2139 )
2140 return cube
2143def spatial_pcolormesh_plot(
2144 cube: iris.cube.Cube,
2145 filename: str = None,
2146 sequence_coordinate: str = "time",
2147 stamp_coordinate: str = "realization",
2148 **kwargs,
2149) -> iris.cube.Cube:
2150 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2152 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2153 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2154 is present then postage stamp plots will be produced.
2156 This function is significantly faster than ``spatial_contour_plot``,
2157 especially at high resolutions, and should be preferred unless contiguous
2158 contour areas are important.
2160 Parameters
2161 ----------
2162 cube: Cube
2163 Iris cube of the data to plot. It should have two spatial dimensions,
2164 such as lat and lon, and may also have a another two dimension to be
2165 plotted sequentially and/or as postage stamp plots.
2166 filename: str, optional
2167 Name of the plot to write, used as a prefix for plot sequences. Defaults
2168 to the recipe name.
2169 sequence_coordinate: str, optional
2170 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2171 This coordinate must exist in the cube.
2172 stamp_coordinate: str, optional
2173 Coordinate about which to plot postage stamp plots. Defaults to
2174 ``"realization"``.
2176 Returns
2177 -------
2178 Cube
2179 The original cube (so further operations can be applied).
2181 Raises
2182 ------
2183 ValueError
2184 If the cube doesn't have the right dimensions.
2185 TypeError
2186 If the cube isn't a single cube.
2187 """
2188 _spatial_plot(
2189 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2190 )
2191 return cube
2194# TODO: Expand function to handle ensemble data.
2195# line_coordinate: str, optional
2196# Coordinate about which to plot multiple lines. Defaults to
2197# ``"realization"``.
2198def plot_line_series(
2199 cube: iris.cube.Cube | iris.cube.CubeList,
2200 filename: str = None,
2201 series_coordinate: str = "time",
2202 # line_coordinate: str = "realization",
2203 **kwargs,
2204) -> iris.cube.Cube | iris.cube.CubeList:
2205 """Plot a line plot for the specified coordinate.
2207 The Cube or CubeList must be 1D.
2209 Parameters
2210 ----------
2211 iris.cube | iris.cube.CubeList
2212 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2213 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2214 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2215 filename: str, optional
2216 Name of the plot to write, used as a prefix for plot sequences. Defaults
2217 to the recipe name.
2218 series_coordinate: str, optional
2219 Coordinate about which to make a series. Defaults to ``"time"``. This
2220 coordinate must exist in the cube.
2222 Returns
2223 -------
2224 iris.cube.Cube | iris.cube.CubeList
2225 The original Cube or CubeList (so further operations can be applied).
2226 plotted data.
2228 Raises
2229 ------
2230 ValueError
2231 If the cubes don't have the right dimensions.
2232 TypeError
2233 If the cube isn't a Cube or CubeList.
2234 """
2235 # Ensure we have a name for the plot file.
2236 title = get_recipe_metadata().get("title", "Untitled")
2238 if filename is None:
2239 filename = slugify(title)
2241 # Add file extension.
2242 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2244 num_models = _get_num_models(cube)
2246 _validate_cube_shape(cube, num_models)
2248 # Iterate over all cubes and extract coordinate to plot.
2249 cubes = iter_maybe(cube)
2250 coords = []
2251 for cube in cubes:
2252 try:
2253 coords.append(cube.coord(series_coordinate))
2254 except iris.exceptions.CoordinateNotFoundError as err:
2255 raise ValueError(
2256 f"Cube must have a {series_coordinate} coordinate."
2257 ) from err
2258 if cube.ndim > 2 or not cube.coords("realization"):
2259 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2261 # Do the actual plotting.
2262 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, title)
2264 # Add list of plots to plot metadata.
2265 plot_index = _append_to_plot_index([plot_filename])
2267 # Make a page to display the plots.
2268 _make_plot_html_page(plot_index)
2270 return cube
2273def plot_vertical_line_series(
2274 cubes: iris.cube.Cube | iris.cube.CubeList,
2275 filename: str = None,
2276 series_coordinate: str = "model_level_number",
2277 sequence_coordinate: str = "time",
2278 # line_coordinate: str = "realization",
2279 **kwargs,
2280) -> iris.cube.Cube | iris.cube.CubeList:
2281 """Plot a line plot against a type of vertical coordinate.
2283 The Cube or CubeList must be 1D.
2285 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2286 then a sequence of plots will be produced.
2288 Parameters
2289 ----------
2290 iris.cube | iris.cube.CubeList
2291 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2292 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2293 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2294 filename: str, optional
2295 Name of the plot to write, used as a prefix for plot sequences. Defaults
2296 to the recipe name.
2297 series_coordinate: str, optional
2298 Coordinate to plot on the y-axis. Can be ``pressure`` or
2299 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2300 for LFRic. Defaults to ``model_level_number``.
2301 This coordinate must exist in the cube.
2302 sequence_coordinate: str, optional
2303 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2304 This coordinate must exist in the cube.
2306 Returns
2307 -------
2308 iris.cube.Cube | iris.cube.CubeList
2309 The original Cube or CubeList (so further operations can be applied).
2310 Plotted data.
2312 Raises
2313 ------
2314 ValueError
2315 If the cubes doesn't have the right dimensions.
2316 TypeError
2317 If the cube isn't a Cube or CubeList.
2318 """
2319 # Ensure we have a name for the plot file.
2320 recipe_title = get_recipe_metadata().get("title", "Untitled")
2322 if filename is None:
2323 filename = slugify(recipe_title)
2325 cubes = iter_maybe(cubes)
2326 # Initialise empty list to hold all data from all cubes in a CubeList
2327 all_data = []
2329 # Store min/max ranges for x range.
2330 x_levels = []
2332 num_models = _get_num_models(cubes)
2334 _validate_cube_shape(cubes, num_models)
2336 # Iterate over all cubes in cube or CubeList and plot.
2337 coords = []
2338 for cube in cubes:
2339 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2340 try:
2341 coords.append(cube.coord(series_coordinate))
2342 except iris.exceptions.CoordinateNotFoundError as err:
2343 raise ValueError(
2344 f"Cube must have a {series_coordinate} coordinate."
2345 ) from err
2347 try:
2348 if cube.ndim > 1 or not cube.coords("realization"): 2348 ↛ 2356line 2348 didn't jump to line 2356 because the condition on line 2348 was always true
2349 cube.coord(sequence_coordinate)
2350 except iris.exceptions.CoordinateNotFoundError as err:
2351 raise ValueError(
2352 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2353 ) from err
2355 # Get minimum and maximum from levels information.
2356 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2357 if levels is not None: 2357 ↛ 2361line 2357 didn't jump to line 2361 because the condition on line 2357 was always true
2358 x_levels.append(min(levels))
2359 x_levels.append(max(levels))
2360 else:
2361 all_data.append(cube.data)
2363 if len(x_levels) == 0: 2363 ↛ 2365line 2363 didn't jump to line 2365 because the condition on line 2363 was never true
2364 # Combine all data into a single NumPy array
2365 combined_data = np.concatenate(all_data)
2367 # Set the lower and upper limit for the x-axis to ensure all plots have
2368 # same range. This needs to read the whole cube over the range of the
2369 # sequence and if applicable postage stamp coordinate.
2370 vmin = np.floor(combined_data.min())
2371 vmax = np.ceil(combined_data.max())
2372 else:
2373 vmin = min(x_levels)
2374 vmax = max(x_levels)
2376 # Matching the slices (matching by seq coord point; it may happen that
2377 # evaluated models do not cover the same seq coord range, hence matching
2378 # necessary)
2379 def filter_cube_iterables(cube_iterables) -> bool:
2380 return len(cube_iterables) == len(coords)
2382 cube_iterables = filter(
2383 filter_cube_iterables,
2384 (
2385 iris.cube.CubeList(
2386 s
2387 for s in itertools.chain.from_iterable(
2388 cb.slices_over(sequence_coordinate) for cb in cubes
2389 )
2390 if s.coord(sequence_coordinate).points[0] == point
2391 )
2392 for point in sorted(
2393 set(
2394 itertools.chain.from_iterable(
2395 cb.coord(sequence_coordinate).points for cb in cubes
2396 )
2397 )
2398 )
2399 ),
2400 )
2402 # Create a plot for each value of the sequence coordinate.
2403 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2404 # similar sequence values. Passing a CubeList into the internal plotting function
2405 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2406 # or an iris.cube.CubeList.
2407 plot_index = []
2408 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2409 for cubes_slice in cube_iterables:
2410 # Use sequence value so multiple sequences can merge.
2411 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2412 sequence_value = seq_coord.points[0]
2413 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2414 # Format the coordinate value in a unit appropriate way.
2415 title = f"{recipe_title}\n [{seq_coord.units.title(sequence_value)}]"
2416 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2417 if nplot == 1 and seq_coord.has_bounds: 2417 ↛ 2418line 2417 didn't jump to line 2418 because the condition on line 2417 was never true
2418 if np.size(seq_coord.bounds) > 1:
2419 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2420 # Do the actual plotting.
2421 _plot_and_save_vertical_line_series(
2422 cubes_slice,
2423 coords,
2424 "realization",
2425 plot_filename,
2426 series_coordinate,
2427 title=title,
2428 vmin=vmin,
2429 vmax=vmax,
2430 )
2431 plot_index.append(plot_filename)
2433 # Add list of plots to plot metadata.
2434 complete_plot_index = _append_to_plot_index(plot_index)
2436 # Make a page to display the plots.
2437 _make_plot_html_page(complete_plot_index)
2439 return cubes
2442def qq_plot(
2443 cubes: iris.cube.CubeList,
2444 coordinates: list[str],
2445 percentiles: list[float],
2446 model_names: list[str],
2447 filename: str = None,
2448 one_to_one: bool = True,
2449 **kwargs,
2450) -> iris.cube.CubeList:
2451 """Plot a Quantile-Quantile plot between two models for common time points.
2453 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2454 collapsed within the operator over all specified coordinates such as
2455 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2456 ensemble members to ensure a 1D cube (array).
2458 Parameters
2459 ----------
2460 cubes: iris.cube.CubeList
2461 Two cubes of the same variable with different models.
2462 coordinate: list[str]
2463 The list of coordinates to collapse over. This list should be
2464 every coordinate within the cube to result in a 1D cube around
2465 the percentile coordinate.
2466 percent: list[float]
2467 A list of percentiles to appear in the plot.
2468 model_names: list[str]
2469 A list of model names to appear on the axis of the plot.
2470 filename: str, optional
2471 Filename of the plot to write.
2472 one_to_one: bool, optional
2473 If True a 1:1 line is plotted; if False it is not. Default is True.
2475 Raises
2476 ------
2477 ValueError
2478 When the cubes are not compatible.
2480 Notes
2481 -----
2482 The quantile-quantile plot is a variant on the scatter plot representing
2483 two datasets by their quantiles (percentiles) for common time points.
2484 This plot does not use a theoretical distribution to compare against, but
2485 compares percentiles of two datasets. This plot does
2486 not use all raw data points, but plots the selected percentiles (quantiles) of
2487 each variable instead for the two datasets, thereby normalising the data for a
2488 direct comparison between the selected percentiles of the two dataset distributions.
2490 Quantile-quantile plots are valuable for comparing against
2491 observations and other models. Identical percentiles between the variables
2492 will lie on the one-to-one line implying the values correspond well to each
2493 other. Where there is a deviation from the one-to-one line a range of
2494 possibilities exist depending on how and where the data is shifted (e.g.,
2495 Wilks 2011 [Wilks2011]_).
2497 For distributions above the one-to-one line the distribution is left-skewed;
2498 below is right-skewed. A distinct break implies a bimodal distribution, and
2499 closer values/values further apart at the tails imply poor representation of
2500 the extremes.
2502 References
2503 ----------
2504 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2505 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2506 """
2507 # Check cubes using same functionality as the difference operator.
2508 if len(cubes) != 2:
2509 raise ValueError("cubes should contain exactly 2 cubes.")
2510 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2511 other: Cube = cubes.extract_cube(
2512 iris.Constraint(
2513 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2514 )
2515 )
2517 # Get spatial coord names.
2518 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2519 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2521 # Ensure cubes to compare are on common differencing grid.
2522 # This is triggered if either
2523 # i) latitude and longitude shapes are not the same. Note grid points
2524 # are not compared directly as these can differ through rounding
2525 # errors.
2526 # ii) or variables are known to often sit on different grid staggering
2527 # in different models (e.g. cell center vs cell edge), as is the case
2528 # for UM and LFRic comparisons.
2529 # In future greater choice of regridding method might be applied depending
2530 # on variable type. Linear regridding can in general be appropriate for smooth
2531 # variables. Care should be taken with interpretation of differences
2532 # given this dependency on regridding.
2533 if (
2534 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2535 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2536 ) or (
2537 base.long_name
2538 in [
2539 "eastward_wind_at_10m",
2540 "northward_wind_at_10m",
2541 "northward_wind_at_cell_centres",
2542 "eastward_wind_at_cell_centres",
2543 "zonal_wind_at_pressure_levels",
2544 "meridional_wind_at_pressure_levels",
2545 "potential_vorticity_at_pressure_levels",
2546 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2547 ]
2548 ):
2549 logging.debug(
2550 "Linear regridding base cube to other grid to compute differences"
2551 )
2552 base = regrid_onto_cube(base, other, method="Linear")
2554 # Extract just common time points.
2555 base, other = _extract_common_time_points(base, other)
2557 # Equalise attributes so we can merge.
2558 fully_equalise_attributes([base, other])
2559 logging.debug("Base: %s\nOther: %s", base, other)
2561 # Collapse cubes.
2562 base = collapse(
2563 base,
2564 coordinate=coordinates,
2565 method="PERCENTILE",
2566 additional_percent=percentiles,
2567 )
2568 other = collapse(
2569 other,
2570 coordinate=coordinates,
2571 method="PERCENTILE",
2572 additional_percent=percentiles,
2573 )
2575 # Ensure we have a name for the plot file.
2576 title = get_recipe_metadata().get("title", "Untitled")
2578 if filename is None:
2579 filename = slugify(title)
2581 # Add file extension.
2582 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2584 # Do the actual plotting on a scatter plot
2585 _plot_and_save_scatter_plot(
2586 base, other, plot_filename, title, one_to_one, model_names
2587 )
2589 # Add list of plots to plot metadata.
2590 plot_index = _append_to_plot_index([plot_filename])
2592 # Make a page to display the plots.
2593 _make_plot_html_page(plot_index)
2595 return iris.cube.CubeList([base, other])
2598def scatter_plot(
2599 cube_x: iris.cube.Cube | iris.cube.CubeList,
2600 cube_y: iris.cube.Cube | iris.cube.CubeList,
2601 filename: str = None,
2602 one_to_one: bool = True,
2603 **kwargs,
2604) -> iris.cube.CubeList:
2605 """Plot a scatter plot between two variables.
2607 Both cubes must be 1D.
2609 Parameters
2610 ----------
2611 cube_x: Cube | CubeList
2612 1 dimensional Cube of the data to plot on y-axis.
2613 cube_y: Cube | CubeList
2614 1 dimensional Cube of the data to plot on x-axis.
2615 filename: str, optional
2616 Filename of the plot to write.
2617 one_to_one: bool, optional
2618 If True a 1:1 line is plotted; if False it is not. Default is True.
2620 Returns
2621 -------
2622 cubes: CubeList
2623 CubeList of the original x and y cubes for further processing.
2625 Raises
2626 ------
2627 ValueError
2628 If the cube doesn't have the right dimensions and cubes not the same
2629 size.
2630 TypeError
2631 If the cube isn't a single cube.
2633 Notes
2634 -----
2635 Scatter plots are used for determining if there is a relationship between
2636 two variables. Positive relations have a slope going from bottom left to top
2637 right; Negative relations have a slope going from top left to bottom right.
2638 """
2639 # Iterate over all cubes in cube or CubeList and plot.
2640 for cube_iter in iter_maybe(cube_x):
2641 # Check cubes are correct shape.
2642 cube_iter = _check_single_cube(cube_iter)
2643 if cube_iter.ndim > 1:
2644 raise ValueError("cube_x must be 1D.")
2646 # Iterate over all cubes in cube or CubeList and plot.
2647 for cube_iter in iter_maybe(cube_y):
2648 # Check cubes are correct shape.
2649 cube_iter = _check_single_cube(cube_iter)
2650 if cube_iter.ndim > 1:
2651 raise ValueError("cube_y must be 1D.")
2653 # Ensure we have a name for the plot file.
2654 title = get_recipe_metadata().get("title", "Untitled")
2656 if filename is None:
2657 filename = slugify(title)
2659 # Add file extension.
2660 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2662 # Do the actual plotting.
2663 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2665 # Add list of plots to plot metadata.
2666 plot_index = _append_to_plot_index([plot_filename])
2668 # Make a page to display the plots.
2669 _make_plot_html_page(plot_index)
2671 return iris.cube.CubeList([cube_x, cube_y])
2674def vector_plot(
2675 cube_u: iris.cube.Cube,
2676 cube_v: iris.cube.Cube,
2677 filename: str = None,
2678 sequence_coordinate: str = "time",
2679 **kwargs,
2680) -> iris.cube.CubeList:
2681 """Plot a vector plot based on the input u and v components."""
2682 recipe_title = get_recipe_metadata().get("title", "Untitled")
2684 # Ensure we have a name for the plot file.
2685 if filename is None: 2685 ↛ 2686line 2685 didn't jump to line 2686 because the condition on line 2685 was never true
2686 filename = slugify(recipe_title)
2688 # Cubes must have a matching sequence coordinate.
2689 try:
2690 # Check that the u and v cubes have the same sequence coordinate.
2691 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2691 ↛ 2692line 2691 didn't jump to line 2692 because the condition on line 2691 was never true
2692 raise ValueError("Coordinates do not match.")
2693 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2694 raise ValueError(
2695 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2696 ) from err
2698 # Create a plot for each value of the sequence coordinate.
2699 plot_index = []
2700 for cube_u_slice, cube_v_slice in zip(
2701 cube_u.slices_over(sequence_coordinate),
2702 cube_v.slices_over(sequence_coordinate),
2703 strict=True,
2704 ):
2705 # Use sequence value so multiple sequences can merge.
2706 sequence_value = cube_u_slice.coord(sequence_coordinate).points[0]
2707 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2708 coord = cube_u_slice.coord(sequence_coordinate)
2709 # Format the coordinate value in a unit appropriate way.
2710 title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
2711 # Do the actual plotting.
2712 _plot_and_save_vector_plot(
2713 cube_u_slice,
2714 cube_v_slice,
2715 filename=plot_filename,
2716 title=title,
2717 method="contourf",
2718 )
2719 plot_index.append(plot_filename)
2721 # Add list of plots to plot metadata.
2722 complete_plot_index = _append_to_plot_index(plot_index)
2724 # Make a page to display the plots.
2725 _make_plot_html_page(complete_plot_index)
2727 return iris.cube.CubeList([cube_u, cube_v])
2730def plot_histogram_series(
2731 cubes: iris.cube.Cube | iris.cube.CubeList,
2732 filename: str = None,
2733 sequence_coordinate: str = "time",
2734 stamp_coordinate: str = "realization",
2735 single_plot: bool = False,
2736 **kwargs,
2737) -> iris.cube.Cube | iris.cube.CubeList:
2738 """Plot a histogram plot for each vertical level provided.
2740 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2741 is present then a sequence of plots will be produced using the time slider
2742 functionality to scroll through histograms against time. If a
2743 stamp_coordinate is present then postage stamp plots will be produced. If
2744 stamp_coordinate and single_plot is True, all postage stamp plots will be
2745 plotted in a single plot instead of separate postage stamp plots.
2747 Parameters
2748 ----------
2749 cubes: Cube | iris.cube.CubeList
2750 Iris cube or CubeList of the data to plot. It should have a single dimension other
2751 than the stamp coordinate.
2752 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2753 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2754 filename: str, optional
2755 Name of the plot to write, used as a prefix for plot sequences. Defaults
2756 to the recipe name.
2757 sequence_coordinate: str, optional
2758 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2759 This coordinate must exist in the cube and will be used for the time
2760 slider.
2761 stamp_coordinate: str, optional
2762 Coordinate about which to plot postage stamp plots. Defaults to
2763 ``"realization"``.
2764 single_plot: bool, optional
2765 If True, all postage stamp plots will be plotted in a single plot. If
2766 False, each postage stamp plot will be plotted separately. Is only valid
2767 if stamp_coordinate exists and has more than a single point.
2769 Returns
2770 -------
2771 iris.cube.Cube | iris.cube.CubeList
2772 The original Cube or CubeList (so further operations can be applied).
2773 Plotted data.
2775 Raises
2776 ------
2777 ValueError
2778 If the cube doesn't have the right dimensions.
2779 TypeError
2780 If the cube isn't a Cube or CubeList.
2781 """
2782 recipe_title = get_recipe_metadata().get("title", "Untitled")
2784 cubes = iter_maybe(cubes)
2786 # Ensure we have a name for the plot file.
2787 if filename is None:
2788 filename = slugify(recipe_title)
2790 # Internal plotting function.
2791 plotting_func = _plot_and_save_histogram_series
2793 num_models = _get_num_models(cubes)
2795 _validate_cube_shape(cubes, num_models)
2797 # If several histograms are plotted with time as sequence_coordinate for the
2798 # time slider option.
2799 for cube in cubes:
2800 try:
2801 cube.coord(sequence_coordinate)
2802 except iris.exceptions.CoordinateNotFoundError as err:
2803 raise ValueError(
2804 f"Cube must have a {sequence_coordinate} coordinate."
2805 ) from err
2807 # Get minimum and maximum from levels information.
2808 levels = None
2809 for cube in cubes: 2809 ↛ 2825line 2809 didn't jump to line 2825 because the loop on line 2809 didn't complete
2810 # First check if user-specified "auto" range variable.
2811 # This maintains the value of levels as None, so proceed.
2812 _, levels, _ = _colorbar_map_levels(cube, axis="y")
2813 if levels is None:
2814 break
2815 # If levels is changed, recheck to use the vmin,vmax or
2816 # levels-based ranges for histogram plots.
2817 _, levels, _ = _colorbar_map_levels(cube)
2818 logging.debug("levels: %s", levels)
2819 if levels is not None: 2819 ↛ 2809line 2819 didn't jump to line 2809 because the condition on line 2819 was always true
2820 vmin = min(levels)
2821 vmax = max(levels)
2822 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
2823 break
2825 if levels is None:
2826 vmin = min(cb.data.min() for cb in cubes)
2827 vmax = max(cb.data.max() for cb in cubes)
2829 # Make postage stamp plots if stamp_coordinate exists and has more than a
2830 # single point. If single_plot is True:
2831 # -- all postage stamp plots will be plotted in a single plot instead of
2832 # separate postage stamp plots.
2833 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2834 # produced per single model only
2835 if num_models == 1: 2835 ↛ 2848line 2835 didn't jump to line 2848 because the condition on line 2835 was always true
2836 if ( 2836 ↛ 2840line 2836 didn't jump to line 2840 because the condition on line 2836 was never true
2837 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2838 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2839 ):
2840 if single_plot:
2841 plotting_func = (
2842 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2843 )
2844 else:
2845 plotting_func = _plot_and_save_postage_stamp_histogram_series
2846 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2847 else:
2848 all_points = sorted(
2849 set(
2850 itertools.chain.from_iterable(
2851 cb.coord(sequence_coordinate).points for cb in cubes
2852 )
2853 )
2854 )
2855 all_slices = list(
2856 itertools.chain.from_iterable(
2857 cb.slices_over(sequence_coordinate) for cb in cubes
2858 )
2859 )
2860 # Matched slices (matched by seq coord point; it may happen that
2861 # evaluated models do not cover the same seq coord range, hence matching
2862 # necessary)
2863 cube_iterables = [
2864 iris.cube.CubeList(
2865 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2866 )
2867 for point in all_points
2868 ]
2870 plot_index = []
2871 nplot = np.size(cube.coord(sequence_coordinate).points)
2872 # Create a plot for each value of the sequence coordinate. Allowing for
2873 # multiple cubes in a CubeList to be plotted in the same plot for similar
2874 # sequence values. Passing a CubeList into the internal plotting function
2875 # for similar values of the sequence coordinate. cube_slice can be an
2876 # iris.cube.Cube or an iris.cube.CubeList.
2877 for cube_slice in cube_iterables:
2878 single_cube = cube_slice
2879 if isinstance(cube_slice, iris.cube.CubeList): 2879 ↛ 2880line 2879 didn't jump to line 2880 because the condition on line 2879 was never true
2880 single_cube = cube_slice[0]
2882 # Use sequence value so multiple sequences can merge.
2883 sequence_value = single_cube.coord(sequence_coordinate).points[0]
2884 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2885 coord = single_cube.coord(sequence_coordinate)
2886 # Format the coordinate value in a unit appropriate way.
2887 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
2888 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2889 if nplot == 1 and coord.has_bounds: 2889 ↛ 2890line 2889 didn't jump to line 2890 because the condition on line 2889 was never true
2890 if np.size(coord.bounds) > 1:
2891 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
2892 # Do the actual plotting.
2893 plotting_func(
2894 cube_slice,
2895 filename=plot_filename,
2896 stamp_coordinate=stamp_coordinate,
2897 title=title,
2898 vmin=vmin,
2899 vmax=vmax,
2900 )
2901 plot_index.append(plot_filename)
2903 # Add list of plots to plot metadata.
2904 complete_plot_index = _append_to_plot_index(plot_index)
2906 # Make a page to display the plots.
2907 _make_plot_html_page(complete_plot_index)
2909 return cubes
2912def plot_power_spectrum_series(
2913 cubes: iris.cube.Cube | iris.cube.CubeList,
2914 filename: str = None,
2915 sequence_coordinate: str = "time",
2916 stamp_coordinate: str = "realization",
2917 single_plot: bool = False,
2918 **kwargs,
2919) -> iris.cube.Cube | iris.cube.CubeList:
2920 """Plot a power spectrum plot for each vertical level provided.
2922 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
2923 is present then a sequence of plots will be produced using the time slider
2924 functionality to scroll through power spectrum against time. If a
2925 stamp_coordinate is present then postage stamp plots will be produced. If
2926 stamp_coordinate and single_plot is True, all postage stamp plots will be
2927 plotted in a single plot instead of separate postage stamp plots.
2929 Parameters
2930 ----------
2931 cubes: Cube | iris.cube.CubeList
2932 Iris cube or CubeList of the data to plot. It should have a single dimension other
2933 than the stamp coordinate.
2934 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2935 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2936 filename: str, optional
2937 Name of the plot to write, used as a prefix for plot sequences. Defaults
2938 to the recipe name.
2939 sequence_coordinate: str, optional
2940 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2941 This coordinate must exist in the cube and will be used for the time
2942 slider.
2943 stamp_coordinate: str, optional
2944 Coordinate about which to plot postage stamp plots. Defaults to
2945 ``"realization"``.
2946 single_plot: bool, optional
2947 If True, all postage stamp plots will be plotted in a single plot. If
2948 False, each postage stamp plot will be plotted separately. Is only valid
2949 if stamp_coordinate exists and has more than a single point.
2951 Returns
2952 -------
2953 iris.cube.Cube | iris.cube.CubeList
2954 The original Cube or CubeList (so further operations can be applied).
2955 Plotted data.
2957 Raises
2958 ------
2959 ValueError
2960 If the cube doesn't have the right dimensions.
2961 TypeError
2962 If the cube isn't a Cube or CubeList.
2963 """
2964 recipe_title = get_recipe_metadata().get("title", "Untitled")
2966 cubes = iter_maybe(cubes)
2967 # Ensure we have a name for the plot file.
2968 if filename is None:
2969 filename = slugify(recipe_title)
2971 # Internal plotting function.
2972 plotting_func = _plot_and_save_power_spectrum_series
2974 num_models = _get_num_models(cubes)
2976 _validate_cube_shape(cubes, num_models)
2978 # If several power spectra are plotted with time as sequence_coordinate for the
2979 # time slider option.
2980 for cube in cubes:
2981 try:
2982 cube.coord(sequence_coordinate)
2983 except iris.exceptions.CoordinateNotFoundError as err:
2984 raise ValueError(
2985 f"Cube must have a {sequence_coordinate} coordinate."
2986 ) from err
2988 # Make postage stamp plots if stamp_coordinate exists and has more than a
2989 # single point. If single_plot is True:
2990 # -- all postage stamp plots will be plotted in a single plot instead of
2991 # separate postage stamp plots.
2992 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2993 # produced per single model only
2994 if num_models == 1: 2994 ↛ 3007line 2994 didn't jump to line 3007 because the condition on line 2994 was always true
2995 if ( 2995 ↛ 2999line 2995 didn't jump to line 2999 because the condition on line 2995 was never true
2996 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2997 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2998 ):
2999 if single_plot:
3000 plotting_func = (
3001 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3002 )
3003 else:
3004 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3005 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3006 else:
3007 all_points = sorted(
3008 set(
3009 itertools.chain.from_iterable(
3010 cb.coord(sequence_coordinate).points for cb in cubes
3011 )
3012 )
3013 )
3014 all_slices = list(
3015 itertools.chain.from_iterable(
3016 cb.slices_over(sequence_coordinate) for cb in cubes
3017 )
3018 )
3019 # Matched slices (matched by seq coord point; it may happen that
3020 # evaluated models do not cover the same seq coord range, hence matching
3021 # necessary)
3022 cube_iterables = [
3023 iris.cube.CubeList(
3024 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3025 )
3026 for point in all_points
3027 ]
3029 plot_index = []
3030 nplot = np.size(cube.coord(sequence_coordinate).points)
3031 # Create a plot for each value of the sequence coordinate. Allowing for
3032 # multiple cubes in a CubeList to be plotted in the same plot for similar
3033 # sequence values. Passing a CubeList into the internal plotting function
3034 # for similar values of the sequence coordinate. cube_slice can be an
3035 # iris.cube.Cube or an iris.cube.CubeList.
3036 for cube_slice in cube_iterables:
3037 single_cube = cube_slice
3038 if isinstance(cube_slice, iris.cube.CubeList): 3038 ↛ 3039line 3038 didn't jump to line 3039 because the condition on line 3038 was never true
3039 single_cube = cube_slice[0]
3041 # Use sequence value so multiple sequences can merge.
3042 sequence_value = single_cube.coord(sequence_coordinate).points[0]
3043 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
3044 coord = single_cube.coord(sequence_coordinate)
3045 # Format the coordinate value in a unit appropriate way.
3046 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
3047 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
3048 if nplot == 1 and coord.has_bounds: 3048 ↛ 3052line 3048 didn't jump to line 3052 because the condition on line 3048 was always true
3049 if np.size(coord.bounds) > 1:
3050 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
3051 # Do the actual plotting.
3052 plotting_func(
3053 cube_slice,
3054 filename=plot_filename,
3055 stamp_coordinate=stamp_coordinate,
3056 title=title,
3057 )
3058 plot_index.append(plot_filename)
3060 # Add list of plots to plot metadata.
3061 complete_plot_index = _append_to_plot_index(plot_index)
3063 # Make a page to display the plots.
3064 _make_plot_html_page(complete_plot_index)
3066 return cubes
3069def _DCT_ps(y_3d):
3070 """Calculate power spectra for regional domains.
3072 Parameters
3073 ----------
3074 y_3d: 3D array
3075 3 dimensional array to calculate spectrum for.
3076 (2D field data with 3rd dimension of time)
3078 Returns: ps_array
3079 Array of power spectra values calculated for input field (for each time)
3081 Method for regional domains:
3082 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3083 as described in Denis et al 2002 [Denis_etal_2002]_.
3085 References
3086 ----------
3087 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3088 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3089 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3090 Monthly Weather Review, Vol. 130, 1812-1828
3091 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3092 """
3093 Nt, Ny, Nx = y_3d.shape
3095 # Max coefficient
3096 Nmin = min(Nx - 1, Ny - 1)
3098 # Create alpha matrix (of wavenumbers)
3099 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3101 # Prepare output array
3102 ps_array = np.zeros((Nt, Nmin))
3104 # Loop over time to get spectrum for each time.
3105 for t in range(Nt):
3106 y_2d = y_3d[t]
3108 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3109 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3110 # cosine basis functions at different spatial frequencies.
3112 # normalise spectrum to allow comparison between models.
3113 fkk = fft.dctn(y_2d, norm="ortho")
3115 # Normalise fkk
3116 fkk = fkk / np.sqrt(Ny * Nx)
3118 # calculate variance of spectral coefficient
3119 sigma_2 = fkk**2 / Nx / Ny
3121 # Group ellipses of alphas into the same wavenumber k/Nmin
3122 for k in range(1, Nmin + 1):
3123 alpha = k / Nmin
3124 alpha_p1 = (k + 1) / Nmin
3126 # Sum up elements matching k
3127 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3128 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3130 return ps_array
3133def _create_alpha_matrix(Ny, Nx):
3134 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3136 Parameters
3137 ----------
3138 Ny, Nx:
3139 Dimensions of the 2D field for which the power spectra is calculated. Used to
3140 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3141 single-scale parameter.
3143 Returns: alpha_matrix
3144 normalisation of 2D wavenumber axes, transforming the spectral domain into
3145 an elliptic coordinate system.
3147 """
3148 # Create x_indices: each row is [1, 2, ..., Nx]
3149 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3151 # Create y_indices: each column is [1, 2, ..., Ny]
3152 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3154 # Compute alpha_matrix
3155 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3157 return alpha_matrix