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