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