Coverage for src / CSET / operators / plot.py: 86%
962 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-16 13:48 +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
1497 ax = plt.gca()
1498 ax.plot(frequency, member.data)
1499 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1501 # Overall figure title.
1502 fig.suptitle(title, fontsize=16)
1504 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1505 logging.info("Saved power spectra postage stamp plot to %s", filename)
1506 plt.close(fig)
1509def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1510 cube: iris.cube.Cube,
1511 filename: str,
1512 title: str,
1513 stamp_coordinate: str,
1514 **kwargs,
1515):
1516 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1517 ax.set_title(title, fontsize=16)
1518 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1519 ax.set_ylabel("Power", fontsize=14)
1520 # Loop over all slices along the stamp_coordinate
1521 for member in cube.slices_over(stamp_coordinate):
1522 frequency = member.coord("frequency").points
1523 ax.plot(
1524 frequency,
1525 member.data,
1526 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1527 )
1529 # Add a legend
1530 ax.legend(fontsize=16)
1532 # Save the figure to a file
1533 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1535 # Close the figure
1536 plt.close(fig)
1539def _spatial_plot(
1540 method: Literal["contourf", "pcolormesh"],
1541 cube: iris.cube.Cube,
1542 filename: str | None,
1543 sequence_coordinate: str,
1544 stamp_coordinate: str,
1545 **kwargs,
1546):
1547 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1549 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1550 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1551 is present then postage stamp plots will be produced.
1553 Parameters
1554 ----------
1555 method: "contourf" | "pcolormesh"
1556 The plotting method to use.
1557 cube: Cube
1558 Iris cube of the data to plot. It should have two spatial dimensions,
1559 such as lat and lon, and may also have a another two dimension to be
1560 plotted sequentially and/or as postage stamp plots.
1561 filename: str | None
1562 Name of the plot to write, used as a prefix for plot sequences. If None
1563 uses the recipe name.
1564 sequence_coordinate: str
1565 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1566 This coordinate must exist in the cube.
1567 stamp_coordinate: str
1568 Coordinate about which to plot postage stamp plots. Defaults to
1569 ``"realization"``.
1571 Raises
1572 ------
1573 ValueError
1574 If the cube doesn't have the right dimensions.
1575 TypeError
1576 If the cube isn't a single cube.
1577 """
1578 recipe_title = get_recipe_metadata().get("title", "Untitled")
1580 # Ensure we have a name for the plot file.
1581 if filename is None:
1582 filename = slugify(recipe_title)
1584 # Ensure we've got a single cube.
1585 cube = _check_single_cube(cube)
1587 # Make postage stamp plots if stamp_coordinate exists and has more than a
1588 # single point.
1589 plotting_func = _plot_and_save_spatial_plot
1590 try:
1591 if cube.coord(stamp_coordinate).shape[0] > 1:
1592 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1593 except iris.exceptions.CoordinateNotFoundError:
1594 pass
1596 # Produce a geographical scatter plot if the data have a
1597 # dimension called observation or model_obs_error
1598 if any( 1598 ↛ 1602line 1598 didn't jump to line 1602 because the condition on line 1598 was never true
1599 crd.var_name == "station" or crd.var_name == "model_obs_error"
1600 for crd in cube.coords()
1601 ):
1602 plotting_func = _plot_and_save_scattermap_plot
1604 # Must have a sequence coordinate.
1605 try:
1606 cube.coord(sequence_coordinate)
1607 except iris.exceptions.CoordinateNotFoundError as err:
1608 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1610 # Create a plot for each value of the sequence coordinate.
1611 plot_index = []
1612 nplot = np.size(cube.coord(sequence_coordinate).points)
1613 for cube_slice in cube.slices_over(sequence_coordinate):
1614 # Use sequence value so multiple sequences can merge.
1615 sequence_value = cube_slice.coord(sequence_coordinate).points[0]
1616 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
1617 coord = cube_slice.coord(sequence_coordinate)
1618 # Format the coordinate value in a unit appropriate way.
1619 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
1620 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
1621 if nplot == 1 and coord.has_bounds:
1622 if np.size(coord.bounds) > 1:
1623 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
1624 # Do the actual plotting.
1625 plotting_func(
1626 cube_slice,
1627 filename=plot_filename,
1628 stamp_coordinate=stamp_coordinate,
1629 title=title,
1630 method=method,
1631 **kwargs,
1632 )
1633 plot_index.append(plot_filename)
1635 # Add list of plots to plot metadata.
1636 complete_plot_index = _append_to_plot_index(plot_index)
1638 # Make a page to display the plots.
1639 _make_plot_html_page(complete_plot_index)
1642def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1643 """Get colourmap for mask.
1645 If "mask_for_" appears anywhere in the name of a cube this function will be called
1646 regardless of the name of the variable to ensure a consistent plot.
1648 Parameters
1649 ----------
1650 cube: Cube
1651 Cube of variable for which the colorbar information is desired.
1652 axis: "x", "y", optional
1653 Select the levels for just this axis of a line plot. The min and max
1654 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1655 setting a universal range is not desirable (e.g. temperature), users
1656 can set ymin/ymax values to "auto" in the colorbar definitions file.
1657 Where no additional xmin/xmax or ymin/ymax values are provided, the
1658 axis bounds default to use the vmin/vmax values provided.
1660 Returns
1661 -------
1662 cmap:
1663 Matplotlib colormap.
1664 levels:
1665 List of levels to use for plotting. For continuous plots the min and max
1666 should be taken as the range.
1667 norm:
1668 BoundaryNorm information.
1669 """
1670 if "difference" not in cube.long_name:
1671 if axis:
1672 levels = [0, 1]
1673 # Complete settings based on levels.
1674 return None, levels, None
1675 else:
1676 # Define the levels and colors.
1677 levels = [0, 1, 2]
1678 colors = ["white", "dodgerblue"]
1679 # Create a custom color map.
1680 cmap = mcolors.ListedColormap(colors)
1681 # Normalize the levels.
1682 norm = mcolors.BoundaryNorm(levels, cmap.N)
1683 logging.debug("Colourmap for %s.", cube.long_name)
1684 return cmap, levels, norm
1685 else:
1686 if axis:
1687 levels = [-1, 1]
1688 return None, levels, None
1689 else:
1690 # Search for if mask difference, set to +/- 0.5 as values plotted <
1691 # not <=.
1692 levels = [-2, -0.5, 0.5, 2]
1693 colors = ["goldenrod", "white", "teal"]
1694 cmap = mcolors.ListedColormap(colors)
1695 norm = mcolors.BoundaryNorm(levels, cmap.N)
1696 logging.debug("Colourmap for %s.", cube.long_name)
1697 return cmap, levels, norm
1700def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1701 """Get a custom colorbar for a cube in the Beaufort Scale.
1703 Specific variable ranges can be separately set in user-supplied definition
1704 for x- or y-axis limits, or indicate where automated range preferred.
1706 Parameters
1707 ----------
1708 cube: Cube
1709 Cube of variable with Beaufort Scale in name.
1710 axis: "x", "y", optional
1711 Select the levels for just this axis of a line plot. The min and max
1712 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1713 setting a universal range is not desirable (e.g. temperature), users
1714 can set ymin/ymax values to "auto" in the colorbar definitions file.
1715 Where no additional xmin/xmax or ymin/ymax values are provided, the
1716 axis bounds default to use the vmin/vmax values provided.
1718 Returns
1719 -------
1720 cmap:
1721 Matplotlib colormap.
1722 levels:
1723 List of levels to use for plotting. For continuous plots the min and max
1724 should be taken as the range.
1725 norm:
1726 BoundaryNorm information.
1727 """
1728 if "difference" not in cube.long_name:
1729 if axis:
1730 levels = [0, 12]
1731 return None, levels, None
1732 else:
1733 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
1734 colors = [
1735 "black",
1736 (0, 0, 0.6),
1737 "blue",
1738 "cyan",
1739 "green",
1740 "yellow",
1741 (1, 0.5, 0),
1742 "red",
1743 "pink",
1744 "magenta",
1745 "purple",
1746 "maroon",
1747 "white",
1748 ]
1749 cmap = mcolors.ListedColormap(colors)
1750 norm = mcolors.BoundaryNorm(levels, cmap.N)
1751 logging.info("change colormap for Beaufort Scale colorbar.")
1752 return cmap, levels, norm
1753 else:
1754 if axis:
1755 levels = [-4, 4]
1756 return None, levels, None
1757 else:
1758 levels = [
1759 -3.5,
1760 -2.5,
1761 -1.5,
1762 -0.5,
1763 0.5,
1764 1.5,
1765 2.5,
1766 3.5,
1767 ]
1768 cmap = plt.get_cmap("bwr", 8)
1769 norm = mcolors.BoundaryNorm(levels, cmap.N)
1770 return cmap, levels, norm
1773def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
1774 """Return altered colourmap for temperature with change in units to Celsius.
1776 If "Celsius" appears anywhere in the name of a cube this function will be called.
1778 Parameters
1779 ----------
1780 cube: Cube
1781 Cube of variable for which the colorbar information is desired.
1782 cmap: Matplotlib colormap.
1783 levels: List
1784 List of levels to use for plotting. For continuous plots the min and max
1785 should be taken as the range.
1786 norm: BoundaryNorm.
1788 Returns
1789 -------
1790 cmap: Matplotlib colormap.
1791 levels: List
1792 List of levels to use for plotting. For continuous plots the min and max
1793 should be taken as the range.
1794 norm: BoundaryNorm.
1795 """
1796 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1797 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
1798 levels = np.array(levels)
1799 levels -= 273
1800 levels = levels.tolist()
1801 else:
1802 # Do nothing keep the existing colourbar attributes
1803 levels = levels
1804 cmap = cmap
1805 norm = norm
1806 return cmap, levels, norm
1809def _custom_colormap_probability(
1810 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
1811):
1812 """Get a custom colorbar for a probability cube.
1814 Specific variable ranges can be separately set in user-supplied definition
1815 for x- or y-axis limits, or indicate where automated range preferred.
1817 Parameters
1818 ----------
1819 cube: Cube
1820 Cube of variable with probability in name.
1821 axis: "x", "y", optional
1822 Select the levels for just this axis of a line plot. The min and max
1823 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1824 setting a universal range is not desirable (e.g. temperature), users
1825 can set ymin/ymax values to "auto" in the colorbar definitions file.
1826 Where no additional xmin/xmax or ymin/ymax values are provided, the
1827 axis bounds default to use the vmin/vmax values provided.
1829 Returns
1830 -------
1831 cmap:
1832 Matplotlib colormap.
1833 levels:
1834 List of levels to use for plotting. For continuous plots the min and max
1835 should be taken as the range.
1836 norm:
1837 BoundaryNorm information.
1838 """
1839 if axis:
1840 levels = [0, 1]
1841 return None, levels, None
1842 else:
1843 cmap = mcolors.ListedColormap(
1844 [
1845 "#FFFFFF",
1846 "#636363",
1847 "#e1dada",
1848 "#B5CAFF",
1849 "#8FB3FF",
1850 "#7F97FF",
1851 "#ABCF63",
1852 "#E8F59E",
1853 "#FFFA14",
1854 "#FFD121",
1855 "#FFA30A",
1856 ]
1857 )
1858 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]
1859 norm = mcolors.BoundaryNorm(levels, cmap.N)
1860 return cmap, levels, norm
1863def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
1864 """Return a custom colourmap for the current recipe."""
1865 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1866 if (
1867 any("surface_microphysical" in name for name in varnames)
1868 and "difference" not in cube.long_name
1869 and "mask" not in cube.long_name
1870 ):
1871 # Define the levels and colors
1872 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
1873 colors = [
1874 "w",
1875 (0, 0, 0.6),
1876 "b",
1877 "c",
1878 "g",
1879 "y",
1880 (1, 0.5, 0),
1881 "r",
1882 "pink",
1883 "m",
1884 "purple",
1885 "maroon",
1886 "gray",
1887 ]
1888 # Create a custom colormap
1889 cmap = mcolors.ListedColormap(colors)
1890 # Normalize the levels
1891 norm = mcolors.BoundaryNorm(levels, cmap.N)
1892 logging.info("change colormap for surface_microphysical variable colorbar.")
1893 else:
1894 # do nothing and keep existing colorbar attributes
1895 cmap = cmap
1896 levels = levels
1897 norm = norm
1898 return cmap, levels, norm
1901def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
1902 """Return custom colourmap for aviation colour state.
1904 If "aviation_colour_state" appears anywhere in the name of a cube
1905 this function will be called.
1907 Parameters
1908 ----------
1909 cube: Cube
1910 Cube of variable for which the colorbar information is desired.
1912 Returns
1913 -------
1914 cmap: Matplotlib colormap.
1915 levels: List
1916 List of levels to use for plotting. For continuous plots the min and max
1917 should be taken as the range.
1918 norm: BoundaryNorm.
1919 """
1920 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
1921 colors = [
1922 "#87ceeb",
1923 "#ffffff",
1924 "#8ced69",
1925 "#ffff00",
1926 "#ffd700",
1927 "#ffa500",
1928 "#fe3620",
1929 ]
1930 # Create a custom colormap
1931 cmap = mcolors.ListedColormap(colors)
1932 # Normalise the levels
1933 norm = mcolors.BoundaryNorm(levels, cmap.N)
1934 return cmap, levels, norm
1937def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
1938 """Return a custom colourmap for the current recipe."""
1939 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
1940 if (
1941 any("visibility_in_air" in name for name in varnames)
1942 and "difference" not in cube.long_name
1943 and "mask" not in cube.long_name
1944 ):
1945 # Define the levels and colors (in km)
1946 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]
1947 norm = mcolors.BoundaryNorm(levels, cmap.N)
1948 colours = [
1949 "#8f00d6",
1950 "#d10000",
1951 "#ff9700",
1952 "#ffff00",
1953 "#00007f",
1954 "#6c9ccd",
1955 "#aae8ff",
1956 "#37a648",
1957 "#8edc64",
1958 "#c5ffc5",
1959 "#dcdcdc",
1960 "#ffffff",
1961 ]
1962 # Create a custom colormap
1963 cmap = mcolors.ListedColormap(colours)
1964 # Normalize the levels
1965 norm = mcolors.BoundaryNorm(levels, cmap.N)
1966 logging.info("change colormap for visibility_in_air variable colorbar.")
1967 else:
1968 # do nothing and keep existing colorbar attributes
1969 cmap = cmap
1970 levels = levels
1971 norm = norm
1972 return cmap, levels, norm
1975def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
1976 """Return number of models based on cube attributes."""
1977 model_names = list(
1978 filter(
1979 lambda x: x is not None,
1980 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
1981 )
1982 )
1983 if not model_names:
1984 logging.debug("Missing model names. Will assume single model.")
1985 return 1
1986 else:
1987 return len(model_names)
1990def _validate_cube_shape(
1991 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
1992) -> None:
1993 """Check all cubes have a model name."""
1994 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 1994 ↛ 1995line 1994 didn't jump to line 1995 because the condition on line 1994 was never true
1995 raise ValueError(
1996 f"The number of model names ({num_models}) should equal the number "
1997 f"of cubes ({len(cube)})."
1998 )
2001def _validate_cubes_coords(
2002 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2003) -> None:
2004 """Check same number of cubes as sequence coordinate for zip functions."""
2005 if len(cubes) != len(coords): 2005 ↛ 2006line 2005 didn't jump to line 2006 because the condition on line 2005 was never true
2006 raise ValueError(
2007 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2008 f"of sequence coordinates ({len(coords)})."
2009 f"Check that number of time entries in input data are consistent if "
2010 f"performing time-averaging steps prior to plotting outputs."
2011 )
2014def _calculate_CFAD(
2015 cube: iris.cube.Cube, vertical_coordinate: str, bin_edges: list[float]
2016) -> iris.cube.Cube:
2017 """Calculate a Contour Frequency by Altitude Diagram (CFAD).
2019 Parameters
2020 ----------
2021 cube: iris.cube.Cube
2022 A cube of the data to be turned into a CFAD. It should be a minimum
2023 of two dimensions with one being a user specified vertical coordinate.
2024 vertical_coordinate: str
2025 The vertical coordinate of the cube for the CFAD to be calculated over.
2026 bin_edges: list[float]
2027 The bin edges for the histogram. The bins need to be specified to
2028 ensure consistency across the CFAD, otherwise it cannot be interpreted.
2030 Notes
2031 -----
2032 Contour Frequency by Altitude Diagrams (CFADs) were first designed by
2033 Yuter and Houze (1995)[YuterandHouze95]. They are calculated by binning the
2034 data by altitude and then by variable bins (e.g. temperature). The variable
2035 bins are then normalised by each altitude. This essenitally creates a
2036 normalised frequency distribution for each altitude. These are then stacked
2037 and combined in a single plot.
2039 References
2040 ----------
2041 .. [YuterandHouze95] Yuter S.E., and Houze, R.A. (1995) "Three-Dimensional
2042 Kinematic and Microphysical Evolution of Florida Cumulonimbus. Part II:
2043 Frequency Distributions of Vertical Velocity, Reflectivity, and
2044 Differential Reflectivity" Monthly Weather Review, vol. 123, 1941-1963,
2045 doi: 10.1175/1520-0493(1995)123<1941:TDKAME>2.0.CO;2
2046 """
2047 # Setup empty array for containing the CFAD data.
2048 CFAD_values = np.zeros(
2049 (len(cube.coord(vertical_coordinate).points), len(bin_edges) - 1)
2050 )
2052 # Calculate the CFAD as a histogram summing to one for each level.
2053 for i, level_cube in enumerate(cube.slices_over(vertical_coordinate)):
2054 # Note setting density to True does not produce the correct
2055 # normalization for a CFAD, where each row must sum to one.
2056 CFAD_values[i, :] = (
2057 np.histogram(level_cube.data.reshape(level_cube.data.size), bins=bin_edges)[
2058 0
2059 ]
2060 / level_cube.data.size
2061 )
2062 # Calculate central points for bins.
2063 bins = (np.array(bin_edges[:-1]) + np.array(bin_edges[1:])) / 2.0
2064 bin_bounds = np.array((bin_edges[:-1], bin_edges[1:])).T
2065 # Now construct the coordinates for the cube.
2066 vert_coord = cube.coord(vertical_coordinate)
2067 bin_coord = iris.coords.DimCoord(
2068 bins, bounds=bin_bounds, standard_name=cube.standard_name, units=cube.units
2069 )
2070 # Now construct the cube that is to be output.
2071 CFAD = iris.cube.Cube(
2072 CFAD_values,
2073 dim_coords_and_dims=[(vert_coord, 0), (bin_coord, 1)],
2074 long_name=f"{cube.name()}_cfad",
2075 units="1",
2076 )
2077 CFAD.rename(f"{cube.name()}_cfad")
2078 return CFAD
2081####################
2082# Public functions #
2083####################
2086def spatial_contour_plot(
2087 cube: iris.cube.Cube,
2088 filename: str = None,
2089 sequence_coordinate: str = "time",
2090 stamp_coordinate: str = "realization",
2091 **kwargs,
2092) -> iris.cube.Cube:
2093 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2095 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2096 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2097 is present then postage stamp plots will be produced.
2099 Parameters
2100 ----------
2101 cube: Cube
2102 Iris cube of the data to plot. It should have two spatial dimensions,
2103 such as lat and lon, and may also have a another two dimension to be
2104 plotted sequentially and/or as postage stamp plots.
2105 filename: str, optional
2106 Name of the plot to write, used as a prefix for plot sequences. Defaults
2107 to the recipe name.
2108 sequence_coordinate: str, optional
2109 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2110 This coordinate must exist in the cube.
2111 stamp_coordinate: str, optional
2112 Coordinate about which to plot postage stamp plots. Defaults to
2113 ``"realization"``.
2115 Returns
2116 -------
2117 Cube
2118 The original cube (so further operations can be applied).
2120 Raises
2121 ------
2122 ValueError
2123 If the cube doesn't have the right dimensions.
2124 TypeError
2125 If the cube isn't a single cube.
2126 """
2127 _spatial_plot(
2128 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2129 )
2130 return cube
2133def spatial_pcolormesh_plot(
2134 cube: iris.cube.Cube,
2135 filename: str = None,
2136 sequence_coordinate: str = "time",
2137 stamp_coordinate: str = "realization",
2138 **kwargs,
2139) -> iris.cube.Cube:
2140 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2142 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2143 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2144 is present then postage stamp plots will be produced.
2146 This function is significantly faster than ``spatial_contour_plot``,
2147 especially at high resolutions, and should be preferred unless contiguous
2148 contour areas are important.
2150 Parameters
2151 ----------
2152 cube: Cube
2153 Iris cube of the data to plot. It should have two spatial dimensions,
2154 such as lat and lon, and may also have a another two dimension to be
2155 plotted sequentially and/or as postage stamp plots.
2156 filename: str, optional
2157 Name of the plot to write, used as a prefix for plot sequences. Defaults
2158 to the recipe name.
2159 sequence_coordinate: str, optional
2160 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2161 This coordinate must exist in the cube.
2162 stamp_coordinate: str, optional
2163 Coordinate about which to plot postage stamp plots. Defaults to
2164 ``"realization"``.
2166 Returns
2167 -------
2168 Cube
2169 The original cube (so further operations can be applied).
2171 Raises
2172 ------
2173 ValueError
2174 If the cube doesn't have the right dimensions.
2175 TypeError
2176 If the cube isn't a single cube.
2177 """
2178 _spatial_plot(
2179 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2180 )
2181 return cube
2184# TODO: Expand function to handle ensemble data.
2185# line_coordinate: str, optional
2186# Coordinate about which to plot multiple lines. Defaults to
2187# ``"realization"``.
2188def plot_line_series(
2189 cube: iris.cube.Cube | iris.cube.CubeList,
2190 filename: str = None,
2191 series_coordinate: str = "time",
2192 # line_coordinate: str = "realization",
2193 **kwargs,
2194) -> iris.cube.Cube | iris.cube.CubeList:
2195 """Plot a line plot for the specified coordinate.
2197 The Cube or CubeList must be 1D.
2199 Parameters
2200 ----------
2201 iris.cube | iris.cube.CubeList
2202 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2203 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2204 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2205 filename: str, optional
2206 Name of the plot to write, used as a prefix for plot sequences. Defaults
2207 to the recipe name.
2208 series_coordinate: str, optional
2209 Coordinate about which to make a series. Defaults to ``"time"``. This
2210 coordinate must exist in the cube.
2212 Returns
2213 -------
2214 iris.cube.Cube | iris.cube.CubeList
2215 The original Cube or CubeList (so further operations can be applied).
2216 plotted data.
2218 Raises
2219 ------
2220 ValueError
2221 If the cubes don't have the right dimensions.
2222 TypeError
2223 If the cube isn't a Cube or CubeList.
2224 """
2225 # Ensure we have a name for the plot file.
2226 title = get_recipe_metadata().get("title", "Untitled")
2228 if filename is None:
2229 filename = slugify(title)
2231 # Add file extension.
2232 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2234 num_models = _get_num_models(cube)
2236 _validate_cube_shape(cube, num_models)
2238 # Iterate over all cubes and extract coordinate to plot.
2239 cubes = iter_maybe(cube)
2240 coords = []
2241 for cube in cubes:
2242 try:
2243 coords.append(cube.coord(series_coordinate))
2244 except iris.exceptions.CoordinateNotFoundError as err:
2245 raise ValueError(
2246 f"Cube must have a {series_coordinate} coordinate."
2247 ) from err
2248 if cube.ndim > 2 or not cube.coords("realization"):
2249 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2251 # Do the actual plotting.
2252 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, title)
2254 # Add list of plots to plot metadata.
2255 plot_index = _append_to_plot_index([plot_filename])
2257 # Make a page to display the plots.
2258 _make_plot_html_page(plot_index)
2260 return cube
2263def plot_vertical_line_series(
2264 cubes: iris.cube.Cube | iris.cube.CubeList,
2265 filename: str = None,
2266 series_coordinate: str = "model_level_number",
2267 sequence_coordinate: str = "time",
2268 # line_coordinate: str = "realization",
2269 **kwargs,
2270) -> iris.cube.Cube | iris.cube.CubeList:
2271 """Plot a line plot against a type of vertical coordinate.
2273 The Cube or CubeList must be 1D.
2275 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2276 then a sequence of plots will be produced.
2278 Parameters
2279 ----------
2280 iris.cube | iris.cube.CubeList
2281 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2282 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2283 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2284 filename: str, optional
2285 Name of the plot to write, used as a prefix for plot sequences. Defaults
2286 to the recipe name.
2287 series_coordinate: str, optional
2288 Coordinate to plot on the y-axis. Can be ``pressure`` or
2289 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2290 for LFRic. Defaults to ``model_level_number``.
2291 This coordinate must exist in the cube.
2292 sequence_coordinate: str, optional
2293 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2294 This coordinate must exist in the cube.
2296 Returns
2297 -------
2298 iris.cube.Cube | iris.cube.CubeList
2299 The original Cube or CubeList (so further operations can be applied).
2300 Plotted data.
2302 Raises
2303 ------
2304 ValueError
2305 If the cubes doesn't have the right dimensions.
2306 TypeError
2307 If the cube isn't a Cube or CubeList.
2308 """
2309 # Ensure we have a name for the plot file.
2310 recipe_title = get_recipe_metadata().get("title", "Untitled")
2312 if filename is None:
2313 filename = slugify(recipe_title)
2315 cubes = iter_maybe(cubes)
2316 # Initialise empty list to hold all data from all cubes in a CubeList
2317 all_data = []
2319 # Store min/max ranges for x range.
2320 x_levels = []
2322 num_models = _get_num_models(cubes)
2324 _validate_cube_shape(cubes, num_models)
2326 # Iterate over all cubes in cube or CubeList and plot.
2327 coords = []
2328 for cube in cubes:
2329 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2330 try:
2331 coords.append(cube.coord(series_coordinate))
2332 except iris.exceptions.CoordinateNotFoundError as err:
2333 raise ValueError(
2334 f"Cube must have a {series_coordinate} coordinate."
2335 ) from err
2337 try:
2338 if cube.ndim > 1 or not cube.coords("realization"): 2338 ↛ 2346line 2338 didn't jump to line 2346 because the condition on line 2338 was always true
2339 cube.coord(sequence_coordinate)
2340 except iris.exceptions.CoordinateNotFoundError as err:
2341 raise ValueError(
2342 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2343 ) from err
2345 # Get minimum and maximum from levels information.
2346 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2347 if levels is not None: 2347 ↛ 2351line 2347 didn't jump to line 2351 because the condition on line 2347 was always true
2348 x_levels.append(min(levels))
2349 x_levels.append(max(levels))
2350 else:
2351 all_data.append(cube.data)
2353 if len(x_levels) == 0: 2353 ↛ 2355line 2353 didn't jump to line 2355 because the condition on line 2353 was never true
2354 # Combine all data into a single NumPy array
2355 combined_data = np.concatenate(all_data)
2357 # Set the lower and upper limit for the x-axis to ensure all plots have
2358 # same range. This needs to read the whole cube over the range of the
2359 # sequence and if applicable postage stamp coordinate.
2360 vmin = np.floor(combined_data.min())
2361 vmax = np.ceil(combined_data.max())
2362 else:
2363 vmin = min(x_levels)
2364 vmax = max(x_levels)
2366 # Matching the slices (matching by seq coord point; it may happen that
2367 # evaluated models do not cover the same seq coord range, hence matching
2368 # necessary)
2369 def filter_cube_iterables(cube_iterables) -> bool:
2370 return len(cube_iterables) == len(coords)
2372 cube_iterables = filter(
2373 filter_cube_iterables,
2374 (
2375 iris.cube.CubeList(
2376 s
2377 for s in itertools.chain.from_iterable(
2378 cb.slices_over(sequence_coordinate) for cb in cubes
2379 )
2380 if s.coord(sequence_coordinate).points[0] == point
2381 )
2382 for point in sorted(
2383 set(
2384 itertools.chain.from_iterable(
2385 cb.coord(sequence_coordinate).points for cb in cubes
2386 )
2387 )
2388 )
2389 ),
2390 )
2392 # Create a plot for each value of the sequence coordinate.
2393 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2394 # similar sequence values. Passing a CubeList into the internal plotting function
2395 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2396 # or an iris.cube.CubeList.
2397 plot_index = []
2398 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2399 for cubes_slice in cube_iterables:
2400 # Use sequence value so multiple sequences can merge.
2401 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2402 sequence_value = seq_coord.points[0]
2403 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2404 # Format the coordinate value in a unit appropriate way.
2405 title = f"{recipe_title}\n [{seq_coord.units.title(sequence_value)}]"
2406 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2407 if nplot == 1 and seq_coord.has_bounds: 2407 ↛ 2408line 2407 didn't jump to line 2408 because the condition on line 2407 was never true
2408 if np.size(seq_coord.bounds) > 1:
2409 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2410 # Do the actual plotting.
2411 _plot_and_save_vertical_line_series(
2412 cubes_slice,
2413 coords,
2414 "realization",
2415 plot_filename,
2416 series_coordinate,
2417 title=title,
2418 vmin=vmin,
2419 vmax=vmax,
2420 )
2421 plot_index.append(plot_filename)
2423 # Add list of plots to plot metadata.
2424 complete_plot_index = _append_to_plot_index(plot_index)
2426 # Make a page to display the plots.
2427 _make_plot_html_page(complete_plot_index)
2429 return cubes
2432def qq_plot(
2433 cubes: iris.cube.CubeList,
2434 coordinates: list[str],
2435 percentiles: list[float],
2436 model_names: list[str],
2437 filename: str = None,
2438 one_to_one: bool = True,
2439 **kwargs,
2440) -> iris.cube.CubeList:
2441 """Plot a Quantile-Quantile plot between two models for common time points.
2443 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2444 collapsed within the operator over all specified coordinates such as
2445 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2446 ensemble members to ensure a 1D cube (array).
2448 Parameters
2449 ----------
2450 cubes: iris.cube.CubeList
2451 Two cubes of the same variable with different models.
2452 coordinate: list[str]
2453 The list of coordinates to collapse over. This list should be
2454 every coordinate within the cube to result in a 1D cube around
2455 the percentile coordinate.
2456 percent: list[float]
2457 A list of percentiles to appear in the plot.
2458 model_names: list[str]
2459 A list of model names to appear on the axis of the plot.
2460 filename: str, optional
2461 Filename of the plot to write.
2462 one_to_one: bool, optional
2463 If True a 1:1 line is plotted; if False it is not. Default is True.
2465 Raises
2466 ------
2467 ValueError
2468 When the cubes are not compatible.
2470 Notes
2471 -----
2472 The quantile-quantile plot is a variant on the scatter plot representing
2473 two datasets by their quantiles (percentiles) for common time points.
2474 This plot does not use a theoretical distribution to compare against, but
2475 compares percentiles of two datasets. This plot does
2476 not use all raw data points, but plots the selected percentiles (quantiles) of
2477 each variable instead for the two datasets, thereby normalising the data for a
2478 direct comparison between the selected percentiles of the two dataset distributions.
2480 Quantile-quantile plots are valuable for comparing against
2481 observations and other models. Identical percentiles between the variables
2482 will lie on the one-to-one line implying the values correspond well to each
2483 other. Where there is a deviation from the one-to-one line a range of
2484 possibilities exist depending on how and where the data is shifted (e.g.,
2485 Wilks 2011 [Wilks2011]_).
2487 For distributions above the one-to-one line the distribution is left-skewed;
2488 below is right-skewed. A distinct break implies a bimodal distribution, and
2489 closer values/values further apart at the tails imply poor representation of
2490 the extremes.
2492 References
2493 ----------
2494 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2495 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2496 """
2497 # Check cubes using same functionality as the difference operator.
2498 if len(cubes) != 2:
2499 raise ValueError("cubes should contain exactly 2 cubes.")
2500 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2501 other: Cube = cubes.extract_cube(
2502 iris.Constraint(
2503 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2504 )
2505 )
2507 # Get spatial coord names.
2508 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2509 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2511 # Ensure cubes to compare are on common differencing grid.
2512 # This is triggered if either
2513 # i) latitude and longitude shapes are not the same. Note grid points
2514 # are not compared directly as these can differ through rounding
2515 # errors.
2516 # ii) or variables are known to often sit on different grid staggering
2517 # in different models (e.g. cell center vs cell edge), as is the case
2518 # for UM and LFRic comparisons.
2519 # In future greater choice of regridding method might be applied depending
2520 # on variable type. Linear regridding can in general be appropriate for smooth
2521 # variables. Care should be taken with interpretation of differences
2522 # given this dependency on regridding.
2523 if (
2524 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2525 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2526 ) or (
2527 base.long_name
2528 in [
2529 "eastward_wind_at_10m",
2530 "northward_wind_at_10m",
2531 "northward_wind_at_cell_centres",
2532 "eastward_wind_at_cell_centres",
2533 "zonal_wind_at_pressure_levels",
2534 "meridional_wind_at_pressure_levels",
2535 "potential_vorticity_at_pressure_levels",
2536 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2537 ]
2538 ):
2539 logging.debug(
2540 "Linear regridding base cube to other grid to compute differences"
2541 )
2542 base = regrid_onto_cube(base, other, method="Linear")
2544 # Extract just common time points.
2545 base, other = _extract_common_time_points(base, other)
2547 # Equalise attributes so we can merge.
2548 fully_equalise_attributes([base, other])
2549 logging.debug("Base: %s\nOther: %s", base, other)
2551 # Collapse cubes.
2552 base = collapse(
2553 base,
2554 coordinate=coordinates,
2555 method="PERCENTILE",
2556 additional_percent=percentiles,
2557 )
2558 other = collapse(
2559 other,
2560 coordinate=coordinates,
2561 method="PERCENTILE",
2562 additional_percent=percentiles,
2563 )
2565 # Ensure we have a name for the plot file.
2566 title = get_recipe_metadata().get("title", "Untitled")
2568 if filename is None:
2569 filename = slugify(title)
2571 # Add file extension.
2572 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2574 # Do the actual plotting on a scatter plot
2575 _plot_and_save_scatter_plot(
2576 base, other, plot_filename, title, one_to_one, model_names
2577 )
2579 # Add list of plots to plot metadata.
2580 plot_index = _append_to_plot_index([plot_filename])
2582 # Make a page to display the plots.
2583 _make_plot_html_page(plot_index)
2585 return iris.cube.CubeList([base, other])
2588def scatter_plot(
2589 cube_x: iris.cube.Cube | iris.cube.CubeList,
2590 cube_y: iris.cube.Cube | iris.cube.CubeList,
2591 filename: str = None,
2592 one_to_one: bool = True,
2593 **kwargs,
2594) -> iris.cube.CubeList:
2595 """Plot a scatter plot between two variables.
2597 Both cubes must be 1D.
2599 Parameters
2600 ----------
2601 cube_x: Cube | CubeList
2602 1 dimensional Cube of the data to plot on y-axis.
2603 cube_y: Cube | CubeList
2604 1 dimensional Cube of the data to plot on x-axis.
2605 filename: str, optional
2606 Filename of the plot to write.
2607 one_to_one: bool, optional
2608 If True a 1:1 line is plotted; if False it is not. Default is True.
2610 Returns
2611 -------
2612 cubes: CubeList
2613 CubeList of the original x and y cubes for further processing.
2615 Raises
2616 ------
2617 ValueError
2618 If the cube doesn't have the right dimensions and cubes not the same
2619 size.
2620 TypeError
2621 If the cube isn't a single cube.
2623 Notes
2624 -----
2625 Scatter plots are used for determining if there is a relationship between
2626 two variables. Positive relations have a slope going from bottom left to top
2627 right; Negative relations have a slope going from top left to bottom right.
2628 """
2629 # Iterate over all cubes in cube or CubeList and plot.
2630 for cube_iter in iter_maybe(cube_x):
2631 # Check cubes are correct shape.
2632 cube_iter = _check_single_cube(cube_iter)
2633 if cube_iter.ndim > 1:
2634 raise ValueError("cube_x must be 1D.")
2636 # Iterate over all cubes in cube or CubeList and plot.
2637 for cube_iter in iter_maybe(cube_y):
2638 # Check cubes are correct shape.
2639 cube_iter = _check_single_cube(cube_iter)
2640 if cube_iter.ndim > 1:
2641 raise ValueError("cube_y must be 1D.")
2643 # Ensure we have a name for the plot file.
2644 title = get_recipe_metadata().get("title", "Untitled")
2646 if filename is None:
2647 filename = slugify(title)
2649 # Add file extension.
2650 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2652 # Do the actual plotting.
2653 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2655 # Add list of plots to plot metadata.
2656 plot_index = _append_to_plot_index([plot_filename])
2658 # Make a page to display the plots.
2659 _make_plot_html_page(plot_index)
2661 return iris.cube.CubeList([cube_x, cube_y])
2664def vector_plot(
2665 cube_u: iris.cube.Cube,
2666 cube_v: iris.cube.Cube,
2667 filename: str = None,
2668 sequence_coordinate: str = "time",
2669 **kwargs,
2670) -> iris.cube.CubeList:
2671 """Plot a vector plot based on the input u and v components."""
2672 recipe_title = get_recipe_metadata().get("title", "Untitled")
2674 # Ensure we have a name for the plot file.
2675 if filename is None: 2675 ↛ 2676line 2675 didn't jump to line 2676 because the condition on line 2675 was never true
2676 filename = slugify(recipe_title)
2678 # Cubes must have a matching sequence coordinate.
2679 try:
2680 # Check that the u and v cubes have the same sequence coordinate.
2681 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2681 ↛ 2682line 2681 didn't jump to line 2682 because the condition on line 2681 was never true
2682 raise ValueError("Coordinates do not match.")
2683 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2684 raise ValueError(
2685 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2686 ) from err
2688 # Create a plot for each value of the sequence coordinate.
2689 plot_index = []
2690 for cube_u_slice, cube_v_slice in zip(
2691 cube_u.slices_over(sequence_coordinate),
2692 cube_v.slices_over(sequence_coordinate),
2693 strict=True,
2694 ):
2695 # Use sequence value so multiple sequences can merge.
2696 sequence_value = cube_u_slice.coord(sequence_coordinate).points[0]
2697 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2698 coord = cube_u_slice.coord(sequence_coordinate)
2699 # Format the coordinate value in a unit appropriate way.
2700 title = f"{recipe_title}\n{coord.units.title(coord.points[0])}"
2701 # Do the actual plotting.
2702 _plot_and_save_vector_plot(
2703 cube_u_slice,
2704 cube_v_slice,
2705 filename=plot_filename,
2706 title=title,
2707 method="contourf",
2708 )
2709 plot_index.append(plot_filename)
2711 # Add list of plots to plot metadata.
2712 complete_plot_index = _append_to_plot_index(plot_index)
2714 # Make a page to display the plots.
2715 _make_plot_html_page(complete_plot_index)
2717 return iris.cube.CubeList([cube_u, cube_v])
2720def plot_histogram_series(
2721 cubes: iris.cube.Cube | iris.cube.CubeList,
2722 filename: str = None,
2723 sequence_coordinate: str = "time",
2724 stamp_coordinate: str = "realization",
2725 single_plot: bool = False,
2726 **kwargs,
2727) -> iris.cube.Cube | iris.cube.CubeList:
2728 """Plot a histogram plot for each vertical level provided.
2730 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2731 is present then a sequence of plots will be produced using the time slider
2732 functionality to scroll through histograms against time. If a
2733 stamp_coordinate is present then postage stamp plots will be produced. If
2734 stamp_coordinate and single_plot is True, all postage stamp plots will be
2735 plotted in a single plot instead of separate postage stamp plots.
2737 Parameters
2738 ----------
2739 cubes: Cube | iris.cube.CubeList
2740 Iris cube or CubeList of the data to plot. It should have a single dimension other
2741 than the stamp coordinate.
2742 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2743 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2744 filename: str, optional
2745 Name of the plot to write, used as a prefix for plot sequences. Defaults
2746 to the recipe name.
2747 sequence_coordinate: str, optional
2748 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2749 This coordinate must exist in the cube and will be used for the time
2750 slider.
2751 stamp_coordinate: str, optional
2752 Coordinate about which to plot postage stamp plots. Defaults to
2753 ``"realization"``.
2754 single_plot: bool, optional
2755 If True, all postage stamp plots will be plotted in a single plot. If
2756 False, each postage stamp plot will be plotted separately. Is only valid
2757 if stamp_coordinate exists and has more than a single point.
2759 Returns
2760 -------
2761 iris.cube.Cube | iris.cube.CubeList
2762 The original Cube or CubeList (so further operations can be applied).
2763 Plotted data.
2765 Raises
2766 ------
2767 ValueError
2768 If the cube doesn't have the right dimensions.
2769 TypeError
2770 If the cube isn't a Cube or CubeList.
2771 """
2772 recipe_title = get_recipe_metadata().get("title", "Untitled")
2774 cubes = iter_maybe(cubes)
2776 # Ensure we have a name for the plot file.
2777 if filename is None:
2778 filename = slugify(recipe_title)
2780 # Internal plotting function.
2781 plotting_func = _plot_and_save_histogram_series
2783 num_models = _get_num_models(cubes)
2785 _validate_cube_shape(cubes, num_models)
2787 # If several histograms are plotted with time as sequence_coordinate for the
2788 # time slider option.
2789 for cube in cubes:
2790 try:
2791 cube.coord(sequence_coordinate)
2792 except iris.exceptions.CoordinateNotFoundError as err:
2793 raise ValueError(
2794 f"Cube must have a {sequence_coordinate} coordinate."
2795 ) from err
2797 # Get minimum and maximum from levels information.
2798 levels = None
2799 for cube in cubes: 2799 ↛ 2815line 2799 didn't jump to line 2815 because the loop on line 2799 didn't complete
2800 # First check if user-specified "auto" range variable.
2801 # This maintains the value of levels as None, so proceed.
2802 _, levels, _ = _colorbar_map_levels(cube, axis="y")
2803 if levels is None:
2804 break
2805 # If levels is changed, recheck to use the vmin,vmax or
2806 # levels-based ranges for histogram plots.
2807 _, levels, _ = _colorbar_map_levels(cube)
2808 logging.debug("levels: %s", levels)
2809 if levels is not None: 2809 ↛ 2799line 2809 didn't jump to line 2799 because the condition on line 2809 was always true
2810 vmin = min(levels)
2811 vmax = max(levels)
2812 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
2813 break
2815 if levels is None:
2816 vmin = min(cb.data.min() for cb in cubes)
2817 vmax = max(cb.data.max() for cb in cubes)
2819 # Make postage stamp plots if stamp_coordinate exists and has more than a
2820 # single point. If single_plot is True:
2821 # -- all postage stamp plots will be plotted in a single plot instead of
2822 # separate postage stamp plots.
2823 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2824 # produced per single model only
2825 if num_models == 1: 2825 ↛ 2838line 2825 didn't jump to line 2838 because the condition on line 2825 was always true
2826 if ( 2826 ↛ 2830line 2826 didn't jump to line 2830 because the condition on line 2826 was never true
2827 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2828 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2829 ):
2830 if single_plot:
2831 plotting_func = (
2832 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2833 )
2834 else:
2835 plotting_func = _plot_and_save_postage_stamp_histogram_series
2836 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2837 else:
2838 all_points = sorted(
2839 set(
2840 itertools.chain.from_iterable(
2841 cb.coord(sequence_coordinate).points for cb in cubes
2842 )
2843 )
2844 )
2845 all_slices = list(
2846 itertools.chain.from_iterable(
2847 cb.slices_over(sequence_coordinate) for cb in cubes
2848 )
2849 )
2850 # Matched slices (matched by seq coord point; it may happen that
2851 # evaluated models do not cover the same seq coord range, hence matching
2852 # necessary)
2853 cube_iterables = [
2854 iris.cube.CubeList(
2855 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2856 )
2857 for point in all_points
2858 ]
2860 plot_index = []
2861 nplot = np.size(cube.coord(sequence_coordinate).points)
2862 # Create a plot for each value of the sequence coordinate. Allowing for
2863 # multiple cubes in a CubeList to be plotted in the same plot for similar
2864 # sequence values. Passing a CubeList into the internal plotting function
2865 # for similar values of the sequence coordinate. cube_slice can be an
2866 # iris.cube.Cube or an iris.cube.CubeList.
2867 for cube_slice in cube_iterables:
2868 single_cube = cube_slice
2869 if isinstance(cube_slice, iris.cube.CubeList): 2869 ↛ 2870line 2869 didn't jump to line 2870 because the condition on line 2869 was never true
2870 single_cube = cube_slice[0]
2872 # Use sequence value so multiple sequences can merge.
2873 sequence_value = single_cube.coord(sequence_coordinate).points[0]
2874 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
2875 coord = single_cube.coord(sequence_coordinate)
2876 # Format the coordinate value in a unit appropriate way.
2877 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
2878 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2879 if nplot == 1 and coord.has_bounds: 2879 ↛ 2880line 2879 didn't jump to line 2880 because the condition on line 2879 was never true
2880 if np.size(coord.bounds) > 1:
2881 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
2882 # Do the actual plotting.
2883 plotting_func(
2884 cube_slice,
2885 filename=plot_filename,
2886 stamp_coordinate=stamp_coordinate,
2887 title=title,
2888 vmin=vmin,
2889 vmax=vmax,
2890 )
2891 plot_index.append(plot_filename)
2893 # Add list of plots to plot metadata.
2894 complete_plot_index = _append_to_plot_index(plot_index)
2896 # Make a page to display the plots.
2897 _make_plot_html_page(complete_plot_index)
2899 return cubes
2902def plot_power_spectrum_series(
2903 cubes: iris.cube.Cube | iris.cube.CubeList,
2904 filename: str = None,
2905 sequence_coordinate: str = "time",
2906 stamp_coordinate: str = "realization",
2907 single_plot: bool = False,
2908 **kwargs,
2909) -> iris.cube.Cube | iris.cube.CubeList:
2910 """Plot a power spectrum plot for each vertical level provided.
2912 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
2913 is present then a sequence of plots will be produced using the time slider
2914 functionality to scroll through power spectrum against time. If a
2915 stamp_coordinate is present then postage stamp plots will be produced. If
2916 stamp_coordinate and single_plot is True, all postage stamp plots will be
2917 plotted in a single plot instead of separate postage stamp plots.
2919 Parameters
2920 ----------
2921 cubes: Cube | iris.cube.CubeList
2922 Iris cube or CubeList of the data to plot. It should have a single dimension other
2923 than the stamp coordinate.
2924 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2925 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2926 filename: str, optional
2927 Name of the plot to write, used as a prefix for plot sequences. Defaults
2928 to the recipe name.
2929 sequence_coordinate: str, optional
2930 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2931 This coordinate must exist in the cube and will be used for the time
2932 slider.
2933 stamp_coordinate: str, optional
2934 Coordinate about which to plot postage stamp plots. Defaults to
2935 ``"realization"``.
2936 single_plot: bool, optional
2937 If True, all postage stamp plots will be plotted in a single plot. If
2938 False, each postage stamp plot will be plotted separately. Is only valid
2939 if stamp_coordinate exists and has more than a single point.
2941 Returns
2942 -------
2943 iris.cube.Cube | iris.cube.CubeList
2944 The original Cube or CubeList (so further operations can be applied).
2945 Plotted data.
2947 Raises
2948 ------
2949 ValueError
2950 If the cube doesn't have the right dimensions.
2951 TypeError
2952 If the cube isn't a Cube or CubeList.
2953 """
2954 recipe_title = get_recipe_metadata().get("title", "Untitled")
2956 cubes = iter_maybe(cubes)
2957 # Ensure we have a name for the plot file.
2958 if filename is None:
2959 filename = slugify(recipe_title)
2961 # Internal plotting function.
2962 plotting_func = _plot_and_save_power_spectrum_series
2964 num_models = _get_num_models(cubes)
2966 _validate_cube_shape(cubes, num_models)
2968 # If several power spectra are plotted with time as sequence_coordinate for the
2969 # time slider option.
2970 for cube in cubes:
2971 try:
2972 cube.coord(sequence_coordinate)
2973 except iris.exceptions.CoordinateNotFoundError as err:
2974 raise ValueError(
2975 f"Cube must have a {sequence_coordinate} coordinate."
2976 ) from err
2978 # Make postage stamp plots if stamp_coordinate exists and has more than a
2979 # single point. If single_plot is True:
2980 # -- all postage stamp plots will be plotted in a single plot instead of
2981 # separate postage stamp plots.
2982 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2983 # produced per single model only
2984 if num_models == 1: 2984 ↛ 2997line 2984 didn't jump to line 2997 because the condition on line 2984 was always true
2985 if ( 2985 ↛ 2989line 2985 didn't jump to line 2989 because the condition on line 2985 was never true
2986 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2987 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2988 ):
2989 if single_plot:
2990 plotting_func = (
2991 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2992 )
2993 else:
2994 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2995 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2996 else:
2997 all_points = sorted(
2998 set(
2999 itertools.chain.from_iterable(
3000 cb.coord(sequence_coordinate).points for cb in cubes
3001 )
3002 )
3003 )
3004 all_slices = list(
3005 itertools.chain.from_iterable(
3006 cb.slices_over(sequence_coordinate) for cb in cubes
3007 )
3008 )
3009 # Matched slices (matched by seq coord point; it may happen that
3010 # evaluated models do not cover the same seq coord range, hence matching
3011 # necessary)
3012 cube_iterables = [
3013 iris.cube.CubeList(
3014 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3015 )
3016 for point in all_points
3017 ]
3019 plot_index = []
3020 nplot = np.size(cube.coord(sequence_coordinate).points)
3021 # Create a plot for each value of the sequence coordinate. Allowing for
3022 # multiple cubes in a CubeList to be plotted in the same plot for similar
3023 # sequence values. Passing a CubeList into the internal plotting function
3024 # for similar values of the sequence coordinate. cube_slice can be an
3025 # iris.cube.Cube or an iris.cube.CubeList.
3026 for cube_slice in cube_iterables:
3027 single_cube = cube_slice
3028 if isinstance(cube_slice, iris.cube.CubeList): 3028 ↛ 3029line 3028 didn't jump to line 3029 because the condition on line 3028 was never true
3029 single_cube = cube_slice[0]
3031 # Use sequence value so multiple sequences can merge.
3032 sequence_value = single_cube.coord(sequence_coordinate).points[0]
3033 plot_filename = f"{filename.rsplit('.', 1)[0]}_{sequence_value}.png"
3034 coord = single_cube.coord(sequence_coordinate)
3035 # Format the coordinate value in a unit appropriate way.
3036 title = f"{recipe_title}\n [{coord.units.title(coord.points[0])}]"
3037 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
3038 if nplot == 1 and coord.has_bounds: 3038 ↛ 3042line 3038 didn't jump to line 3042 because the condition on line 3038 was always true
3039 if np.size(coord.bounds) > 1:
3040 title = f"{recipe_title}\n [{coord.units.title(coord.bounds[0][0])} to {coord.units.title(coord.bounds[0][1])}]"
3041 # Do the actual plotting.
3042 plotting_func(
3043 cube_slice,
3044 filename=plot_filename,
3045 stamp_coordinate=stamp_coordinate,
3046 title=title,
3047 )
3048 plot_index.append(plot_filename)
3050 # Add list of plots to plot metadata.
3051 complete_plot_index = _append_to_plot_index(plot_index)
3053 # Make a page to display the plots.
3054 _make_plot_html_page(complete_plot_index)
3056 return cubes
3059def _DCT_ps(y_3d):
3060 """Calculate power spectra for regional domains.
3062 Parameters
3063 ----------
3064 y_3d: 3D array
3065 3 dimensional array to calculate spectrum for.
3066 (2D field data with 3rd dimension of time)
3068 Returns: ps_array
3069 Array of power spectra values calculated for input field (for each time)
3071 Method for regional domains:
3072 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3073 as described in Denis et al 2002 [Denis_etal_2002]_.
3075 References
3076 ----------
3077 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3078 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3079 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3080 Monthly Weather Review, Vol. 130, 1812-1828
3081 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3082 """
3083 Nt, Ny, Nx = y_3d.shape
3085 # Max coefficient
3086 Nmin = min(Nx - 1, Ny - 1)
3088 # Create alpha matrix (of wavenumbers)
3089 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3091 # Prepare output array
3092 ps_array = np.zeros((Nt, Nmin))
3094 # Loop over time to get spectrum for each time.
3095 for t in range(Nt):
3096 y_2d = y_3d[t]
3098 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3099 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3100 # cosine basis functions at different spatial frequencies.
3102 # normalise spectrum to allow comparison between models.
3103 fkk = fft.dctn(y_2d, norm="ortho")
3105 # Normalise fkk
3106 fkk = fkk / np.sqrt(Ny * Nx)
3108 # calculate variance of spectral coefficient
3109 sigma_2 = fkk**2 / Nx / Ny
3111 # Group ellipses of alphas into the same wavenumber k/Nmin
3112 for k in range(1, Nmin + 1):
3113 alpha = k / Nmin
3114 alpha_p1 = (k + 1) / Nmin
3116 # Sum up elements matching k
3117 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3118 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3120 return ps_array
3123def _create_alpha_matrix(Ny, Nx):
3124 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3126 Parameters
3127 ----------
3128 Ny, Nx:
3129 Dimensions of the 2D field for which the power spectra is calculated. Used to
3130 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3131 single-scale parameter.
3133 Returns: alpha_matrix
3134 normalisation of 2D wavenumber axes, transforming the spectral domain into
3135 an elliptic coordinate system.
3137 """
3138 # Create x_indices: each row is [1, 2, ..., Nx]
3139 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3141 # Create y_indices: each column is [1, 2, ..., Ny]
3142 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3144 # Compute alpha_matrix
3145 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3147 return alpha_matrix