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