Coverage for src / CSET / operators / plot.py: 86%
953 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-21 14:08 +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 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 410 ↛ 412line 410 didn't jump to line 412 because the condition on line 410 was never true
411 # Define Transverse Mercator projection for TM inputs.
412 projection = ccrs.TransverseMercator(
413 central_longitude=coord_system.longitude_of_central_meridian,
414 central_latitude=coord_system.latitude_of_projection_origin,
415 false_easting=coord_system.false_easting,
416 false_northing=coord_system.false_northing,
417 scale_factor=coord_system.scale_factor_at_central_meridian,
418 )
419 crs = projection
420 else:
421 # Define regular map projection for non-rotated pole inputs.
422 # Alternatives might include e.g. for global model outputs:
423 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
424 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
425 projection = ccrs.PlateCarree(central_longitude=central_longitude)
426 crs = ccrs.PlateCarree()
428 # Define axes for plot (or subplot) with required map projection.
429 if subplot is not None:
430 axes = figure.add_subplot(
431 grid_size, grid_size, subplot, projection=projection
432 )
433 else:
434 axes = figure.add_subplot(projection=projection)
436 # Add coastlines if cube contains x and y map coordinates.
437 if cmap.name in ["viridis", "Greys"]:
438 coastcol = "magenta"
439 else:
440 coastcol = "black"
441 logging.debug("Plotting coastlines in colour %s.", coastcol)
442 axes.coastlines(resolution="10m", color=coastcol)
444 # If is lat/lon spatial map, fix extent to keep plot tight.
445 # Specifying crs within set_extent helps ensure only data region is shown.
446 if isinstance(coord_system, iris.coord_systems.GeogCS):
447 axes.set_extent([x1, x2, y1, y2], crs=crs)
449 except ValueError:
450 # Skip if not both x and y map coordinates.
451 axes = figure.gca()
452 pass
454 return axes
457def _get_plot_resolution() -> int:
458 """Get resolution of rasterised plots in pixels per inch."""
459 return get_recipe_metadata().get("plot_resolution", 100)
462def _plot_and_save_spatial_plot(
463 cube: iris.cube.Cube,
464 filename: str,
465 title: str,
466 method: Literal["contourf", "pcolormesh"],
467 **kwargs,
468):
469 """Plot and save a spatial plot.
471 Parameters
472 ----------
473 cube: Cube
474 2 dimensional (lat and lon) Cube of the data to plot.
475 filename: str
476 Filename of the plot to write.
477 title: str
478 Plot title.
479 method: "contourf" | "pcolormesh"
480 The plotting method to use.
481 """
482 # Setup plot details, size, resolution, etc.
483 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
485 # Specify the color bar
486 cmap, levels, norm = _colorbar_map_levels(cube)
488 # Setup plot map projection, extent and coastlines.
489 axes = _setup_spatial_map(cube, fig, cmap)
491 # Plot the field.
492 if method == "contourf":
493 # Filled contour plot of the field.
494 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
495 elif method == "pcolormesh":
496 try:
497 vmin = min(levels)
498 vmax = max(levels)
499 except TypeError:
500 vmin, vmax = None, None
501 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
502 # if levels are defined.
503 if norm is not None:
504 vmin = None
505 vmax = None
506 logging.debug("Plotting using defined levels.")
507 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
508 else:
509 raise ValueError(f"Unknown plotting method: {method}")
511 # Check to see if transect, and if so, adjust y axis.
512 if is_transect(cube):
513 if "pressure" in [coord.name() for coord in cube.coords()]:
514 axes.invert_yaxis()
515 axes.set_yscale("log")
516 axes.set_ylim(1100, 100)
517 # If both model_level_number and level_height exists, iplt can construct
518 # plot as a function of height above orography (NOT sea level).
519 elif {"model_level_number", "level_height"}.issubset( 519 ↛ 524line 519 didn't jump to line 524 because the condition on line 519 was always true
520 {coord.name() for coord in cube.coords()}
521 ):
522 axes.set_yscale("log")
524 axes.set_title(
525 f"{title}\n"
526 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
527 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
528 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
529 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
530 fontsize=16,
531 )
533 else:
534 # Add title.
535 axes.set_title(title, fontsize=16)
537 # Add watermark with min/max/mean. Currently not user togglable.
538 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
539 axes.annotate(
540 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
541 xy=(1, -0.05),
542 xycoords="axes fraction",
543 xytext=(-5, 5),
544 textcoords="offset points",
545 ha="right",
546 va="bottom",
547 size=11,
548 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
549 )
551 # Add colour bar.
552 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
553 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
554 # add ticks and tick_labels for every levels if less than 20 levels exist
555 if levels is not None and len(levels) < 20:
556 cbar.set_ticks(levels)
557 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
558 if "visibility" in cube.name(): 558 ↛ 559line 558 didn't jump to line 559 because the condition on line 558 was never true
559 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
560 logging.debug("Set colorbar ticks and labels.")
562 # Save plot.
563 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
564 logging.info("Saved spatial plot to %s", filename)
565 plt.close(fig)
568def _plot_and_save_postage_stamp_spatial_plot(
569 cube: iris.cube.Cube,
570 filename: str,
571 stamp_coordinate: str,
572 title: str,
573 method: Literal["contourf", "pcolormesh"],
574 **kwargs,
575):
576 """Plot postage stamp spatial plots from an ensemble.
578 Parameters
579 ----------
580 cube: Cube
581 Iris cube of data to be plotted. It must have the stamp coordinate.
582 filename: str
583 Filename of the plot to write.
584 stamp_coordinate: str
585 Coordinate that becomes different plots.
586 method: "contourf" | "pcolormesh"
587 The plotting method to use.
589 Raises
590 ------
591 ValueError
592 If the cube doesn't have the right dimensions.
593 """
594 # Use the smallest square grid that will fit the members.
595 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
597 fig = plt.figure(figsize=(10, 10))
599 # Specify the color bar
600 cmap, levels, norm = _colorbar_map_levels(cube)
602 # Make a subplot for each member.
603 for member, subplot in zip(
604 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
605 ):
606 # Setup subplot map projection, extent and coastlines.
607 axes = _setup_spatial_map(
608 member, fig, cmap, grid_size=grid_size, subplot=subplot
609 )
610 if method == "contourf":
611 # Filled contour plot of the field.
612 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
613 elif method == "pcolormesh":
614 if levels is not None:
615 vmin = min(levels)
616 vmax = max(levels)
617 else:
618 raise TypeError("Unknown vmin and vmax range.")
619 vmin, vmax = None, None
620 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
621 # if levels are defined.
622 if norm is not None: 622 ↛ 623line 622 didn't jump to line 623 because the condition on line 622 was never true
623 vmin = None
624 vmax = None
625 # pcolormesh plot of the field.
626 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
627 else:
628 raise ValueError(f"Unknown plotting method: {method}")
629 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
630 axes.set_axis_off()
632 # Put the shared colorbar in its own axes.
633 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])
634 colorbar = fig.colorbar(
635 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
636 )
637 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
639 # Overall figure title.
640 fig.suptitle(title, fontsize=16)
642 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
643 logging.info("Saved contour postage stamp plot to %s", filename)
644 plt.close(fig)
647def _plot_and_save_line_series(
648 cubes: iris.cube.CubeList,
649 coords: list[iris.coords.Coord],
650 ensemble_coord: str,
651 filename: str,
652 title: str,
653 **kwargs,
654):
655 """Plot and save a 1D line series.
657 Parameters
658 ----------
659 cubes: Cube or CubeList
660 Cube or CubeList containing the cubes to plot on the y-axis.
661 coords: list[Coord]
662 Coordinates to plot on the x-axis, one per cube.
663 ensemble_coord: str
664 Ensemble coordinate in the cube.
665 filename: str
666 Filename of the plot to write.
667 title: str
668 Plot title.
669 """
670 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
672 model_colors_map = _get_model_colors_map(cubes)
674 # Store min/max ranges.
675 y_levels = []
677 # Check match-up across sequence coords gives consistent sizes
678 _validate_cubes_coords(cubes, coords)
680 for cube, coord in zip(cubes, coords, strict=True):
681 label = None
682 color = "black"
683 if model_colors_map:
684 label = cube.attributes.get("model_name")
685 color = model_colors_map.get(label)
686 for cube_slice in cube.slices_over(ensemble_coord):
687 # Label with (control) if part of an ensemble or not otherwise.
688 if cube_slice.coord(ensemble_coord).points == [0]:
689 iplt.plot(
690 coord,
691 cube_slice,
692 color=color,
693 marker="o",
694 ls="-",
695 lw=3,
696 label=f"{label} (control)"
697 if len(cube.coord(ensemble_coord).points) > 1
698 else label,
699 )
700 # Label with (perturbed) if part of an ensemble and not the control.
701 else:
702 iplt.plot(
703 coord,
704 cube_slice,
705 color=color,
706 ls="-",
707 lw=1.5,
708 alpha=0.75,
709 label=f"{label} (member)",
710 )
712 # Calculate the global min/max if multiple cubes are given.
713 _, levels, _ = _colorbar_map_levels(cube, axis="y")
714 if levels is not None: 714 ↛ 715line 714 didn't jump to line 715 because the condition on line 714 was never true
715 y_levels.append(min(levels))
716 y_levels.append(max(levels))
718 # Get the current axes.
719 ax = plt.gca()
721 # Add some labels and tweak the style.
722 # check if cubes[0] works for single cube if not CubeList
723 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
724 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
725 ax.set_title(title, fontsize=16)
727 ax.ticklabel_format(axis="y", useOffset=False)
728 ax.tick_params(axis="x", labelrotation=15)
729 ax.tick_params(axis="both", labelsize=12)
731 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
732 if y_levels: 732 ↛ 733line 732 didn't jump to line 733 because the condition on line 732 was never true
733 ax.set_ylim(min(y_levels), max(y_levels))
734 # Add zero line.
735 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
736 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
737 logging.debug(
738 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
739 )
740 else:
741 ax.autoscale()
743 # Add gridlines
744 ax.grid(linestyle="--", color="grey", linewidth=1)
745 # Ientify unique labels for legend
746 handles = list(
747 {
748 label: handle
749 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
750 }.values()
751 )
752 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
754 # Save plot.
755 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
756 logging.info("Saved line plot to %s", filename)
757 plt.close(fig)
760def _plot_and_save_vertical_line_series(
761 cubes: iris.cube.CubeList,
762 coords: list[iris.coords.Coord],
763 ensemble_coord: str,
764 filename: str,
765 series_coordinate: str,
766 title: str,
767 vmin: float,
768 vmax: float,
769 **kwargs,
770):
771 """Plot and save a 1D line series in vertical.
773 Parameters
774 ----------
775 cubes: CubeList
776 1 dimensional Cube or CubeList of the data to plot on x-axis.
777 coord: list[Coord]
778 Coordinates to plot on the y-axis, one per cube.
779 ensemble_coord: str
780 Ensemble coordinate in the cube.
781 filename: str
782 Filename of the plot to write.
783 series_coordinate: str
784 Coordinate to use as vertical axis.
785 title: str
786 Plot title.
787 vmin: float
788 Minimum value for the x-axis.
789 vmax: float
790 Maximum value for the x-axis.
791 """
792 # plot the vertical pressure axis using log scale
793 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
795 model_colors_map = _get_model_colors_map(cubes)
797 # Check match-up across sequence coords gives consistent sizes
798 _validate_cubes_coords(cubes, coords)
800 for cube, coord in zip(cubes, coords, strict=True):
801 label = None
802 color = "black"
803 if model_colors_map: 803 ↛ 804line 803 didn't jump to line 804 because the condition on line 803 was never true
804 label = cube.attributes.get("model_name")
805 color = model_colors_map.get(label)
807 for cube_slice in cube.slices_over(ensemble_coord):
808 # If ensemble data given plot control member with (control)
809 # unless single forecast.
810 if cube_slice.coord(ensemble_coord).points == [0]:
811 iplt.plot(
812 cube_slice,
813 coord,
814 color=color,
815 marker="o",
816 ls="-",
817 lw=3,
818 label=f"{label} (control)"
819 if len(cube.coord(ensemble_coord).points) > 1
820 else label,
821 )
822 # If ensemble data given plot perturbed members with (perturbed).
823 else:
824 iplt.plot(
825 cube_slice,
826 coord,
827 color=color,
828 ls="-",
829 lw=1.5,
830 alpha=0.75,
831 label=f"{label} (member)",
832 )
834 # Get the current axis
835 ax = plt.gca()
837 # Special handling for pressure level data.
838 if series_coordinate == "pressure": 838 ↛ 860line 838 didn't jump to line 860 because the condition on line 838 was always true
839 # Invert y-axis and set to log scale.
840 ax.invert_yaxis()
841 ax.set_yscale("log")
843 # Define y-ticks and labels for pressure log axis.
844 y_tick_labels = [
845 "1000",
846 "850",
847 "700",
848 "500",
849 "300",
850 "200",
851 "100",
852 ]
853 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
855 # Set y-axis limits and ticks.
856 ax.set_ylim(1100, 100)
858 # Test if series_coordinate is model level data. The UM data uses
859 # model_level_number and lfric uses full_levels as coordinate.
860 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
861 # Define y-ticks and labels for vertical axis.
862 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
863 y_tick_labels = [str(int(i)) for i in y_ticks]
864 ax.set_ylim(min(y_ticks), max(y_ticks))
866 ax.set_yticks(y_ticks)
867 ax.set_yticklabels(y_tick_labels)
869 # Set x-axis limits.
870 ax.set_xlim(vmin, vmax)
871 # Mark y=0 if present in plot.
872 if vmin < 0.0 and vmax > 0.0: 872 ↛ 873line 872 didn't jump to line 873 because the condition on line 872 was never true
873 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
875 # Add some labels and tweak the style.
876 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
877 ax.set_xlabel(
878 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
879 )
880 ax.set_title(title, fontsize=16)
881 ax.ticklabel_format(axis="x")
882 ax.tick_params(axis="y")
883 ax.tick_params(axis="both", labelsize=12)
885 # Add gridlines
886 ax.grid(linestyle="--", color="grey", linewidth=1)
887 # Ientify unique labels for legend
888 handles = list(
889 {
890 label: handle
891 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
892 }.values()
893 )
894 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
896 # Save plot.
897 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
898 logging.info("Saved line plot to %s", filename)
899 plt.close(fig)
902def _plot_and_save_scatter_plot(
903 cube_x: iris.cube.Cube | iris.cube.CubeList,
904 cube_y: iris.cube.Cube | iris.cube.CubeList,
905 filename: str,
906 title: str,
907 one_to_one: bool,
908 model_names: list[str] = None,
909 **kwargs,
910):
911 """Plot and save a 2D scatter plot.
913 Parameters
914 ----------
915 cube_x: Cube | CubeList
916 1 dimensional Cube or CubeList of the data to plot on x-axis.
917 cube_y: Cube | CubeList
918 1 dimensional Cube or CubeList of the data to plot on y-axis.
919 filename: str
920 Filename of the plot to write.
921 title: str
922 Plot title.
923 one_to_one: bool
924 Whether a 1:1 line is plotted.
925 """
926 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
927 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
928 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
929 # over the pairs simultaneously.
931 # Ensure cube_x and cube_y are iterable
932 cube_x_iterable = iter_maybe(cube_x)
933 cube_y_iterable = iter_maybe(cube_y)
935 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
936 iplt.scatter(cube_x_iter, cube_y_iter)
937 if one_to_one is True:
938 plt.plot(
939 [
940 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
941 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
942 ],
943 [
944 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
945 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
946 ],
947 "k",
948 linestyle="--",
949 )
950 ax = plt.gca()
952 # Add some labels and tweak the style.
953 if model_names is None:
954 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
955 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
956 else:
957 # Add the model names, these should be order of base (x) and other (y).
958 ax.set_xlabel(
959 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
960 )
961 ax.set_ylabel(
962 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
963 )
964 ax.set_title(title, fontsize=16)
965 ax.ticklabel_format(axis="y", useOffset=False)
966 ax.tick_params(axis="x", labelrotation=15)
967 ax.tick_params(axis="both", labelsize=12)
968 ax.autoscale()
970 # Save plot.
971 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
972 logging.info("Saved scatter plot to %s", filename)
973 plt.close(fig)
976def _plot_and_save_vector_plot(
977 cube_u: iris.cube.Cube,
978 cube_v: iris.cube.Cube,
979 filename: str,
980 title: str,
981 method: Literal["contourf", "pcolormesh"],
982 **kwargs,
983):
984 """Plot and save a 2D vector plot.
986 Parameters
987 ----------
988 cube_u: Cube
989 2 dimensional Cube of u component of the data.
990 cube_v: Cube
991 2 dimensional Cube of v component of the data.
992 filename: str
993 Filename of the plot to write.
994 title: str
995 Plot title.
996 """
997 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
999 # Create a cube containing the magnitude of the vector field.
1000 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1001 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1003 # Specify the color bar
1004 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1006 # Setup plot map projection, extent and coastlines.
1007 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1009 if method == "contourf": 1009 ↛ 1012line 1009 didn't jump to line 1012 because the condition on line 1009 was always true
1010 # Filled contour plot of the field.
1011 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1012 elif method == "pcolormesh":
1013 try:
1014 vmin = min(levels)
1015 vmax = max(levels)
1016 except TypeError:
1017 vmin, vmax = None, None
1018 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1019 # if levels are defined.
1020 if norm is not None:
1021 vmin = None
1022 vmax = None
1023 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1024 else:
1025 raise ValueError(f"Unknown plotting method: {method}")
1027 # Check to see if transect, and if so, adjust y axis.
1028 if is_transect(cube_vec_mag): 1028 ↛ 1029line 1028 didn't jump to line 1029 because the condition on line 1028 was never true
1029 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1030 axes.invert_yaxis()
1031 axes.set_yscale("log")
1032 axes.set_ylim(1100, 100)
1033 # If both model_level_number and level_height exists, iplt can construct
1034 # plot as a function of height above orography (NOT sea level).
1035 elif {"model_level_number", "level_height"}.issubset(
1036 {coord.name() for coord in cube_vec_mag.coords()}
1037 ):
1038 axes.set_yscale("log")
1040 axes.set_title(
1041 f"{title}\n"
1042 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1043 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1044 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1045 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1046 fontsize=16,
1047 )
1049 else:
1050 # Add title.
1051 axes.set_title(title, fontsize=16)
1053 # Add watermark with min/max/mean. Currently not user togglable.
1054 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1055 axes.annotate(
1056 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}",
1057 xy=(1, -0.05),
1058 xycoords="axes fraction",
1059 xytext=(-5, 5),
1060 textcoords="offset points",
1061 ha="right",
1062 va="bottom",
1063 size=11,
1064 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1065 )
1067 # Add colour bar.
1068 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1069 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1070 # add ticks and tick_labels for every levels if less than 20 levels exist
1071 if levels is not None and len(levels) < 20: 1071 ↛ 1072line 1071 didn't jump to line 1072 because the condition on line 1071 was never true
1072 cbar.set_ticks(levels)
1073 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1075 # 30 barbs along the longest axis of the plot, or a barb per point for data
1076 # with less than 30 points.
1077 step = max(max(cube_u.shape) // 30, 1)
1078 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1080 # Save plot.
1081 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1082 logging.info("Saved vector plot to %s", filename)
1083 plt.close(fig)
1086def _plot_and_save_histogram_series(
1087 cubes: iris.cube.Cube | iris.cube.CubeList,
1088 filename: str,
1089 title: str,
1090 vmin: float,
1091 vmax: float,
1092 **kwargs,
1093):
1094 """Plot and save a histogram series.
1096 Parameters
1097 ----------
1098 cubes: Cube or CubeList
1099 2 dimensional Cube or CubeList of the data to plot as histogram.
1100 filename: str
1101 Filename of the plot to write.
1102 title: str
1103 Plot title.
1104 vmin: float
1105 minimum for colorbar
1106 vmax: float
1107 maximum for colorbar
1108 """
1109 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1110 ax = plt.gca()
1112 model_colors_map = _get_model_colors_map(cubes)
1114 # Set default that histograms will produce probability density function
1115 # at each bin (integral over range sums to 1).
1116 density = True
1118 for cube in iter_maybe(cubes):
1119 # Easier to check title (where var name originates)
1120 # than seeing if long names exist etc.
1121 # Exception case, where distribution better fits log scales/bins.
1122 if "surface_microphysical" in title:
1123 if "amount" in title: 1123 ↛ 1125line 1123 didn't jump to line 1125 because the condition on line 1123 was never true
1124 # Compute histogram following Klingaman et al. (2017): ASoP
1125 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1126 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1127 density = False
1128 else:
1129 bins = 10.0 ** (
1130 np.arange(-10, 27, 1) / 10.0
1131 ) # Suggestion from RMED toolbox.
1132 bins = np.insert(bins, 0, 0)
1133 ax.set_yscale("log")
1134 vmin = bins[1]
1135 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1136 ax.set_xscale("log")
1137 elif "lightning" in title:
1138 bins = [0, 1, 2, 3, 4, 5]
1139 else:
1140 bins = np.linspace(vmin, vmax, 51)
1141 logging.debug(
1142 "Plotting histogram with %s bins %s - %s.",
1143 np.size(bins),
1144 np.min(bins),
1145 np.max(bins),
1146 )
1148 # Reshape cube data into a single array to allow for a single histogram.
1149 # Otherwise we plot xdim histograms stacked.
1150 cube_data_1d = (cube.data).flatten()
1152 label = None
1153 color = "black"
1154 if model_colors_map: 1154 ↛ 1155line 1154 didn't jump to line 1155 because the condition on line 1154 was never true
1155 label = cube.attributes.get("model_name")
1156 color = model_colors_map[label]
1157 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1159 # Compute area under curve.
1160 if "surface_microphysical" in title and "amount" in title: 1160 ↛ 1161line 1160 didn't jump to line 1161 because the condition on line 1160 was never true
1161 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1162 x = x * bin_mean / x.sum()
1163 x = x[1:]
1164 y = y[1:]
1166 ax.plot(
1167 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1168 )
1170 # Add some labels and tweak the style.
1171 ax.set_title(title, fontsize=16)
1172 ax.set_xlabel(
1173 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1174 )
1175 ax.set_ylabel("Normalised probability density", fontsize=14)
1176 if "surface_microphysical" in title and "amount" in title: 1176 ↛ 1177line 1176 didn't jump to line 1177 because the condition on line 1176 was never true
1177 ax.set_ylabel(
1178 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1179 )
1180 ax.set_xlim(vmin, vmax)
1181 ax.tick_params(axis="both", labelsize=12)
1183 # Overlay grid-lines onto histogram plot.
1184 ax.grid(linestyle="--", color="grey", linewidth=1)
1185 if model_colors_map: 1185 ↛ 1186line 1185 didn't jump to line 1186 because the condition on line 1185 was never true
1186 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1188 # Save plot.
1189 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1190 logging.info("Saved line plot to %s", filename)
1191 plt.close(fig)
1194def _plot_and_save_postage_stamp_histogram_series(
1195 cube: iris.cube.Cube,
1196 filename: str,
1197 title: str,
1198 stamp_coordinate: str,
1199 vmin: float,
1200 vmax: float,
1201 **kwargs,
1202):
1203 """Plot and save postage (ensemble members) stamps for a histogram series.
1205 Parameters
1206 ----------
1207 cube: Cube
1208 2 dimensional Cube of the data to plot as histogram.
1209 filename: str
1210 Filename of the plot to write.
1211 title: str
1212 Plot title.
1213 stamp_coordinate: str
1214 Coordinate that becomes different plots.
1215 vmin: float
1216 minimum for pdf x-axis
1217 vmax: float
1218 maximum for pdf x-axis
1219 """
1220 # Use the smallest square grid that will fit the members.
1221 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1223 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1224 # Make a subplot for each member.
1225 for member, subplot in zip(
1226 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1227 ):
1228 # Implicit interface is much easier here, due to needing to have the
1229 # cartopy GeoAxes generated.
1230 plt.subplot(grid_size, grid_size, subplot)
1231 # Reshape cube data into a single array to allow for a single histogram.
1232 # Otherwise we plot xdim histograms stacked.
1233 member_data_1d = (member.data).flatten()
1234 plt.hist(member_data_1d, density=True, stacked=True)
1235 ax = plt.gca()
1236 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1237 ax.set_xlim(vmin, vmax)
1239 # Overall figure title.
1240 fig.suptitle(title, fontsize=16)
1242 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1243 logging.info("Saved histogram postage stamp plot to %s", filename)
1244 plt.close(fig)
1247def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1248 cube: iris.cube.Cube,
1249 filename: str,
1250 title: str,
1251 stamp_coordinate: str,
1252 vmin: float,
1253 vmax: float,
1254 **kwargs,
1255):
1256 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1257 ax.set_title(title, fontsize=16)
1258 ax.set_xlim(vmin, vmax)
1259 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1260 ax.set_ylabel("normalised probability density", fontsize=14)
1261 # Loop over all slices along the stamp_coordinate
1262 for member in cube.slices_over(stamp_coordinate):
1263 # Flatten the member data to 1D
1264 member_data_1d = member.data.flatten()
1265 # Plot the histogram using plt.hist
1266 plt.hist(
1267 member_data_1d,
1268 density=True,
1269 stacked=True,
1270 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1271 )
1273 # Add a legend
1274 ax.legend(fontsize=16)
1276 # Save the figure to a file
1277 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1279 # Close the figure
1280 plt.close(fig)
1283def _plot_and_save_scattermap_plot(
1284 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1285):
1286 """Plot and save a geographical scatter plot.
1288 Parameters
1289 ----------
1290 cube: Cube
1291 1 dimensional Cube of the data points with auxiliary latitude and
1292 longitude coordinates,
1293 filename: str
1294 Filename of the plot to write.
1295 title: str
1296 Plot title.
1297 projection: str
1298 Mapping projection to be used by cartopy.
1299 """
1300 # Setup plot details, size, resolution, etc.
1301 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1302 if projection is not None:
1303 # Apart from the default, the only projection we currently support is
1304 # a stereographic projection over the North Pole.
1305 if projection == "NP_Stereo":
1306 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1307 else:
1308 raise ValueError(f"Unknown projection: {projection}")
1309 else:
1310 axes = plt.axes(projection=ccrs.PlateCarree())
1312 # Scatter plot of the field. The marker size is chosen to give
1313 # symbols that decrease in size as the number of observations
1314 # increases, although the fraction of the figure covered by
1315 # symbols increases roughly as N^(1/2), disregarding overlaps,
1316 # and has been selected for the default figure size of (10, 10).
1317 # Should this be changed, the marker size should be adjusted in
1318 # proportion to the area of the figure.
1319 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1320 klon = None
1321 klat = None
1322 for kc in range(len(cube.aux_coords)):
1323 if cube.aux_coords[kc].standard_name == "latitude":
1324 klat = kc
1325 elif cube.aux_coords[kc].standard_name == "longitude":
1326 klon = kc
1327 scatter_map = iplt.scatter(
1328 cube.aux_coords[klon],
1329 cube.aux_coords[klat],
1330 c=cube.data[:],
1331 s=mrk_size,
1332 cmap="jet",
1333 edgecolors="k",
1334 )
1336 # Add coastlines.
1337 try:
1338 axes.coastlines(resolution="10m")
1339 except AttributeError:
1340 pass
1342 # Add title.
1343 axes.set_title(title, fontsize=16)
1345 # Add colour bar.
1346 cbar = fig.colorbar(scatter_map)
1347 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1349 # Save plot.
1350 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1351 logging.info("Saved geographical scatter plot to %s", filename)
1352 plt.close(fig)
1355def _plot_and_save_power_spectrum_series(
1356 cubes: iris.cube.Cube | iris.cube.CubeList,
1357 filename: str,
1358 title: str,
1359 **kwargs,
1360):
1361 """Plot and save a power spectrum series.
1363 Parameters
1364 ----------
1365 cubes: Cube or CubeList
1366 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1367 filename: str
1368 Filename of the plot to write.
1369 title: str
1370 Plot title.
1371 """
1372 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1373 ax = plt.gca()
1375 model_colors_map = _get_model_colors_map(cubes)
1377 for cube in iter_maybe(cubes):
1378 # Calculate power spectrum
1380 # Extract time coordinate and convert to datetime
1381 time_coord = cube.coord("time")
1382 time_points = time_coord.units.num2date(time_coord.points)
1384 # Choose one time point (e.g., the first one)
1385 target_time = time_points[0]
1387 # Bind target_time inside the lambda using a default argument
1388 time_constraint = iris.Constraint(
1389 time=lambda cell, target_time=target_time: cell.point == target_time
1390 )
1392 cube = cube.extract(time_constraint)
1394 if cube.ndim == 2:
1395 cube_3d = cube.data[np.newaxis, :, :]
1396 logging.debug("Adding in new axis for a 2 dimensional cube.")
1397 elif cube.ndim == 3: 1397 ↛ 1398line 1397 didn't jump to line 1398 because the condition on line 1397 was never true
1398 cube_3d = cube.data
1399 else:
1400 raise ValueError("Cube dimensions unsuitable for power spectra code")
1401 raise ValueError(
1402 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1403 )
1405 # Calculate spectra
1406 ps_array = _DCT_ps(cube_3d)
1408 ps_cube = iris.cube.Cube(
1409 ps_array,
1410 long_name="power_spectra",
1411 )
1413 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1415 # Create a frequency/wavelength array for coordinate
1416 ps_len = ps_cube.data.shape[1]
1417 freqs = np.arange(1, ps_len + 1)
1418 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1420 # Convert datetime to numeric time using original units
1421 numeric_time = time_coord.units.date2num(time_points)
1422 # Create a new DimCoord with numeric time
1423 new_time_coord = iris.coords.DimCoord(
1424 numeric_time, standard_name="time", units=time_coord.units
1425 )
1427 # Add time and frequency coordinate to spectra cube.
1428 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1429 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1431 # Extract data from the cube
1432 frequency = ps_cube.coord("frequency").points
1433 power_spectrum = ps_cube.data
1435 label = None
1436 color = "black"
1437 if model_colors_map: 1437 ↛ 1438line 1437 didn't jump to line 1438 because the condition on line 1437 was never true
1438 label = ps_cube.attributes.get("model_name")
1439 color = model_colors_map[label]
1440 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1442 # Add some labels and tweak the style.
1443 ax.set_title(title, fontsize=16)
1444 ax.set_xlabel("Wavenumber", fontsize=14)
1445 ax.set_ylabel("Power", fontsize=14)
1446 ax.tick_params(axis="both", labelsize=12)
1448 # Set log-log scale
1449 ax.set_xscale("log")
1450 ax.set_yscale("log")
1452 # Overlay grid-lines onto power spectrum plot.
1453 ax.grid(linestyle="--", color="grey", linewidth=1)
1454 if model_colors_map: 1454 ↛ 1455line 1454 didn't jump to line 1455 because the condition on line 1454 was never true
1455 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1457 # Save plot.
1458 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1459 logging.info("Saved line plot to %s", filename)
1460 plt.close(fig)
1463def _plot_and_save_postage_stamp_power_spectrum_series(
1464 cube: iris.cube.Cube,
1465 filename: str,
1466 title: str,
1467 stamp_coordinate: str,
1468 **kwargs,
1469):
1470 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1472 Parameters
1473 ----------
1474 cube: Cube
1475 2 dimensional Cube of the data to plot as power spectrum.
1476 filename: str
1477 Filename of the plot to write.
1478 title: str
1479 Plot title.
1480 stamp_coordinate: str
1481 Coordinate that becomes different plots.
1482 """
1483 # Use the smallest square grid that will fit the members.
1484 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1486 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1487 # Make a subplot for each member.
1488 for member, subplot in zip(
1489 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1490 ):
1491 # Implicit interface is much easier here, due to needing to have the
1492 # cartopy GeoAxes generated.
1493 plt.subplot(grid_size, grid_size, subplot)
1495 frequency = member.coord("frequency").points
1496 power_spectrum = member.data
1498 ax = plt.gca()
1499 ax.plot(frequency, power_spectrum[0])
1500 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1502 # Overall figure title.
1503 fig.suptitle(title, fontsize=16)
1505 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1506 logging.info("Saved power spectra postage stamp plot to %s", filename)
1507 plt.close(fig)
1510def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1511 cube: iris.cube.Cube,
1512 filename: str,
1513 title: str,
1514 stamp_coordinate: str,
1515 **kwargs,
1516):
1517 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1518 ax.set_title(title, fontsize=16)
1519 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1520 ax.set_ylabel("Power", fontsize=14)
1521 # Loop over all slices along the stamp_coordinate
1522 for member in cube.slices_over(stamp_coordinate):
1523 frequency = member.coord("frequency").points
1524 power_spectrum = member.data
1525 ax.plot(
1526 frequency,
1527 power_spectrum[0],
1528 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1529 )
1531 # Add a legend
1532 ax.legend(fontsize=16)
1534 # Save the figure to a file
1535 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1537 # Close the figure
1538 plt.close(fig)
1541def _spatial_plot(
1542 method: Literal["contourf", "pcolormesh"],
1543 cube: iris.cube.Cube,
1544 filename: str | None,
1545 sequence_coordinate: str,
1546 stamp_coordinate: str,
1547 **kwargs,
1548):
1549 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1551 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1552 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1553 is present then postage stamp plots will be produced.
1555 Parameters
1556 ----------
1557 method: "contourf" | "pcolormesh"
1558 The plotting method to use.
1559 cube: Cube
1560 Iris cube of the data to plot. It should have two spatial dimensions,
1561 such as lat and lon, and may also have a another two dimension to be
1562 plotted sequentially and/or as postage stamp plots.
1563 filename: str | None
1564 Name of the plot to write, used as a prefix for plot sequences. If None
1565 uses the recipe name.
1566 sequence_coordinate: str
1567 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1568 This coordinate must exist in the cube.
1569 stamp_coordinate: str
1570 Coordinate about which to plot postage stamp plots. Defaults to
1571 ``"realization"``.
1573 Raises
1574 ------
1575 ValueError
1576 If the cube doesn't have the right dimensions.
1577 TypeError
1578 If the cube isn't a single cube.
1579 """
1580 recipe_title = get_recipe_metadata().get("title", "Untitled")
1582 # Ensure we have a name for the plot file.
1583 if filename is None:
1584 filename = slugify(recipe_title)
1586 # Ensure we've got a single cube.
1587 cube = _check_single_cube(cube)
1589 # Make postage stamp plots if stamp_coordinate exists and has more than a
1590 # single point.
1591 plotting_func = _plot_and_save_spatial_plot
1592 try:
1593 if cube.coord(stamp_coordinate).shape[0] > 1:
1594 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1595 except iris.exceptions.CoordinateNotFoundError:
1596 pass
1598 # Produce a geographical scatter plot if the data have a
1599 # dimension called observation or model_obs_error
1600 if any( 1600 ↛ 1604line 1600 didn't jump to line 1604 because the condition on line 1600 was never true
1601 crd.var_name == "station" or crd.var_name == "model_obs_error"
1602 for crd in cube.coords()
1603 ):
1604 plotting_func = _plot_and_save_scattermap_plot
1606 # Must have a sequence coordinate.
1607 try:
1608 cube.coord(sequence_coordinate)
1609 except iris.exceptions.CoordinateNotFoundError as err:
1610 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1612 # Create a plot for each value of the sequence coordinate.
1613 plot_index = []
1614 nplot = np.size(cube.coord(sequence_coordinate).points)
1615 for cube_slice in cube.slices_over(sequence_coordinate):
1616 # Use sequence value so multiple sequences can merge.
1617 sequence_value = cube_slice.coord(sequence_coordinate).points[0]
1618 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
1619 coord = cube_slice.coord(sequence_coordinate)
1620 # Format the coordinate value in a unit appropriate way.
1621 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
1622 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
1623 if nplot == 1 and coord.has_bounds:
1624 if np.size(coord.bounds) > 1:
1625 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
1626 # Do the actual plotting.
1627 plotting_func(
1628 cube_slice,
1629 filename=plot_filename,
1630 stamp_coordinate=stamp_coordinate,
1631 title=title,
1632 method=method,
1633 **kwargs,
1634 )
1635 plot_index.append(plot_filename)
1637 # Add list of plots to plot metadata.
1638 complete_plot_index = _append_to_plot_index(plot_index)
1640 # Make a page to display the plots.
1641 _make_plot_html_page(complete_plot_index)
1644def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1645 """Get colourmap for mask.
1647 If "mask_for_" appears anywhere in the name of a cube this function will be called
1648 regardless of the name of the variable to ensure a consistent plot.
1650 Parameters
1651 ----------
1652 cube: Cube
1653 Cube of variable for which the colorbar information is desired.
1654 axis: "x", "y", optional
1655 Select the levels for just this axis of a line plot. The min and max
1656 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1657 setting a universal range is not desirable (e.g. temperature), users
1658 can set ymin/ymax values to "auto" in the colorbar definitions file.
1659 Where no additional xmin/xmax or ymin/ymax values are provided, the
1660 axis bounds default to use the vmin/vmax values provided.
1662 Returns
1663 -------
1664 cmap:
1665 Matplotlib colormap.
1666 levels:
1667 List of levels to use for plotting. For continuous plots the min and max
1668 should be taken as the range.
1669 norm:
1670 BoundaryNorm information.
1671 """
1672 if "difference" not in cube.long_name:
1673 if axis:
1674 levels = [0, 1]
1675 # Complete settings based on levels.
1676 return None, levels, None
1677 else:
1678 # Define the levels and colors.
1679 levels = [0, 1, 2]
1680 colors = ["white", "dodgerblue"]
1681 # Create a custom color map.
1682 cmap = mcolors.ListedColormap(colors)
1683 # Normalize the levels.
1684 norm = mcolors.BoundaryNorm(levels, cmap.N)
1685 logging.debug("Colourmap for %s.", cube.long_name)
1686 return cmap, levels, norm
1687 else:
1688 if axis:
1689 levels = [-1, 1]
1690 return None, levels, None
1691 else:
1692 # Search for if mask difference, set to +/- 0.5 as values plotted <
1693 # not <=.
1694 levels = [-2, -0.5, 0.5, 2]
1695 colors = ["goldenrod", "white", "teal"]
1696 cmap = mcolors.ListedColormap(colors)
1697 norm = mcolors.BoundaryNorm(levels, cmap.N)
1698 logging.debug("Colourmap for %s.", cube.long_name)
1699 return cmap, levels, norm
1702def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1703 """Get a custom colorbar for a cube in the Beaufort Scale.
1705 Specific variable ranges can be separately set in user-supplied definition
1706 for x- or y-axis limits, or indicate where automated range preferred.
1708 Parameters
1709 ----------
1710 cube: Cube
1711 Cube of variable with Beaufort Scale in name.
1712 axis: "x", "y", optional
1713 Select the levels for just this axis of a line plot. The min and max
1714 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1715 setting a universal range is not desirable (e.g. temperature), users
1716 can set ymin/ymax values to "auto" in the colorbar definitions file.
1717 Where no additional xmin/xmax or ymin/ymax values are provided, the
1718 axis bounds default to use the vmin/vmax values provided.
1720 Returns
1721 -------
1722 cmap:
1723 Matplotlib colormap.
1724 levels:
1725 List of levels to use for plotting. For continuous plots the min and max
1726 should be taken as the range.
1727 norm:
1728 BoundaryNorm information.
1729 """
1730 if "difference" not in cube.long_name:
1731 if axis:
1732 levels = [0, 12]
1733 return None, levels, None
1734 else:
1735 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
1736 colors = [
1737 "black",
1738 (0, 0, 0.6),
1739 "blue",
1740 "cyan",
1741 "green",
1742 "yellow",
1743 (1, 0.5, 0),
1744 "red",
1745 "pink",
1746 "magenta",
1747 "purple",
1748 "maroon",
1749 "white",
1750 ]
1751 cmap = mcolors.ListedColormap(colors)
1752 norm = mcolors.BoundaryNorm(levels, cmap.N)
1753 logging.info("change colormap for Beaufort Scale colorbar.")
1754 return cmap, levels, norm
1755 else:
1756 if axis:
1757 levels = [-4, 4]
1758 return None, levels, None
1759 else:
1760 levels = [
1761 -3.5,
1762 -2.5,
1763 -1.5,
1764 -0.5,
1765 0.5,
1766 1.5,
1767 2.5,
1768 3.5,
1769 ]
1770 cmap = plt.get_cmap("bwr", 8)
1771 norm = mcolors.BoundaryNorm(levels, cmap.N)
1772 return cmap, levels, norm
1775def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
1776 """Return altered colourmap for temperature with change in units to Celsius.
1778 If "Celsius" appears anywhere in the name of a cube this function will be called.
1780 Parameters
1781 ----------
1782 cube: Cube
1783 Cube of variable for which the colorbar information is desired.
1784 cmap: Matplotlib colormap.
1785 levels: List
1786 List of levels to use for plotting. For continuous plots the min and max
1787 should be taken as the range.
1788 norm: BoundaryNorm.
1790 Returns
1791 -------
1792 cmap: Matplotlib colormap.
1793 levels: List
1794 List of levels to use for plotting. For continuous plots the min and max
1795 should be taken as the range.
1796 norm: BoundaryNorm.
1797 """
1798 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1799 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
1800 levels = np.array(levels)
1801 levels -= 273
1802 levels = levels.tolist()
1803 else:
1804 # Do nothing keep the existing colourbar attributes
1805 levels = levels
1806 cmap = cmap
1807 norm = norm
1808 return cmap, levels, norm
1811def _custom_colormap_probability(
1812 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
1813):
1814 """Get a custom colorbar for a probability cube.
1816 Specific variable ranges can be separately set in user-supplied definition
1817 for x- or y-axis limits, or indicate where automated range preferred.
1819 Parameters
1820 ----------
1821 cube: Cube
1822 Cube of variable with probability in name.
1823 axis: "x", "y", optional
1824 Select the levels for just this axis of a line plot. The min and max
1825 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1826 setting a universal range is not desirable (e.g. temperature), users
1827 can set ymin/ymax values to "auto" in the colorbar definitions file.
1828 Where no additional xmin/xmax or ymin/ymax values are provided, the
1829 axis bounds default to use the vmin/vmax values provided.
1831 Returns
1832 -------
1833 cmap:
1834 Matplotlib colormap.
1835 levels:
1836 List of levels to use for plotting. For continuous plots the min and max
1837 should be taken as the range.
1838 norm:
1839 BoundaryNorm information.
1840 """
1841 if axis:
1842 levels = [0, 1]
1843 return None, levels, None
1844 else:
1845 cmap = mcolors.ListedColormap(
1846 [
1847 "#FFFFFF",
1848 "#636363",
1849 "#e1dada",
1850 "#B5CAFF",
1851 "#8FB3FF",
1852 "#7F97FF",
1853 "#ABCF63",
1854 "#E8F59E",
1855 "#FFFA14",
1856 "#FFD121",
1857 "#FFA30A",
1858 ]
1859 )
1860 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]
1861 norm = mcolors.BoundaryNorm(levels, cmap.N)
1862 return cmap, levels, norm
1865def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
1866 """Return a custom colourmap for the current recipe."""
1867 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1868 if (
1869 any("surface_microphysical" in name for name in varnames)
1870 and "difference" not in cube.long_name
1871 and "mask" not in cube.long_name
1872 ):
1873 # Define the levels and colors
1874 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
1875 colors = [
1876 "w",
1877 (0, 0, 0.6),
1878 "b",
1879 "c",
1880 "g",
1881 "y",
1882 (1, 0.5, 0),
1883 "r",
1884 "pink",
1885 "m",
1886 "purple",
1887 "maroon",
1888 "gray",
1889 ]
1890 # Create a custom colormap
1891 cmap = mcolors.ListedColormap(colors)
1892 # Normalize the levels
1893 norm = mcolors.BoundaryNorm(levels, cmap.N)
1894 logging.info("change colormap for surface_microphysical variable colorbar.")
1895 else:
1896 # do nothing and keep existing colorbar attributes
1897 cmap = cmap
1898 levels = levels
1899 norm = norm
1900 return cmap, levels, norm
1903def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
1904 """Return custom colourmap for aviation colour state.
1906 If "aviation_colour_state" appears anywhere in the name of a cube
1907 this function will be called.
1909 Parameters
1910 ----------
1911 cube: Cube
1912 Cube of variable for which the colorbar information is desired.
1914 Returns
1915 -------
1916 cmap: Matplotlib colormap.
1917 levels: List
1918 List of levels to use for plotting. For continuous plots the min and max
1919 should be taken as the range.
1920 norm: BoundaryNorm.
1921 """
1922 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
1923 colors = [
1924 "#87ceeb",
1925 "#ffffff",
1926 "#8ced69",
1927 "#ffff00",
1928 "#ffd700",
1929 "#ffa500",
1930 "#fe3620",
1931 ]
1932 # Create a custom colormap
1933 cmap = mcolors.ListedColormap(colors)
1934 # Normalise the levels
1935 norm = mcolors.BoundaryNorm(levels, cmap.N)
1936 return cmap, levels, norm
1939def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
1940 """Return a custom colourmap for the current recipe."""
1941 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1942 if (
1943 any("visibility_in_air" in name for name in varnames)
1944 and "difference" not in cube.long_name
1945 and "mask" not in cube.long_name
1946 ):
1947 # Define the levels and colors (in km)
1948 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]
1949 norm = mcolors.BoundaryNorm(levels, cmap.N)
1950 colours = [
1951 "#8f00d6",
1952 "#d10000",
1953 "#ff9700",
1954 "#ffff00",
1955 "#00007f",
1956 "#6c9ccd",
1957 "#aae8ff",
1958 "#37a648",
1959 "#8edc64",
1960 "#c5ffc5",
1961 "#dcdcdc",
1962 "#ffffff",
1963 ]
1964 # Create a custom colormap
1965 cmap = mcolors.ListedColormap(colours)
1966 # Normalize the levels
1967 norm = mcolors.BoundaryNorm(levels, cmap.N)
1968 logging.info("change colormap for visibility_in_air variable colorbar.")
1969 else:
1970 # do nothing and keep existing colorbar attributes
1971 cmap = cmap
1972 levels = levels
1973 norm = norm
1974 return cmap, levels, norm
1977def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
1978 """Return number of models based on cube attributes."""
1979 model_names = list(
1980 filter(
1981 lambda x: x is not None,
1982 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
1983 )
1984 )
1985 if not model_names:
1986 logging.debug("Missing model names. Will assume single model.")
1987 return 1
1988 else:
1989 return len(model_names)
1992def _validate_cube_shape(
1993 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
1994) -> None:
1995 """Check all cubes have a model name."""
1996 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 1996 ↛ 1997line 1996 didn't jump to line 1997 because the condition on line 1996 was never true
1997 raise ValueError(
1998 f"The number of model names ({num_models}) should equal the number "
1999 f"of cubes ({len(cube)})."
2000 )
2003def _validate_cubes_coords(
2004 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2005) -> None:
2006 """Check same number of cubes as sequence coordinate for zip functions."""
2007 if len(cubes) != len(coords): 2007 ↛ 2008line 2007 didn't jump to line 2008 because the condition on line 2007 was never true
2008 raise ValueError(
2009 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2010 f"of sequence coordinates ({len(coords)})."
2011 f"Check that number of time entries in input data are consistent if "
2012 f"performing time-averaging steps prior to plotting outputs."
2013 )
2016####################
2017# Public functions #
2018####################
2021def spatial_contour_plot(
2022 cube: iris.cube.Cube,
2023 filename: str = None,
2024 sequence_coordinate: str = "time",
2025 stamp_coordinate: str = "realization",
2026 **kwargs,
2027) -> iris.cube.Cube:
2028 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2030 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2031 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2032 is present then postage stamp plots will be produced.
2034 Parameters
2035 ----------
2036 cube: Cube
2037 Iris cube of the data to plot. It should have two spatial dimensions,
2038 such as lat and lon, and may also have a another two dimension to be
2039 plotted sequentially and/or as postage stamp plots.
2040 filename: str, optional
2041 Name of the plot to write, used as a prefix for plot sequences. Defaults
2042 to the recipe name.
2043 sequence_coordinate: str, optional
2044 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2045 This coordinate must exist in the cube.
2046 stamp_coordinate: str, optional
2047 Coordinate about which to plot postage stamp plots. Defaults to
2048 ``"realization"``.
2050 Returns
2051 -------
2052 Cube
2053 The original cube (so further operations can be applied).
2055 Raises
2056 ------
2057 ValueError
2058 If the cube doesn't have the right dimensions.
2059 TypeError
2060 If the cube isn't a single cube.
2061 """
2062 _spatial_plot(
2063 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2064 )
2065 return cube
2068def spatial_pcolormesh_plot(
2069 cube: iris.cube.Cube,
2070 filename: str = None,
2071 sequence_coordinate: str = "time",
2072 stamp_coordinate: str = "realization",
2073 **kwargs,
2074) -> iris.cube.Cube:
2075 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2077 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2078 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2079 is present then postage stamp plots will be produced.
2081 This function is significantly faster than ``spatial_contour_plot``,
2082 especially at high resolutions, and should be preferred unless contiguous
2083 contour areas are important.
2085 Parameters
2086 ----------
2087 cube: Cube
2088 Iris cube of the data to plot. It should have two spatial dimensions,
2089 such as lat and lon, and may also have a another two dimension to be
2090 plotted sequentially and/or as postage stamp plots.
2091 filename: str, optional
2092 Name of the plot to write, used as a prefix for plot sequences. Defaults
2093 to the recipe name.
2094 sequence_coordinate: str, optional
2095 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2096 This coordinate must exist in the cube.
2097 stamp_coordinate: str, optional
2098 Coordinate about which to plot postage stamp plots. Defaults to
2099 ``"realization"``.
2101 Returns
2102 -------
2103 Cube
2104 The original cube (so further operations can be applied).
2106 Raises
2107 ------
2108 ValueError
2109 If the cube doesn't have the right dimensions.
2110 TypeError
2111 If the cube isn't a single cube.
2112 """
2113 _spatial_plot(
2114 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2115 )
2116 return cube
2119# TODO: Expand function to handle ensemble data.
2120# line_coordinate: str, optional
2121# Coordinate about which to plot multiple lines. Defaults to
2122# ``"realization"``.
2123def plot_line_series(
2124 cube: iris.cube.Cube | iris.cube.CubeList,
2125 filename: str = None,
2126 series_coordinate: str = "time",
2127 # line_coordinate: str = "realization",
2128 **kwargs,
2129) -> iris.cube.Cube | iris.cube.CubeList:
2130 """Plot a line plot for the specified coordinate.
2132 The Cube or CubeList must be 1D.
2134 Parameters
2135 ----------
2136 iris.cube | iris.cube.CubeList
2137 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2138 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2139 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2140 filename: str, optional
2141 Name of the plot to write, used as a prefix for plot sequences. Defaults
2142 to the recipe name.
2143 series_coordinate: str, optional
2144 Coordinate about which to make a series. Defaults to ``"time"``. This
2145 coordinate must exist in the cube.
2147 Returns
2148 -------
2149 iris.cube.Cube | iris.cube.CubeList
2150 The original Cube or CubeList (so further operations can be applied).
2151 plotted data.
2153 Raises
2154 ------
2155 ValueError
2156 If the cubes don't have the right dimensions.
2157 TypeError
2158 If the cube isn't a Cube or CubeList.
2159 """
2160 # Ensure we have a name for the plot file.
2161 title = get_recipe_metadata().get("title", "Untitled")
2163 if filename is None:
2164 filename = slugify(title)
2166 # Add file extension.
2167 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2169 num_models = _get_num_models(cube)
2171 _validate_cube_shape(cube, num_models)
2173 # Iterate over all cubes and extract coordinate to plot.
2174 cubes = iter_maybe(cube)
2175 coords = []
2176 for cube in cubes:
2177 try:
2178 coords.append(cube.coord(series_coordinate))
2179 except iris.exceptions.CoordinateNotFoundError as err:
2180 raise ValueError(
2181 f"Cube must have a {series_coordinate} coordinate."
2182 ) from err
2183 if cube.ndim > 2 or not cube.coords("realization"):
2184 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2186 # Do the actual plotting.
2187 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, title)
2189 # Add list of plots to plot metadata.
2190 plot_index = _append_to_plot_index([plot_filename])
2192 # Make a page to display the plots.
2193 _make_plot_html_page(plot_index)
2195 return cube
2198def plot_vertical_line_series(
2199 cubes: iris.cube.Cube | iris.cube.CubeList,
2200 filename: str = None,
2201 series_coordinate: str = "model_level_number",
2202 sequence_coordinate: str = "time",
2203 # line_coordinate: str = "realization",
2204 **kwargs,
2205) -> iris.cube.Cube | iris.cube.CubeList:
2206 """Plot a line plot against a type of vertical coordinate.
2208 The Cube or CubeList must be 1D.
2210 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2211 then a sequence of plots will be produced.
2213 Parameters
2214 ----------
2215 iris.cube | iris.cube.CubeList
2216 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2217 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2218 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2219 filename: str, optional
2220 Name of the plot to write, used as a prefix for plot sequences. Defaults
2221 to the recipe name.
2222 series_coordinate: str, optional
2223 Coordinate to plot on the y-axis. Can be ``pressure`` or
2224 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2225 for LFRic. Defaults to ``model_level_number``.
2226 This coordinate must exist in the cube.
2227 sequence_coordinate: str, optional
2228 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2229 This coordinate must exist in the cube.
2231 Returns
2232 -------
2233 iris.cube.Cube | iris.cube.CubeList
2234 The original Cube or CubeList (so further operations can be applied).
2235 Plotted data.
2237 Raises
2238 ------
2239 ValueError
2240 If the cubes doesn't have the right dimensions.
2241 TypeError
2242 If the cube isn't a Cube or CubeList.
2243 """
2244 # Ensure we have a name for the plot file.
2245 recipe_title = get_recipe_metadata().get("title", "Untitled")
2247 if filename is None:
2248 filename = slugify(recipe_title)
2250 cubes = iter_maybe(cubes)
2251 # Initialise empty list to hold all data from all cubes in a CubeList
2252 all_data = []
2254 # Store min/max ranges for x range.
2255 x_levels = []
2257 num_models = _get_num_models(cubes)
2259 _validate_cube_shape(cubes, num_models)
2261 # Iterate over all cubes in cube or CubeList and plot.
2262 coords = []
2263 for cube in cubes:
2264 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2265 try:
2266 coords.append(cube.coord(series_coordinate))
2267 except iris.exceptions.CoordinateNotFoundError as err:
2268 raise ValueError(
2269 f"Cube must have a {series_coordinate} coordinate."
2270 ) from err
2272 try:
2273 if cube.ndim > 1 or not cube.coords("realization"): 2273 ↛ 2281line 2273 didn't jump to line 2281 because the condition on line 2273 was always true
2274 cube.coord(sequence_coordinate)
2275 except iris.exceptions.CoordinateNotFoundError as err:
2276 raise ValueError(
2277 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2278 ) from err
2280 # Get minimum and maximum from levels information.
2281 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2282 if levels is not None: 2282 ↛ 2286line 2282 didn't jump to line 2286 because the condition on line 2282 was always true
2283 x_levels.append(min(levels))
2284 x_levels.append(max(levels))
2285 else:
2286 all_data.append(cube.data)
2288 if len(x_levels) == 0: 2288 ↛ 2290line 2288 didn't jump to line 2290 because the condition on line 2288 was never true
2289 # Combine all data into a single NumPy array
2290 combined_data = np.concatenate(all_data)
2292 # Set the lower and upper limit for the x-axis to ensure all plots have
2293 # same range. This needs to read the whole cube over the range of the
2294 # sequence and if applicable postage stamp coordinate.
2295 vmin = np.floor(combined_data.min())
2296 vmax = np.ceil(combined_data.max())
2297 else:
2298 vmin = min(x_levels)
2299 vmax = max(x_levels)
2301 # Matching the slices (matching by seq coord point; it may happen that
2302 # evaluated models do not cover the same seq coord range, hence matching
2303 # necessary)
2304 def filter_cube_iterables(cube_iterables) -> bool:
2305 return len(cube_iterables) == len(coords)
2307 cube_iterables = filter(
2308 filter_cube_iterables,
2309 (
2310 iris.cube.CubeList(
2311 s
2312 for s in itertools.chain.from_iterable(
2313 cb.slices_over(sequence_coordinate) for cb in cubes
2314 )
2315 if s.coord(sequence_coordinate).points[0] == point
2316 )
2317 for point in sorted(
2318 set(
2319 itertools.chain.from_iterable(
2320 cb.coord(sequence_coordinate).points for cb in cubes
2321 )
2322 )
2323 )
2324 ),
2325 )
2327 # Create a plot for each value of the sequence coordinate.
2328 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2329 # similar sequence values. Passing a CubeList into the internal plotting function
2330 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2331 # or an iris.cube.CubeList.
2332 plot_index = []
2333 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2334 for cubes_slice in cube_iterables:
2335 # Use sequence value so multiple sequences can merge.
2336 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2337 sequence_value = seq_coord.points[0]
2338 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2339 # Format the coordinate value in a unit appropriate way.
2340 title = f"{recipe_title}\n [{seq_coord.units.title(sequence_value)}]"
2341 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2342 if nplot == 1 and seq_coord.has_bounds: 2342 ↛ 2343line 2342 didn't jump to line 2343 because the condition on line 2342 was never true
2343 if np.size(seq_coord.bounds) > 1:
2344 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2345 # Do the actual plotting.
2346 _plot_and_save_vertical_line_series(
2347 cubes_slice,
2348 coords,
2349 "realization",
2350 plot_filename,
2351 series_coordinate,
2352 title=title,
2353 vmin=vmin,
2354 vmax=vmax,
2355 )
2356 plot_index.append(plot_filename)
2358 # Add list of plots to plot metadata.
2359 complete_plot_index = _append_to_plot_index(plot_index)
2361 # Make a page to display the plots.
2362 _make_plot_html_page(complete_plot_index)
2364 return cubes
2367def qq_plot(
2368 cubes: iris.cube.CubeList,
2369 coordinates: list[str],
2370 percentiles: list[float],
2371 model_names: list[str],
2372 filename: str = None,
2373 one_to_one: bool = True,
2374 **kwargs,
2375) -> iris.cube.CubeList:
2376 """Plot a Quantile-Quantile plot between two models for common time points.
2378 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2379 collapsed within the operator over all specified coordinates such as
2380 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2381 ensemble members to ensure a 1D cube (array).
2383 Parameters
2384 ----------
2385 cubes: iris.cube.CubeList
2386 Two cubes of the same variable with different models.
2387 coordinate: list[str]
2388 The list of coordinates to collapse over. This list should be
2389 every coordinate within the cube to result in a 1D cube around
2390 the percentile coordinate.
2391 percent: list[float]
2392 A list of percentiles to appear in the plot.
2393 model_names: list[str]
2394 A list of model names to appear on the axis of the plot.
2395 filename: str, optional
2396 Filename of the plot to write.
2397 one_to_one: bool, optional
2398 If True a 1:1 line is plotted; if False it is not. Default is True.
2400 Raises
2401 ------
2402 ValueError
2403 When the cubes are not compatible.
2405 Notes
2406 -----
2407 The quantile-quantile plot is a variant on the scatter plot representing
2408 two datasets by their quantiles (percentiles) for common time points.
2409 This plot does not use a theoretical distribution to compare against, but
2410 compares percentiles of two datasets. This plot does
2411 not use all raw data points, but plots the selected percentiles (quantiles) of
2412 each variable instead for the two datasets, thereby normalising the data for a
2413 direct comparison between the selected percentiles of the two dataset distributions.
2415 Quantile-quantile plots are valuable for comparing against
2416 observations and other models. Identical percentiles between the variables
2417 will lie on the one-to-one line implying the values correspond well to each
2418 other. Where there is a deviation from the one-to-one line a range of
2419 possibilities exist depending on how and where the data is shifted (e.g.,
2420 Wilks 2011 [Wilks2011]_).
2422 For distributions above the one-to-one line the distribution is left-skewed;
2423 below is right-skewed. A distinct break implies a bimodal distribution, and
2424 closer values/values further apart at the tails imply poor representation of
2425 the extremes.
2427 References
2428 ----------
2429 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2430 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2431 """
2432 # Check cubes using same functionality as the difference operator.
2433 if len(cubes) != 2:
2434 raise ValueError("cubes should contain exactly 2 cubes.")
2435 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2436 other: Cube = cubes.extract_cube(
2437 iris.Constraint(
2438 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2439 )
2440 )
2442 # Get spatial coord names.
2443 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2444 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2446 # Ensure cubes to compare are on common differencing grid.
2447 # This is triggered if either
2448 # i) latitude and longitude shapes are not the same. Note grid points
2449 # are not compared directly as these can differ through rounding
2450 # errors.
2451 # ii) or variables are known to often sit on different grid staggering
2452 # in different models (e.g. cell center vs cell edge), as is the case
2453 # for UM and LFRic comparisons.
2454 # In future greater choice of regridding method might be applied depending
2455 # on variable type. Linear regridding can in general be appropriate for smooth
2456 # variables. Care should be taken with interpretation of differences
2457 # given this dependency on regridding.
2458 if (
2459 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2460 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2461 ) or (
2462 base.long_name
2463 in [
2464 "eastward_wind_at_10m",
2465 "northward_wind_at_10m",
2466 "northward_wind_at_cell_centres",
2467 "eastward_wind_at_cell_centres",
2468 "zonal_wind_at_pressure_levels",
2469 "meridional_wind_at_pressure_levels",
2470 "potential_vorticity_at_pressure_levels",
2471 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2472 ]
2473 ):
2474 logging.debug(
2475 "Linear regridding base cube to other grid to compute differences"
2476 )
2477 base = regrid_onto_cube(base, other, method="Linear")
2479 # Extract just common time points.
2480 base, other = _extract_common_time_points(base, other)
2482 # Equalise attributes so we can merge.
2483 fully_equalise_attributes([base, other])
2484 logging.debug("Base: %s\nOther: %s", base, other)
2486 # Collapse cubes.
2487 base = collapse(
2488 base,
2489 coordinate=coordinates,
2490 method="PERCENTILE",
2491 additional_percent=percentiles,
2492 )
2493 other = collapse(
2494 other,
2495 coordinate=coordinates,
2496 method="PERCENTILE",
2497 additional_percent=percentiles,
2498 )
2500 # Ensure we have a name for the plot file.
2501 title = get_recipe_metadata().get("title", "Untitled")
2503 if filename is None:
2504 filename = slugify(title)
2506 # Add file extension.
2507 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2509 # Do the actual plotting on a scatter plot
2510 _plot_and_save_scatter_plot(
2511 base, other, plot_filename, title, one_to_one, model_names
2512 )
2514 # Add list of plots to plot metadata.
2515 plot_index = _append_to_plot_index([plot_filename])
2517 # Make a page to display the plots.
2518 _make_plot_html_page(plot_index)
2520 return iris.cube.CubeList([base, other])
2523def scatter_plot(
2524 cube_x: iris.cube.Cube | iris.cube.CubeList,
2525 cube_y: iris.cube.Cube | iris.cube.CubeList,
2526 filename: str = None,
2527 one_to_one: bool = True,
2528 **kwargs,
2529) -> iris.cube.CubeList:
2530 """Plot a scatter plot between two variables.
2532 Both cubes must be 1D.
2534 Parameters
2535 ----------
2536 cube_x: Cube | CubeList
2537 1 dimensional Cube of the data to plot on y-axis.
2538 cube_y: Cube | CubeList
2539 1 dimensional Cube of the data to plot on x-axis.
2540 filename: str, optional
2541 Filename of the plot to write.
2542 one_to_one: bool, optional
2543 If True a 1:1 line is plotted; if False it is not. Default is True.
2545 Returns
2546 -------
2547 cubes: CubeList
2548 CubeList of the original x and y cubes for further processing.
2550 Raises
2551 ------
2552 ValueError
2553 If the cube doesn't have the right dimensions and cubes not the same
2554 size.
2555 TypeError
2556 If the cube isn't a single cube.
2558 Notes
2559 -----
2560 Scatter plots are used for determining if there is a relationship between
2561 two variables. Positive relations have a slope going from bottom left to top
2562 right; Negative relations have a slope going from top left to bottom right.
2563 """
2564 # Iterate over all cubes in cube or CubeList and plot.
2565 for cube_iter in iter_maybe(cube_x):
2566 # Check cubes are correct shape.
2567 cube_iter = _check_single_cube(cube_iter)
2568 if cube_iter.ndim > 1:
2569 raise ValueError("cube_x must be 1D.")
2571 # Iterate over all cubes in cube or CubeList and plot.
2572 for cube_iter in iter_maybe(cube_y):
2573 # Check cubes are correct shape.
2574 cube_iter = _check_single_cube(cube_iter)
2575 if cube_iter.ndim > 1:
2576 raise ValueError("cube_y must be 1D.")
2578 # Ensure we have a name for the plot file.
2579 title = get_recipe_metadata().get("title", "Untitled")
2581 if filename is None:
2582 filename = slugify(title)
2584 # Add file extension.
2585 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2587 # Do the actual plotting.
2588 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2590 # Add list of plots to plot metadata.
2591 plot_index = _append_to_plot_index([plot_filename])
2593 # Make a page to display the plots.
2594 _make_plot_html_page(plot_index)
2596 return iris.cube.CubeList([cube_x, cube_y])
2599def vector_plot(
2600 cube_u: iris.cube.Cube,
2601 cube_v: iris.cube.Cube,
2602 filename: str = None,
2603 sequence_coordinate: str = "time",
2604 **kwargs,
2605) -> iris.cube.CubeList:
2606 """Plot a vector plot based on the input u and v components."""
2607 recipe_title = get_recipe_metadata().get("title", "Untitled")
2609 # Ensure we have a name for the plot file.
2610 if filename is None: 2610 ↛ 2611line 2610 didn't jump to line 2611 because the condition on line 2610 was never true
2611 filename = slugify(recipe_title)
2613 # Cubes must have a matching sequence coordinate.
2614 try:
2615 # Check that the u and v cubes have the same sequence coordinate.
2616 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2616 ↛ 2617line 2616 didn't jump to line 2617 because the condition on line 2616 was never true
2617 raise ValueError("Coordinates do not match.")
2618 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2619 raise ValueError(
2620 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2621 ) from err
2623 # Create a plot for each value of the sequence coordinate.
2624 plot_index = []
2625 for cube_u_slice, cube_v_slice in zip(
2626 cube_u.slices_over(sequence_coordinate),
2627 cube_v.slices_over(sequence_coordinate),
2628 strict=True,
2629 ):
2630 # Use sequence value so multiple sequences can merge.
2631 sequence_value = cube_u_slice.coord(sequence_coordinate).points[0]
2632 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2633 coord = cube_u_slice.coord(sequence_coordinate)
2634 # Format the coordinate value in a unit appropriate way.
2635 title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
2636 # Do the actual plotting.
2637 _plot_and_save_vector_plot(
2638 cube_u_slice,
2639 cube_v_slice,
2640 filename=plot_filename,
2641 title=title,
2642 method="contourf",
2643 )
2644 plot_index.append(plot_filename)
2646 # Add list of plots to plot metadata.
2647 complete_plot_index = _append_to_plot_index(plot_index)
2649 # Make a page to display the plots.
2650 _make_plot_html_page(complete_plot_index)
2652 return iris.cube.CubeList([cube_u, cube_v])
2655def plot_histogram_series(
2656 cubes: iris.cube.Cube | iris.cube.CubeList,
2657 filename: str = None,
2658 sequence_coordinate: str = "time",
2659 stamp_coordinate: str = "realization",
2660 single_plot: bool = False,
2661 **kwargs,
2662) -> iris.cube.Cube | iris.cube.CubeList:
2663 """Plot a histogram plot for each vertical level provided.
2665 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2666 is present then a sequence of plots will be produced using the time slider
2667 functionality to scroll through histograms against time. If a
2668 stamp_coordinate is present then postage stamp plots will be produced. If
2669 stamp_coordinate and single_plot is True, all postage stamp plots will be
2670 plotted in a single plot instead of separate postage stamp plots.
2672 Parameters
2673 ----------
2674 cubes: Cube | iris.cube.CubeList
2675 Iris cube or CubeList of the data to plot. It should have a single dimension other
2676 than the stamp coordinate.
2677 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2678 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2679 filename: str, optional
2680 Name of the plot to write, used as a prefix for plot sequences. Defaults
2681 to the recipe name.
2682 sequence_coordinate: str, optional
2683 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2684 This coordinate must exist in the cube and will be used for the time
2685 slider.
2686 stamp_coordinate: str, optional
2687 Coordinate about which to plot postage stamp plots. Defaults to
2688 ``"realization"``.
2689 single_plot: bool, optional
2690 If True, all postage stamp plots will be plotted in a single plot. If
2691 False, each postage stamp plot will be plotted separately. Is only valid
2692 if stamp_coordinate exists and has more than a single point.
2694 Returns
2695 -------
2696 iris.cube.Cube | iris.cube.CubeList
2697 The original Cube or CubeList (so further operations can be applied).
2698 Plotted data.
2700 Raises
2701 ------
2702 ValueError
2703 If the cube doesn't have the right dimensions.
2704 TypeError
2705 If the cube isn't a Cube or CubeList.
2706 """
2707 recipe_title = get_recipe_metadata().get("title", "Untitled")
2709 cubes = iter_maybe(cubes)
2711 # Ensure we have a name for the plot file.
2712 if filename is None:
2713 filename = slugify(recipe_title)
2715 # Internal plotting function.
2716 plotting_func = _plot_and_save_histogram_series
2718 num_models = _get_num_models(cubes)
2720 _validate_cube_shape(cubes, num_models)
2722 # If several histograms are plotted with time as sequence_coordinate for the
2723 # time slider option.
2724 for cube in cubes:
2725 try:
2726 cube.coord(sequence_coordinate)
2727 except iris.exceptions.CoordinateNotFoundError as err:
2728 raise ValueError(
2729 f"Cube must have a {sequence_coordinate} coordinate."
2730 ) from err
2732 # Get minimum and maximum from levels information.
2733 levels = None
2734 for cube in cubes: 2734 ↛ 2750line 2734 didn't jump to line 2750 because the loop on line 2734 didn't complete
2735 # First check if user-specified "auto" range variable.
2736 # This maintains the value of levels as None, so proceed.
2737 _, levels, _ = _colorbar_map_levels(cube, axis="y")
2738 if levels is None:
2739 break
2740 # If levels is changed, recheck to use the vmin,vmax or
2741 # levels-based ranges for histogram plots.
2742 _, levels, _ = _colorbar_map_levels(cube)
2743 logging.debug("levels: %s", levels)
2744 if levels is not None: 2744 ↛ 2734line 2744 didn't jump to line 2734 because the condition on line 2744 was always true
2745 vmin = min(levels)
2746 vmax = max(levels)
2747 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
2748 break
2750 if levels is None:
2751 vmin = min(cb.data.min() for cb in cubes)
2752 vmax = max(cb.data.max() for cb in cubes)
2754 # Make postage stamp plots if stamp_coordinate exists and has more than a
2755 # single point. If single_plot is True:
2756 # -- all postage stamp plots will be plotted in a single plot instead of
2757 # separate postage stamp plots.
2758 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2759 # produced per single model only
2760 if num_models == 1: 2760 ↛ 2773line 2760 didn't jump to line 2773 because the condition on line 2760 was always true
2761 if ( 2761 ↛ 2765line 2761 didn't jump to line 2765 because the condition on line 2761 was never true
2762 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2763 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2764 ):
2765 if single_plot:
2766 plotting_func = (
2767 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2768 )
2769 else:
2770 plotting_func = _plot_and_save_postage_stamp_histogram_series
2771 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2772 else:
2773 all_points = sorted(
2774 set(
2775 itertools.chain.from_iterable(
2776 cb.coord(sequence_coordinate).points for cb in cubes
2777 )
2778 )
2779 )
2780 all_slices = list(
2781 itertools.chain.from_iterable(
2782 cb.slices_over(sequence_coordinate) for cb in cubes
2783 )
2784 )
2785 # Matched slices (matched by seq coord point; it may happen that
2786 # evaluated models do not cover the same seq coord range, hence matching
2787 # necessary)
2788 cube_iterables = [
2789 iris.cube.CubeList(
2790 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2791 )
2792 for point in all_points
2793 ]
2795 plot_index = []
2796 nplot = np.size(cube.coord(sequence_coordinate).points)
2797 # Create a plot for each value of the sequence coordinate. Allowing for
2798 # multiple cubes in a CubeList to be plotted in the same plot for similar
2799 # sequence values. Passing a CubeList into the internal plotting function
2800 # for similar values of the sequence coordinate. cube_slice can be an
2801 # iris.cube.Cube or an iris.cube.CubeList.
2802 for cube_slice in cube_iterables:
2803 single_cube = cube_slice
2804 if isinstance(cube_slice, iris.cube.CubeList): 2804 ↛ 2805line 2804 didn't jump to line 2805 because the condition on line 2804 was never true
2805 single_cube = cube_slice[0]
2807 # Use sequence value so multiple sequences can merge.
2808 sequence_value = single_cube.coord(sequence_coordinate).points[0]
2809 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2810 coord = single_cube.coord(sequence_coordinate)
2811 # Format the coordinate value in a unit appropriate way.
2812 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
2813 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2814 if nplot == 1 and coord.has_bounds: 2814 ↛ 2815line 2814 didn't jump to line 2815 because the condition on line 2814 was never true
2815 if np.size(coord.bounds) > 1:
2816 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
2817 # Do the actual plotting.
2818 plotting_func(
2819 cube_slice,
2820 filename=plot_filename,
2821 stamp_coordinate=stamp_coordinate,
2822 title=title,
2823 vmin=vmin,
2824 vmax=vmax,
2825 )
2826 plot_index.append(plot_filename)
2828 # Add list of plots to plot metadata.
2829 complete_plot_index = _append_to_plot_index(plot_index)
2831 # Make a page to display the plots.
2832 _make_plot_html_page(complete_plot_index)
2834 return cubes
2837def plot_power_spectrum_series(
2838 cubes: iris.cube.Cube | iris.cube.CubeList,
2839 filename: str = None,
2840 sequence_coordinate: str = "time",
2841 stamp_coordinate: str = "realization",
2842 single_plot: bool = False,
2843 **kwargs,
2844) -> iris.cube.Cube | iris.cube.CubeList:
2845 """Plot a power spectrum plot for each vertical level provided.
2847 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
2848 is present then a sequence of plots will be produced using the time slider
2849 functionality to scroll through power spectrum against time. If a
2850 stamp_coordinate is present then postage stamp plots will be produced. If
2851 stamp_coordinate and single_plot is True, all postage stamp plots will be
2852 plotted in a single plot instead of separate postage stamp plots.
2854 Parameters
2855 ----------
2856 cubes: Cube | iris.cube.CubeList
2857 Iris cube or CubeList of the data to plot. It should have a single dimension other
2858 than the stamp coordinate.
2859 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2860 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2861 filename: str, optional
2862 Name of the plot to write, used as a prefix for plot sequences. Defaults
2863 to the recipe name.
2864 sequence_coordinate: str, optional
2865 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2866 This coordinate must exist in the cube and will be used for the time
2867 slider.
2868 stamp_coordinate: str, optional
2869 Coordinate about which to plot postage stamp plots. Defaults to
2870 ``"realization"``.
2871 single_plot: bool, optional
2872 If True, all postage stamp plots will be plotted in a single plot. If
2873 False, each postage stamp plot will be plotted separately. Is only valid
2874 if stamp_coordinate exists and has more than a single point.
2876 Returns
2877 -------
2878 iris.cube.Cube | iris.cube.CubeList
2879 The original Cube or CubeList (so further operations can be applied).
2880 Plotted data.
2882 Raises
2883 ------
2884 ValueError
2885 If the cube doesn't have the right dimensions.
2886 TypeError
2887 If the cube isn't a Cube or CubeList.
2888 """
2889 recipe_title = get_recipe_metadata().get("title", "Untitled")
2891 cubes = iter_maybe(cubes)
2892 # Ensure we have a name for the plot file.
2893 if filename is None:
2894 filename = slugify(recipe_title)
2896 # Internal plotting function.
2897 plotting_func = _plot_and_save_power_spectrum_series
2899 num_models = _get_num_models(cubes)
2901 _validate_cube_shape(cubes, num_models)
2903 # If several power spectra are plotted with time as sequence_coordinate for the
2904 # time slider option.
2905 for cube in cubes:
2906 try:
2907 cube.coord(sequence_coordinate)
2908 except iris.exceptions.CoordinateNotFoundError as err:
2909 raise ValueError(
2910 f"Cube must have a {sequence_coordinate} coordinate."
2911 ) from err
2913 # Make postage stamp plots if stamp_coordinate exists and has more than a
2914 # single point. If single_plot is True:
2915 # -- all postage stamp plots will be plotted in a single plot instead of
2916 # separate postage stamp plots.
2917 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2918 # produced per single model only
2919 if num_models == 1: 2919 ↛ 2932line 2919 didn't jump to line 2932 because the condition on line 2919 was always true
2920 if ( 2920 ↛ 2924line 2920 didn't jump to line 2924 because the condition on line 2920 was never true
2921 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2922 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2923 ):
2924 if single_plot:
2925 plotting_func = (
2926 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2927 )
2928 else:
2929 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2930 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2931 else:
2932 all_points = sorted(
2933 set(
2934 itertools.chain.from_iterable(
2935 cb.coord(sequence_coordinate).points for cb in cubes
2936 )
2937 )
2938 )
2939 all_slices = list(
2940 itertools.chain.from_iterable(
2941 cb.slices_over(sequence_coordinate) for cb in cubes
2942 )
2943 )
2944 # Matched slices (matched by seq coord point; it may happen that
2945 # evaluated models do not cover the same seq coord range, hence matching
2946 # necessary)
2947 cube_iterables = [
2948 iris.cube.CubeList(
2949 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2950 )
2951 for point in all_points
2952 ]
2954 plot_index = []
2955 nplot = np.size(cube.coord(sequence_coordinate).points)
2956 # Create a plot for each value of the sequence coordinate. Allowing for
2957 # multiple cubes in a CubeList to be plotted in the same plot for similar
2958 # sequence values. Passing a CubeList into the internal plotting function
2959 # for similar values of the sequence coordinate. cube_slice can be an
2960 # iris.cube.Cube or an iris.cube.CubeList.
2961 for cube_slice in cube_iterables:
2962 single_cube = cube_slice
2963 if isinstance(cube_slice, iris.cube.CubeList): 2963 ↛ 2964line 2963 didn't jump to line 2964 because the condition on line 2963 was never true
2964 single_cube = cube_slice[0]
2966 # Use sequence value so multiple sequences can merge.
2967 sequence_value = single_cube.coord(sequence_coordinate).points[0]
2968 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2969 coord = single_cube.coord(sequence_coordinate)
2970 # Format the coordinate value in a unit appropriate way.
2971 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
2972 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2973 if nplot == 1 and coord.has_bounds: 2973 ↛ 2977line 2973 didn't jump to line 2977 because the condition on line 2973 was always true
2974 if np.size(coord.bounds) > 1:
2975 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
2976 # Do the actual plotting.
2977 plotting_func(
2978 cube_slice,
2979 filename=plot_filename,
2980 stamp_coordinate=stamp_coordinate,
2981 title=title,
2982 )
2983 plot_index.append(plot_filename)
2985 # Add list of plots to plot metadata.
2986 complete_plot_index = _append_to_plot_index(plot_index)
2988 # Make a page to display the plots.
2989 _make_plot_html_page(complete_plot_index)
2991 return cubes
2994def _DCT_ps(y_3d):
2995 """Calculate power spectra for regional domains.
2997 Parameters
2998 ----------
2999 y_3d: 3D array
3000 3 dimensional array to calculate spectrum for.
3001 (2D field data with 3rd dimension of time)
3003 Returns: ps_array
3004 Array of power spectra values calculated for input field (for each time)
3006 Method for regional domains:
3007 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3008 as described in Denis et al 2002 [Denis_etal_2002]_.
3010 References
3011 ----------
3012 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3013 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3014 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3015 Monthly Weather Review, Vol. 130, 1812-1828
3016 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3017 """
3018 Nt, Ny, Nx = y_3d.shape
3020 # Max coefficient
3021 Nmin = min(Nx - 1, Ny - 1)
3023 # Create alpha matrix (of wavenumbers)
3024 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3026 # Prepare output array
3027 ps_array = np.zeros((Nt, Nmin))
3029 # Loop over time to get spectrum for each time.
3030 for t in range(Nt):
3031 y_2d = y_3d[t]
3033 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3034 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3035 # cosine basis functions at different spatial frequencies.
3037 # normalise spectrum to allow comparison between models.
3038 fkk = fft.dctn(y_2d, norm="ortho")
3040 # Normalise fkk
3041 fkk = fkk / np.sqrt(Ny * Nx)
3043 # calculate variance of spectral coefficient
3044 sigma_2 = fkk**2 / Nx / Ny
3046 # Group ellipses of alphas into the same wavenumber k/Nmin
3047 for k in range(1, Nmin + 1):
3048 alpha = k / Nmin
3049 alpha_p1 = (k + 1) / Nmin
3051 # Sum up elements matching k
3052 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3053 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3055 return ps_array
3058def _create_alpha_matrix(Ny, Nx):
3059 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3061 Parameters
3062 ----------
3063 Ny, Nx:
3064 Dimensions of the 2D field for which the power spectra is calculated. Used to
3065 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3066 single-scale parameter.
3068 Returns: alpha_matrix
3069 normalisation of 2D wavenumber axes, transforming the spectral domain into
3070 an elliptic coordinate system.
3072 """
3073 # Create x_indices: each row is [1, 2, ..., Nx]
3074 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3076 # Create y_indices: each column is [1, 2, ..., Ny]
3077 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3079 # Compute alpha_matrix
3080 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3082 return alpha_matrix