Coverage for src / CSET / operators / plot.py: 86%
961 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-24 08:36 +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 )
2006def _calculate_CFAD(
2007 cube: iris.cube.Cube, vertical_coordinate: str, bin_edges: list[float]
2008) -> iris.cube.Cube:
2009 """Calculate a Contour Frequency by Altitude Diagram (CFAD).
2011 Parameters
2012 ----------
2013 cube: iris.cube.Cube
2014 A cube of the data to be turned into a CFAD. It should be a minimum
2015 of two dimensions with one being a user specified vertical coordinate.
2016 vertical_coordinate: str
2017 The vertical coordinate of the cube for the CFAD to be calculated over.
2018 bin_edges: list[float]
2019 The bin edges for the histogram. The bins need to be specified to
2020 ensure consistency across the CFAD, otherwise it cannot be interpreted.
2022 Notes
2023 -----
2024 Contour Frequency by Altitude Diagrams (CFADs) were first designed by
2025 Yuter and Houze (1995)[YuterandHouze95]. They are calculated by binning the
2026 data by altitude and then by variable bins (e.g. temperature). The variable
2027 bins are then normalised by each altitude. This essenitally creates a
2028 normalised frequency distribution for each altitude. These are then stacked
2029 and combined in a single plot.
2031 References
2032 ----------
2033 .. [YuterandHouze95] Yuter S.E., and Houze, R.A. (1995) "Three-Dimensional
2034 Kinematic and Microphysical Evolution of Florida Cumulonimbus. Part II:
2035 Frequency Distributions of Vertical Velocity, Reflectivity, and
2036 Differential Reflectivity" Monthly Weather Review, vol. 123, 1941-1963,
2037 doi: 10.1175/1520-0493(1995)123<1941:TDKAME>2.0.CO;2
2038 """
2039 # Setup empty array for containing the CFAD data.
2040 CFAD_values = np.zeros(
2041 (len(cube.coord(vertical_coordinate).points), len(bin_edges) - 1)
2042 )
2044 # Calculate the CFAD as a histogram summing to one for each level.
2045 for i, level_cube in enumerate(cube.slices_over(vertical_coordinate)):
2046 # Note setting density to True does not produce the correct
2047 # normalization for a CFAD, where each row must sum to one.
2048 CFAD_values[i, :] = (
2049 np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bin_edges)[
2050 0
2051 ]
2052 / level_cube.data.size
2053 )
2054 # Calculate central points for bins.
2055 bins = (np.array(bin_edges[:-1]) + np.array(bin_edges[1:])) / 2.0
2056 bin_bounds = np.array((bin_edges[:-1], bin_edges[1:])).T
2057 # Now construct the coordinates for the cube.
2058 vert_coord = cube.coord(vertical_coordinate)
2059 bin_coord = iris.coords.DimCoord(
2060 bins, bounds=bin_bounds, standard_name=cube.standard_name, units=cube.units
2061 )
2062 # Now construct the cube that is to be output.
2063 CFAD = iris.cube.Cube(
2064 CFAD_values,
2065 dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)],
2066 long_name=f"{cube.name()}_cfad",
2067 units="1",
2068 )
2069 CFAD.rename(f"{cube.name()}_cfad")
2070 return CFAD
2073####################
2074# Public functions #
2075####################
2078def spatial_contour_plot(
2079 cube: iris.cube.Cube,
2080 filename: str = None,
2081 sequence_coordinate: str = "time",
2082 stamp_coordinate: str = "realization",
2083 **kwargs,
2084) -> iris.cube.Cube:
2085 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2087 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2088 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2089 is present then postage stamp plots will be produced.
2091 Parameters
2092 ----------
2093 cube: Cube
2094 Iris cube of the data to plot. It should have two spatial dimensions,
2095 such as lat and lon, and may also have a another two dimension to be
2096 plotted sequentially and/or as postage stamp plots.
2097 filename: str, optional
2098 Name of the plot to write, used as a prefix for plot sequences. Defaults
2099 to the recipe name.
2100 sequence_coordinate: str, optional
2101 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2102 This coordinate must exist in the cube.
2103 stamp_coordinate: str, optional
2104 Coordinate about which to plot postage stamp plots. Defaults to
2105 ``"realization"``.
2107 Returns
2108 -------
2109 Cube
2110 The original cube (so further operations can be applied).
2112 Raises
2113 ------
2114 ValueError
2115 If the cube doesn't have the right dimensions.
2116 TypeError
2117 If the cube isn't a single cube.
2118 """
2119 _spatial_plot(
2120 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2121 )
2122 return cube
2125def spatial_pcolormesh_plot(
2126 cube: iris.cube.Cube,
2127 filename: str = None,
2128 sequence_coordinate: str = "time",
2129 stamp_coordinate: str = "realization",
2130 **kwargs,
2131) -> iris.cube.Cube:
2132 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2134 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2135 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2136 is present then postage stamp plots will be produced.
2138 This function is significantly faster than ``spatial_contour_plot``,
2139 especially at high resolutions, and should be preferred unless contiguous
2140 contour areas are important.
2142 Parameters
2143 ----------
2144 cube: Cube
2145 Iris cube of the data to plot. It should have two spatial dimensions,
2146 such as lat and lon, and may also have a another two dimension to be
2147 plotted sequentially and/or as postage stamp plots.
2148 filename: str, optional
2149 Name of the plot to write, used as a prefix for plot sequences. Defaults
2150 to the recipe name.
2151 sequence_coordinate: str, optional
2152 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2153 This coordinate must exist in the cube.
2154 stamp_coordinate: str, optional
2155 Coordinate about which to plot postage stamp plots. Defaults to
2156 ``"realization"``.
2158 Returns
2159 -------
2160 Cube
2161 The original cube (so further operations can be applied).
2163 Raises
2164 ------
2165 ValueError
2166 If the cube doesn't have the right dimensions.
2167 TypeError
2168 If the cube isn't a single cube.
2169 """
2170 _spatial_plot(
2171 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2172 )
2173 return cube
2176# TODO: Expand function to handle ensemble data.
2177# line_coordinate: str, optional
2178# Coordinate about which to plot multiple lines. Defaults to
2179# ``"realization"``.
2180def plot_line_series(
2181 cube: iris.cube.Cube | iris.cube.CubeList,
2182 filename: str = None,
2183 series_coordinate: str = "time",
2184 # line_coordinate: str = "realization",
2185 **kwargs,
2186) -> iris.cube.Cube | iris.cube.CubeList:
2187 """Plot a line plot for the specified coordinate.
2189 The Cube or CubeList must be 1D.
2191 Parameters
2192 ----------
2193 iris.cube | iris.cube.CubeList
2194 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2195 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2196 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2197 filename: str, optional
2198 Name of the plot to write, used as a prefix for plot sequences. Defaults
2199 to the recipe name.
2200 series_coordinate: str, optional
2201 Coordinate about which to make a series. Defaults to ``"time"``. This
2202 coordinate must exist in the cube.
2204 Returns
2205 -------
2206 iris.cube.Cube | iris.cube.CubeList
2207 The original Cube or CubeList (so further operations can be applied).
2208 plotted data.
2210 Raises
2211 ------
2212 ValueError
2213 If the cubes don't have the right dimensions.
2214 TypeError
2215 If the cube isn't a Cube or CubeList.
2216 """
2217 # Ensure we have a name for the plot file.
2218 title = get_recipe_metadata().get("title", "Untitled")
2220 if filename is None:
2221 filename = slugify(title)
2223 # Add file extension.
2224 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2226 num_models = _get_num_models(cube)
2228 _validate_cube_shape(cube, num_models)
2230 # Iterate over all cubes and extract coordinate to plot.
2231 cubes = iter_maybe(cube)
2232 coords = []
2233 for cube in cubes:
2234 try:
2235 coords.append(cube.coord(series_coordinate))
2236 except iris.exceptions.CoordinateNotFoundError as err:
2237 raise ValueError(
2238 f"Cube must have a {series_coordinate} coordinate."
2239 ) from err
2240 if cube.ndim > 2 or not cube.coords("realization"):
2241 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2243 # Do the actual plotting.
2244 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, title)
2246 # Add list of plots to plot metadata.
2247 plot_index = _append_to_plot_index([plot_filename])
2249 # Make a page to display the plots.
2250 _make_plot_html_page(plot_index)
2252 return cube
2255def plot_vertical_line_series(
2256 cubes: iris.cube.Cube | iris.cube.CubeList,
2257 filename: str = None,
2258 series_coordinate: str = "model_level_number",
2259 sequence_coordinate: str = "time",
2260 # line_coordinate: str = "realization",
2261 **kwargs,
2262) -> iris.cube.Cube | iris.cube.CubeList:
2263 """Plot a line plot against a type of vertical coordinate.
2265 The Cube or CubeList must be 1D.
2267 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2268 then a sequence of plots will be produced.
2270 Parameters
2271 ----------
2272 iris.cube | iris.cube.CubeList
2273 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2274 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2275 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2276 filename: str, optional
2277 Name of the plot to write, used as a prefix for plot sequences. Defaults
2278 to the recipe name.
2279 series_coordinate: str, optional
2280 Coordinate to plot on the y-axis. Can be ``pressure`` or
2281 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2282 for LFRic. Defaults to ``model_level_number``.
2283 This coordinate must exist in the cube.
2284 sequence_coordinate: str, optional
2285 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2286 This coordinate must exist in the cube.
2288 Returns
2289 -------
2290 iris.cube.Cube | iris.cube.CubeList
2291 The original Cube or CubeList (so further operations can be applied).
2292 Plotted data.
2294 Raises
2295 ------
2296 ValueError
2297 If the cubes doesn't have the right dimensions.
2298 TypeError
2299 If the cube isn't a Cube or CubeList.
2300 """
2301 # Ensure we have a name for the plot file.
2302 recipe_title = get_recipe_metadata().get("title", "Untitled")
2304 if filename is None:
2305 filename = slugify(recipe_title)
2307 cubes = iter_maybe(cubes)
2308 # Initialise empty list to hold all data from all cubes in a CubeList
2309 all_data = []
2311 # Store min/max ranges for x range.
2312 x_levels = []
2314 num_models = _get_num_models(cubes)
2316 _validate_cube_shape(cubes, num_models)
2318 # Iterate over all cubes in cube or CubeList and plot.
2319 coords = []
2320 for cube in cubes:
2321 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2322 try:
2323 coords.append(cube.coord(series_coordinate))
2324 except iris.exceptions.CoordinateNotFoundError as err:
2325 raise ValueError(
2326 f"Cube must have a {series_coordinate} coordinate."
2327 ) from err
2329 try:
2330 if cube.ndim > 1 or not cube.coords("realization"): 2330 ↛ 2338line 2330 didn't jump to line 2338 because the condition on line 2330 was always true
2331 cube.coord(sequence_coordinate)
2332 except iris.exceptions.CoordinateNotFoundError as err:
2333 raise ValueError(
2334 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2335 ) from err
2337 # Get minimum and maximum from levels information.
2338 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2339 if levels is not None: 2339 ↛ 2343line 2339 didn't jump to line 2343 because the condition on line 2339 was always true
2340 x_levels.append(min(levels))
2341 x_levels.append(max(levels))
2342 else:
2343 all_data.append(cube.data)
2345 if len(x_levels) == 0: 2345 ↛ 2347line 2345 didn't jump to line 2347 because the condition on line 2345 was never true
2346 # Combine all data into a single NumPy array
2347 combined_data = np.concatenate(all_data)
2349 # Set the lower and upper limit for the x-axis to ensure all plots have
2350 # same range. This needs to read the whole cube over the range of the
2351 # sequence and if applicable postage stamp coordinate.
2352 vmin = np.floor(combined_data.min())
2353 vmax = np.ceil(combined_data.max())
2354 else:
2355 vmin = min(x_levels)
2356 vmax = max(x_levels)
2358 # Matching the slices (matching by seq coord point; it may happen that
2359 # evaluated models do not cover the same seq coord range, hence matching
2360 # necessary)
2361 def filter_cube_iterables(cube_iterables) -> bool:
2362 return len(cube_iterables) == len(coords)
2364 cube_iterables = filter(
2365 filter_cube_iterables,
2366 (
2367 iris.cube.CubeList(
2368 s
2369 for s in itertools.chain.from_iterable(
2370 cb.slices_over(sequence_coordinate) for cb in cubes
2371 )
2372 if s.coord(sequence_coordinate).points[0] == point
2373 )
2374 for point in sorted(
2375 set(
2376 itertools.chain.from_iterable(
2377 cb.coord(sequence_coordinate).points for cb in cubes
2378 )
2379 )
2380 )
2381 ),
2382 )
2384 # Create a plot for each value of the sequence coordinate.
2385 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2386 # similar sequence values. Passing a CubeList into the internal plotting function
2387 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2388 # or an iris.cube.CubeList.
2389 plot_index = []
2390 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2391 for cubes_slice in cube_iterables:
2392 # Use sequence value so multiple sequences can merge.
2393 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2394 sequence_value = seq_coord.points[0]
2395 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2396 # Format the coordinate value in a unit appropriate way.
2397 title = f"{recipe_title}\n [{seq_coord.units.title(sequence_value)}]"
2398 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2399 if nplot == 1 and seq_coord.has_bounds: 2399 ↛ 2400line 2399 didn't jump to line 2400 because the condition on line 2399 was never true
2400 if np.size(seq_coord.bounds) > 1:
2401 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2402 # Do the actual plotting.
2403 _plot_and_save_vertical_line_series(
2404 cubes_slice,
2405 coords,
2406 "realization",
2407 plot_filename,
2408 series_coordinate,
2409 title=title,
2410 vmin=vmin,
2411 vmax=vmax,
2412 )
2413 plot_index.append(plot_filename)
2415 # Add list of plots to plot metadata.
2416 complete_plot_index = _append_to_plot_index(plot_index)
2418 # Make a page to display the plots.
2419 _make_plot_html_page(complete_plot_index)
2421 return cubes
2424def qq_plot(
2425 cubes: iris.cube.CubeList,
2426 coordinates: list[str],
2427 percentiles: list[float],
2428 model_names: list[str],
2429 filename: str = None,
2430 one_to_one: bool = True,
2431 **kwargs,
2432) -> iris.cube.CubeList:
2433 """Plot a Quantile-Quantile plot between two models for common time points.
2435 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2436 collapsed within the operator over all specified coordinates such as
2437 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2438 ensemble members to ensure a 1D cube (array).
2440 Parameters
2441 ----------
2442 cubes: iris.cube.CubeList
2443 Two cubes of the same variable with different models.
2444 coordinate: list[str]
2445 The list of coordinates to collapse over. This list should be
2446 every coordinate within the cube to result in a 1D cube around
2447 the percentile coordinate.
2448 percent: list[float]
2449 A list of percentiles to appear in the plot.
2450 model_names: list[str]
2451 A list of model names to appear on the axis of the plot.
2452 filename: str, optional
2453 Filename of the plot to write.
2454 one_to_one: bool, optional
2455 If True a 1:1 line is plotted; if False it is not. Default is True.
2457 Raises
2458 ------
2459 ValueError
2460 When the cubes are not compatible.
2462 Notes
2463 -----
2464 The quantile-quantile plot is a variant on the scatter plot representing
2465 two datasets by their quantiles (percentiles) for common time points.
2466 This plot does not use a theoretical distribution to compare against, but
2467 compares percentiles of two datasets. This plot does
2468 not use all raw data points, but plots the selected percentiles (quantiles) of
2469 each variable instead for the two datasets, thereby normalising the data for a
2470 direct comparison between the selected percentiles of the two dataset distributions.
2472 Quantile-quantile plots are valuable for comparing against
2473 observations and other models. Identical percentiles between the variables
2474 will lie on the one-to-one line implying the values correspond well to each
2475 other. Where there is a deviation from the one-to-one line a range of
2476 possibilities exist depending on how and where the data is shifted (e.g.,
2477 Wilks 2011 [Wilks2011]_).
2479 For distributions above the one-to-one line the distribution is left-skewed;
2480 below is right-skewed. A distinct break implies a bimodal distribution, and
2481 closer values/values further apart at the tails imply poor representation of
2482 the extremes.
2484 References
2485 ----------
2486 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2487 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2488 """
2489 # Check cubes using same functionality as the difference operator.
2490 if len(cubes) != 2:
2491 raise ValueError("cubes should contain exactly 2 cubes.")
2492 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2493 other: Cube = cubes.extract_cube(
2494 iris.Constraint(
2495 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2496 )
2497 )
2499 # Get spatial coord names.
2500 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2501 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2503 # Ensure cubes to compare are on common differencing grid.
2504 # This is triggered if either
2505 # i) latitude and longitude shapes are not the same. Note grid points
2506 # are not compared directly as these can differ through rounding
2507 # errors.
2508 # ii) or variables are known to often sit on different grid staggering
2509 # in different models (e.g. cell center vs cell edge), as is the case
2510 # for UM and LFRic comparisons.
2511 # In future greater choice of regridding method might be applied depending
2512 # on variable type. Linear regridding can in general be appropriate for smooth
2513 # variables. Care should be taken with interpretation of differences
2514 # given this dependency on regridding.
2515 if (
2516 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2517 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2518 ) or (
2519 base.long_name
2520 in [
2521 "eastward_wind_at_10m",
2522 "northward_wind_at_10m",
2523 "northward_wind_at_cell_centres",
2524 "eastward_wind_at_cell_centres",
2525 "zonal_wind_at_pressure_levels",
2526 "meridional_wind_at_pressure_levels",
2527 "potential_vorticity_at_pressure_levels",
2528 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2529 ]
2530 ):
2531 logging.debug(
2532 "Linear regridding base cube to other grid to compute differences"
2533 )
2534 base = regrid_onto_cube(base, other, method="Linear")
2536 # Extract just common time points.
2537 base, other = _extract_common_time_points(base, other)
2539 # Equalise attributes so we can merge.
2540 fully_equalise_attributes([base, other])
2541 logging.debug("Base: %s\nOther: %s", base, other)
2543 # Collapse cubes.
2544 base = collapse(
2545 base,
2546 coordinate=coordinates,
2547 method="PERCENTILE",
2548 additional_percent=percentiles,
2549 )
2550 other = collapse(
2551 other,
2552 coordinate=coordinates,
2553 method="PERCENTILE",
2554 additional_percent=percentiles,
2555 )
2557 # Ensure we have a name for the plot file.
2558 title = get_recipe_metadata().get("title", "Untitled")
2560 if filename is None:
2561 filename = slugify(title)
2563 # Add file extension.
2564 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2566 # Do the actual plotting on a scatter plot
2567 _plot_and_save_scatter_plot(
2568 base, other, plot_filename, title, one_to_one, model_names
2569 )
2571 # Add list of plots to plot metadata.
2572 plot_index = _append_to_plot_index([plot_filename])
2574 # Make a page to display the plots.
2575 _make_plot_html_page(plot_index)
2577 return iris.cube.CubeList([base, other])
2580def scatter_plot(
2581 cube_x: iris.cube.Cube | iris.cube.CubeList,
2582 cube_y: iris.cube.Cube | iris.cube.CubeList,
2583 filename: str = None,
2584 one_to_one: bool = True,
2585 **kwargs,
2586) -> iris.cube.CubeList:
2587 """Plot a scatter plot between two variables.
2589 Both cubes must be 1D.
2591 Parameters
2592 ----------
2593 cube_x: Cube | CubeList
2594 1 dimensional Cube of the data to plot on y-axis.
2595 cube_y: Cube | CubeList
2596 1 dimensional Cube of the data to plot on x-axis.
2597 filename: str, optional
2598 Filename of the plot to write.
2599 one_to_one: bool, optional
2600 If True a 1:1 line is plotted; if False it is not. Default is True.
2602 Returns
2603 -------
2604 cubes: CubeList
2605 CubeList of the original x and y cubes for further processing.
2607 Raises
2608 ------
2609 ValueError
2610 If the cube doesn't have the right dimensions and cubes not the same
2611 size.
2612 TypeError
2613 If the cube isn't a single cube.
2615 Notes
2616 -----
2617 Scatter plots are used for determining if there is a relationship between
2618 two variables. Positive relations have a slope going from bottom left to top
2619 right; Negative relations have a slope going from top left to bottom right.
2620 """
2621 # Iterate over all cubes in cube or CubeList and plot.
2622 for cube_iter in iter_maybe(cube_x):
2623 # Check cubes are correct shape.
2624 cube_iter = _check_single_cube(cube_iter)
2625 if cube_iter.ndim > 1:
2626 raise ValueError("cube_x must be 1D.")
2628 # Iterate over all cubes in cube or CubeList and plot.
2629 for cube_iter in iter_maybe(cube_y):
2630 # Check cubes are correct shape.
2631 cube_iter = _check_single_cube(cube_iter)
2632 if cube_iter.ndim > 1:
2633 raise ValueError("cube_y must be 1D.")
2635 # Ensure we have a name for the plot file.
2636 title = get_recipe_metadata().get("title", "Untitled")
2638 if filename is None:
2639 filename = slugify(title)
2641 # Add file extension.
2642 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2644 # Do the actual plotting.
2645 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2647 # Add list of plots to plot metadata.
2648 plot_index = _append_to_plot_index([plot_filename])
2650 # Make a page to display the plots.
2651 _make_plot_html_page(plot_index)
2653 return iris.cube.CubeList([cube_x, cube_y])
2656def vector_plot(
2657 cube_u: iris.cube.Cube,
2658 cube_v: iris.cube.Cube,
2659 filename: str = None,
2660 sequence_coordinate: str = "time",
2661 **kwargs,
2662) -> iris.cube.CubeList:
2663 """Plot a vector plot based on the input u and v components."""
2664 recipe_title = get_recipe_metadata().get("title", "Untitled")
2666 # Ensure we have a name for the plot file.
2667 if filename is None: 2667 ↛ 2668line 2667 didn't jump to line 2668 because the condition on line 2667 was never true
2668 filename = slugify(recipe_title)
2670 # Cubes must have a matching sequence coordinate.
2671 try:
2672 # Check that the u and v cubes have the same sequence coordinate.
2673 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2673 ↛ 2674line 2673 didn't jump to line 2674 because the condition on line 2673 was never true
2674 raise ValueError("Coordinates do not match.")
2675 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2676 raise ValueError(
2677 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2678 ) from err
2680 # Create a plot for each value of the sequence coordinate.
2681 plot_index = []
2682 for cube_u_slice, cube_v_slice in zip(
2683 cube_u.slices_over(sequence_coordinate),
2684 cube_v.slices_over(sequence_coordinate),
2685 strict=True,
2686 ):
2687 # Use sequence value so multiple sequences can merge.
2688 sequence_value = cube_u_slice.coord(sequence_coordinate).points[0]
2689 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2690 coord = cube_u_slice.coord(sequence_coordinate)
2691 # Format the coordinate value in a unit appropriate way.
2692 title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
2693 # Do the actual plotting.
2694 _plot_and_save_vector_plot(
2695 cube_u_slice,
2696 cube_v_slice,
2697 filename=plot_filename,
2698 title=title,
2699 method="contourf",
2700 )
2701 plot_index.append(plot_filename)
2703 # Add list of plots to plot metadata.
2704 complete_plot_index = _append_to_plot_index(plot_index)
2706 # Make a page to display the plots.
2707 _make_plot_html_page(complete_plot_index)
2709 return iris.cube.CubeList([cube_u, cube_v])
2712def plot_histogram_series(
2713 cubes: iris.cube.Cube | iris.cube.CubeList,
2714 filename: str = None,
2715 sequence_coordinate: str = "time",
2716 stamp_coordinate: str = "realization",
2717 single_plot: bool = False,
2718 **kwargs,
2719) -> iris.cube.Cube | iris.cube.CubeList:
2720 """Plot a histogram plot for each vertical level provided.
2722 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2723 is present then a sequence of plots will be produced using the time slider
2724 functionality to scroll through histograms against time. If a
2725 stamp_coordinate is present then postage stamp plots will be produced. If
2726 stamp_coordinate and single_plot is True, all postage stamp plots will be
2727 plotted in a single plot instead of separate postage stamp plots.
2729 Parameters
2730 ----------
2731 cubes: Cube | iris.cube.CubeList
2732 Iris cube or CubeList of the data to plot. It should have a single dimension other
2733 than the stamp coordinate.
2734 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2735 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2736 filename: str, optional
2737 Name of the plot to write, used as a prefix for plot sequences. Defaults
2738 to the recipe name.
2739 sequence_coordinate: str, optional
2740 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2741 This coordinate must exist in the cube and will be used for the time
2742 slider.
2743 stamp_coordinate: str, optional
2744 Coordinate about which to plot postage stamp plots. Defaults to
2745 ``"realization"``.
2746 single_plot: bool, optional
2747 If True, all postage stamp plots will be plotted in a single plot. If
2748 False, each postage stamp plot will be plotted separately. Is only valid
2749 if stamp_coordinate exists and has more than a single point.
2751 Returns
2752 -------
2753 iris.cube.Cube | iris.cube.CubeList
2754 The original Cube or CubeList (so further operations can be applied).
2755 Plotted data.
2757 Raises
2758 ------
2759 ValueError
2760 If the cube doesn't have the right dimensions.
2761 TypeError
2762 If the cube isn't a Cube or CubeList.
2763 """
2764 recipe_title = get_recipe_metadata().get("title", "Untitled")
2766 cubes = iter_maybe(cubes)
2768 # Ensure we have a name for the plot file.
2769 if filename is None:
2770 filename = slugify(recipe_title)
2772 # Internal plotting function.
2773 plotting_func = _plot_and_save_histogram_series
2775 num_models = _get_num_models(cubes)
2777 _validate_cube_shape(cubes, num_models)
2779 # If several histograms are plotted with time as sequence_coordinate for the
2780 # time slider option.
2781 for cube in cubes:
2782 try:
2783 cube.coord(sequence_coordinate)
2784 except iris.exceptions.CoordinateNotFoundError as err:
2785 raise ValueError(
2786 f"Cube must have a {sequence_coordinate} coordinate."
2787 ) from err
2789 # Get minimum and maximum from levels information.
2790 levels = None
2791 for cube in cubes: 2791 ↛ 2807line 2791 didn't jump to line 2807 because the loop on line 2791 didn't complete
2792 # First check if user-specified "auto" range variable.
2793 # This maintains the value of levels as None, so proceed.
2794 _, levels, _ = _colorbar_map_levels(cube, axis="y")
2795 if levels is None:
2796 break
2797 # If levels is changed, recheck to use the vmin,vmax or
2798 # levels-based ranges for histogram plots.
2799 _, levels, _ = _colorbar_map_levels(cube)
2800 logging.debug("levels: %s", levels)
2801 if levels is not None: 2801 ↛ 2791line 2801 didn't jump to line 2791 because the condition on line 2801 was always true
2802 vmin = min(levels)
2803 vmax = max(levels)
2804 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
2805 break
2807 if levels is None:
2808 vmin = min(cb.data.min() for cb in cubes)
2809 vmax = max(cb.data.max() for cb in cubes)
2811 # Make postage stamp plots if stamp_coordinate exists and has more than a
2812 # single point. If single_plot is True:
2813 # -- all postage stamp plots will be plotted in a single plot instead of
2814 # separate postage stamp plots.
2815 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2816 # produced per single model only
2817 if num_models == 1: 2817 ↛ 2830line 2817 didn't jump to line 2830 because the condition on line 2817 was always true
2818 if ( 2818 ↛ 2822line 2818 didn't jump to line 2822 because the condition on line 2818 was never true
2819 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2820 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2821 ):
2822 if single_plot:
2823 plotting_func = (
2824 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2825 )
2826 else:
2827 plotting_func = _plot_and_save_postage_stamp_histogram_series
2828 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2829 else:
2830 all_points = sorted(
2831 set(
2832 itertools.chain.from_iterable(
2833 cb.coord(sequence_coordinate).points for cb in cubes
2834 )
2835 )
2836 )
2837 all_slices = list(
2838 itertools.chain.from_iterable(
2839 cb.slices_over(sequence_coordinate) for cb in cubes
2840 )
2841 )
2842 # Matched slices (matched by seq coord point; it may happen that
2843 # evaluated models do not cover the same seq coord range, hence matching
2844 # necessary)
2845 cube_iterables = [
2846 iris.cube.CubeList(
2847 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2848 )
2849 for point in all_points
2850 ]
2852 plot_index = []
2853 nplot = np.size(cube.coord(sequence_coordinate).points)
2854 # Create a plot for each value of the sequence coordinate. Allowing for
2855 # multiple cubes in a CubeList to be plotted in the same plot for similar
2856 # sequence values. Passing a CubeList into the internal plotting function
2857 # for similar values of the sequence coordinate. cube_slice can be an
2858 # iris.cube.Cube or an iris.cube.CubeList.
2859 for cube_slice in cube_iterables:
2860 single_cube = cube_slice
2861 if isinstance(cube_slice, iris.cube.CubeList): 2861 ↛ 2862line 2861 didn't jump to line 2862 because the condition on line 2861 was never true
2862 single_cube = cube_slice[0]
2864 # Use sequence value so multiple sequences can merge.
2865 sequence_value = single_cube.coord(sequence_coordinate).points[0]
2866 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2867 coord = single_cube.coord(sequence_coordinate)
2868 # Format the coordinate value in a unit appropriate way.
2869 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
2870 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2871 if nplot == 1 and coord.has_bounds: 2871 ↛ 2872line 2871 didn't jump to line 2872 because the condition on line 2871 was never true
2872 if np.size(coord.bounds) > 1:
2873 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
2874 # Do the actual plotting.
2875 plotting_func(
2876 cube_slice,
2877 filename=plot_filename,
2878 stamp_coordinate=stamp_coordinate,
2879 title=title,
2880 vmin=vmin,
2881 vmax=vmax,
2882 )
2883 plot_index.append(plot_filename)
2885 # Add list of plots to plot metadata.
2886 complete_plot_index = _append_to_plot_index(plot_index)
2888 # Make a page to display the plots.
2889 _make_plot_html_page(complete_plot_index)
2891 return cubes
2894def plot_power_spectrum_series(
2895 cubes: iris.cube.Cube | iris.cube.CubeList,
2896 filename: str = None,
2897 sequence_coordinate: str = "time",
2898 stamp_coordinate: str = "realization",
2899 single_plot: bool = False,
2900 **kwargs,
2901) -> iris.cube.Cube | iris.cube.CubeList:
2902 """Plot a power spectrum plot for each vertical level provided.
2904 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
2905 is present then a sequence of plots will be produced using the time slider
2906 functionality to scroll through power spectrum against time. If a
2907 stamp_coordinate is present then postage stamp plots will be produced. If
2908 stamp_coordinate and single_plot is True, all postage stamp plots will be
2909 plotted in a single plot instead of separate postage stamp plots.
2911 Parameters
2912 ----------
2913 cubes: Cube | iris.cube.CubeList
2914 Iris cube or CubeList of the data to plot. It should have a single dimension other
2915 than the stamp coordinate.
2916 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2917 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2918 filename: str, optional
2919 Name of the plot to write, used as a prefix for plot sequences. Defaults
2920 to the recipe name.
2921 sequence_coordinate: str, optional
2922 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2923 This coordinate must exist in the cube and will be used for the time
2924 slider.
2925 stamp_coordinate: str, optional
2926 Coordinate about which to plot postage stamp plots. Defaults to
2927 ``"realization"``.
2928 single_plot: bool, optional
2929 If True, all postage stamp plots will be plotted in a single plot. If
2930 False, each postage stamp plot will be plotted separately. Is only valid
2931 if stamp_coordinate exists and has more than a single point.
2933 Returns
2934 -------
2935 iris.cube.Cube | iris.cube.CubeList
2936 The original Cube or CubeList (so further operations can be applied).
2937 Plotted data.
2939 Raises
2940 ------
2941 ValueError
2942 If the cube doesn't have the right dimensions.
2943 TypeError
2944 If the cube isn't a Cube or CubeList.
2945 """
2946 recipe_title = get_recipe_metadata().get("title", "Untitled")
2948 cubes = iter_maybe(cubes)
2949 # Ensure we have a name for the plot file.
2950 if filename is None:
2951 filename = slugify(recipe_title)
2953 # Internal plotting function.
2954 plotting_func = _plot_and_save_power_spectrum_series
2956 num_models = _get_num_models(cubes)
2958 _validate_cube_shape(cubes, num_models)
2960 # If several power spectra are plotted with time as sequence_coordinate for the
2961 # time slider option.
2962 for cube in cubes:
2963 try:
2964 cube.coord(sequence_coordinate)
2965 except iris.exceptions.CoordinateNotFoundError as err:
2966 raise ValueError(
2967 f"Cube must have a {sequence_coordinate} coordinate."
2968 ) from err
2970 # Make postage stamp plots if stamp_coordinate exists and has more than a
2971 # single point. If single_plot is True:
2972 # -- all postage stamp plots will be plotted in a single plot instead of
2973 # separate postage stamp plots.
2974 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2975 # produced per single model only
2976 if num_models == 1: 2976 ↛ 2989line 2976 didn't jump to line 2989 because the condition on line 2976 was always true
2977 if ( 2977 ↛ 2981line 2977 didn't jump to line 2981 because the condition on line 2977 was never true
2978 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2979 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2980 ):
2981 if single_plot:
2982 plotting_func = (
2983 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2984 )
2985 else:
2986 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2987 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2988 else:
2989 all_points = sorted(
2990 set(
2991 itertools.chain.from_iterable(
2992 cb.coord(sequence_coordinate).points for cb in cubes
2993 )
2994 )
2995 )
2996 all_slices = list(
2997 itertools.chain.from_iterable(
2998 cb.slices_over(sequence_coordinate) for cb in cubes
2999 )
3000 )
3001 # Matched slices (matched by seq coord point; it may happen that
3002 # evaluated models do not cover the same seq coord range, hence matching
3003 # necessary)
3004 cube_iterables = [
3005 iris.cube.CubeList(
3006 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3007 )
3008 for point in all_points
3009 ]
3011 plot_index = []
3012 nplot = np.size(cube.coord(sequence_coordinate).points)
3013 # Create a plot for each value of the sequence coordinate. Allowing for
3014 # multiple cubes in a CubeList to be plotted in the same plot for similar
3015 # sequence values. Passing a CubeList into the internal plotting function
3016 # for similar values of the sequence coordinate. cube_slice can be an
3017 # iris.cube.Cube or an iris.cube.CubeList.
3018 for cube_slice in cube_iterables:
3019 single_cube = cube_slice
3020 if isinstance(cube_slice, iris.cube.CubeList): 3020 ↛ 3021line 3020 didn't jump to line 3021 because the condition on line 3020 was never true
3021 single_cube = cube_slice[0]
3023 # Use sequence value so multiple sequences can merge.
3024 sequence_value = single_cube.coord(sequence_coordinate).points[0]
3025 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
3026 coord = single_cube.coord(sequence_coordinate)
3027 # Format the coordinate value in a unit appropriate way.
3028 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
3029 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
3030 if nplot == 1 and coord.has_bounds: 3030 ↛ 3034line 3030 didn't jump to line 3034 because the condition on line 3030 was always true
3031 if np.size(coord.bounds) > 1:
3032 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
3033 # Do the actual plotting.
3034 plotting_func(
3035 cube_slice,
3036 filename=plot_filename,
3037 stamp_coordinate=stamp_coordinate,
3038 title=title,
3039 )
3040 plot_index.append(plot_filename)
3042 # Add list of plots to plot metadata.
3043 complete_plot_index = _append_to_plot_index(plot_index)
3045 # Make a page to display the plots.
3046 _make_plot_html_page(complete_plot_index)
3048 return cubes
3051def _DCT_ps(y_3d):
3052 """Calculate power spectra for regional domains.
3054 Parameters
3055 ----------
3056 y_3d: 3D array
3057 3 dimensional array to calculate spectrum for.
3058 (2D field data with 3rd dimension of time)
3060 Returns: ps_array
3061 Array of power spectra values calculated for input field (for each time)
3063 Method for regional domains:
3064 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3065 as described in Denis et al 2002 [Denis_etal_2002]_.
3067 References
3068 ----------
3069 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3070 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3071 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3072 Monthly Weather Review, Vol. 130, 1812-1828
3073 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3074 """
3075 Nt, Ny, Nx = y_3d.shape
3077 # Max coefficient
3078 Nmin = min(Nx - 1, Ny - 1)
3080 # Create alpha matrix (of wavenumbers)
3081 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3083 # Prepare output array
3084 ps_array = np.zeros((Nt, Nmin))
3086 # Loop over time to get spectrum for each time.
3087 for t in range(Nt):
3088 y_2d = y_3d[t]
3090 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3091 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3092 # cosine basis functions at different spatial frequencies.
3094 # normalise spectrum to allow comparison between models.
3095 fkk = fft.dctn(y_2d, norm="ortho")
3097 # Normalise fkk
3098 fkk = fkk / np.sqrt(Ny * Nx)
3100 # calculate variance of spectral coefficient
3101 sigma_2 = fkk**2 / Nx / Ny
3103 # Group ellipses of alphas into the same wavenumber k/Nmin
3104 for k in range(1, Nmin + 1):
3105 alpha = k / Nmin
3106 alpha_p1 = (k + 1) / Nmin
3108 # Sum up elements matching k
3109 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3110 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3112 return ps_array
3115def _create_alpha_matrix(Ny, Nx):
3116 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3118 Parameters
3119 ----------
3120 Ny, Nx:
3121 Dimensions of the 2D field for which the power spectra is calculated. Used to
3122 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3123 single-scale parameter.
3125 Returns: alpha_matrix
3126 normalisation of 2D wavenumber axes, transforming the spectral domain into
3127 an elliptic coordinate system.
3129 """
3130 # Create x_indices: each row is [1, 2, ..., Nx]
3131 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3133 # Create y_indices: each column is [1, 2, ..., Ny]
3134 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3136 # Compute alpha_matrix
3137 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3139 return alpha_matrix