Coverage for src / CSET / operators / plot.py: 83%
1054 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-23 15:04 +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 cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.colors as mcolors
36import matplotlib.pyplot as plt
37import numpy as np
38import scipy.fft as fft
39from cartopy.mpl.geoaxes import GeoAxes
40from iris.cube import Cube
41from markdown_it import MarkdownIt
42from mpl_toolkits.axes_grid1.inset_locator import inset_axes
44from CSET._common import (
45 combine_dicts,
46 filename_slugify,
47 get_recipe_metadata,
48 iter_maybe,
49 render_file,
50 slugify,
51)
52from CSET.operators._utils import (
53 fully_equalise_attributes,
54 get_cube_yxcoordname,
55 is_transect,
56 slice_over_maybe,
57)
58from CSET.operators.collapse import collapse
59from CSET.operators.misc import _extract_common_time_points
60from CSET.operators.regrid import regrid_onto_cube
62# Use a non-interactive plotting backend.
63mpl.use("agg")
65DEFAULT_DISCRETE_COLORS = mpl.colormaps["tab10"].colors + mpl.colormaps["Accent"].colors
67############################
68# Private helper functions #
69############################
72def _append_to_plot_index(plot_index: list) -> list:
73 """Add plots into the plot index, returning the complete plot index."""
74 with open("meta.json", "r+t", encoding="UTF-8") as fp:
75 fcntl.flock(fp, fcntl.LOCK_EX)
76 fp.seek(0)
77 meta = json.load(fp)
78 complete_plot_index = meta.get("plots", [])
79 complete_plot_index = complete_plot_index + plot_index
80 meta["plots"] = complete_plot_index
81 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
82 os.getenv("DO_CASE_AGGREGATION")
83 ):
84 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
85 fp.seek(0)
86 fp.truncate()
87 json.dump(meta, fp, indent=2)
88 return complete_plot_index
91def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:
92 """Ensure a single cube is given.
94 If a CubeList of length one is given that the contained cube is returned,
95 otherwise an error is raised.
97 Parameters
98 ----------
99 cube: Cube | CubeList
100 The cube to check.
102 Returns
103 -------
104 cube: Cube
105 The checked cube.
107 Raises
108 ------
109 TypeError
110 If the input cube is not a Cube or CubeList of a single Cube.
111 """
112 if isinstance(cube, iris.cube.Cube):
113 return cube
114 if isinstance(cube, iris.cube.CubeList):
115 if len(cube) == 1:
116 return cube[0]
117 raise TypeError("Must have a single cube", cube)
120def _make_plot_html_page(plots: list):
121 """Create a HTML page to display a plot image."""
122 # Debug check that plots actually contains some strings.
123 assert isinstance(plots[0], str)
125 # Load HTML template file.
126 operator_files = importlib.resources.files()
127 template_file = operator_files.joinpath("_plot_page_template.html")
129 # Get some metadata.
130 meta = get_recipe_metadata()
131 title = meta.get("title", "Untitled")
132 description = MarkdownIt().render(meta.get("description", "*No description.*"))
134 # Prepare template variables.
135 variables = {
136 "title": title,
137 "description": description,
138 "initial_plot": plots[0],
139 "plots": plots,
140 "title_slug": slugify(title),
141 }
143 # Render template.
144 html = render_file(template_file, **variables)
146 # Save completed HTML.
147 with open("index.html", "wt", encoding="UTF-8") as fp:
148 fp.write(html)
151@functools.cache
152def _load_colorbar_map(user_colorbar_file: str = None) -> dict:
153 """Load the colorbar definitions from a file.
155 This is a separate function to make it cacheable.
156 """
157 colorbar_file = importlib.resources.files().joinpath("_colorbar_definition.json")
158 with open(colorbar_file, "rt", encoding="UTF-8") as fp:
159 colorbar = json.load(fp)
161 logging.debug("User colour bar file: %s", user_colorbar_file)
162 override_colorbar = {}
163 if user_colorbar_file:
164 try:
165 with open(user_colorbar_file, "rt", encoding="UTF-8") as fp:
166 override_colorbar = json.load(fp)
167 except FileNotFoundError:
168 logging.warning("Colorbar file does not exist. Using default values.")
170 # Overwrite values with the user supplied colorbar definition.
171 colorbar = combine_dicts(colorbar, override_colorbar)
172 return colorbar
175def _get_model_colors_map(cubes: iris.cube.CubeList | iris.cube.Cube) -> dict:
176 """Get an appropriate colors for model lines in line plots.
178 For each model in the list of cubes colors either from user provided
179 color definition file (so-called style file) or from default colors are mapped
180 to model_name attribute.
182 Parameters
183 ----------
184 cubes: CubeList or Cube
185 Cubes with model_name attribute
187 Returns
188 -------
189 model_colors_map:
190 Dictionary mapping model_name attribute to colors
191 """
192 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
193 colorbar = _load_colorbar_map(user_colorbar_file)
194 model_names = sorted(
195 filter(
196 lambda x: x is not None,
197 (cube.attributes.get("model_name", None) for cube in iter_maybe(cubes)),
198 )
199 )
200 if not model_names:
201 return {}
202 use_user_colors = all(mname in colorbar.keys() for mname in model_names)
203 if use_user_colors: 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true
204 return {mname: colorbar[mname] for mname in model_names}
206 color_list = itertools.cycle(DEFAULT_DISCRETE_COLORS)
207 return {mname: color for mname, color in zip(model_names, color_list, strict=False)}
210def _colorbar_map_levels(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
211 """Get an appropriate colorbar for the given cube.
213 For the given variable the appropriate colorbar is looked up from a
214 combination of the built-in CSET colorbar definitions, and any user supplied
215 definitions. As well as varying on variables, these definitions may also
216 exist for specific pressure levels to account for variables with
217 significantly different ranges at different heights. The colorbars also exist
218 for masks and mask differences for considering variable presence diagnostics.
219 Specific variable ranges can be separately set in user-supplied definition
220 for x- or y-axis limits, or indicate where automated range preferred.
222 Parameters
223 ----------
224 cube: Cube
225 Cube of variable for which the colorbar information is desired.
226 axis: "x", "y", optional
227 Select the levels for just this axis of a line plot. The min and max
228 can be set by xmin/xmax or ymin/ymax respectively. For variables where
229 setting a universal range is not desirable (e.g. temperature), users
230 can set ymin/ymax values to "auto" in the colorbar definitions file.
231 Where no additional xmin/xmax or ymin/ymax values are provided, the
232 axis bounds default to use the vmin/vmax values provided.
234 Returns
235 -------
236 cmap:
237 Matplotlib colormap.
238 levels:
239 List of levels to use for plotting. For continuous plots the min and max
240 should be taken as the range.
241 norm:
242 BoundaryNorm information.
243 """
244 # Grab the colorbar file from the recipe global metadata.
245 user_colorbar_file = get_recipe_metadata().get("style_file_path", None)
246 colorbar = _load_colorbar_map(user_colorbar_file)
247 cmap = None
249 try:
250 # We assume that pressure is a scalar coordinate here.
251 pressure_level_raw = cube.coord("pressure").points[0]
252 # Ensure pressure_level is a string, as it is used as a JSON key.
253 pressure_level = str(int(pressure_level_raw))
254 except iris.exceptions.CoordinateNotFoundError:
255 pressure_level = None
257 # First try long name, then standard name, then var name. This order is used
258 # as long name is the one we correct between models, so it most likely to be
259 # consistent.
260 varnames = list(filter(None, [cube.long_name, cube.standard_name, cube.var_name]))
261 for varname in varnames:
262 # Get the colormap for this variable.
263 try:
264 var_colorbar = colorbar[varname]
265 cmap = plt.get_cmap(colorbar[varname]["cmap"], 51)
266 varname_key = varname
267 break
268 except KeyError:
269 logging.debug("Cube name %s has no colorbar definition.", varname)
271 # Get colormap if it is a mask.
272 if any("mask_for_" in name for name in varnames):
273 cmap, levels, norm = _custom_colormap_mask(cube, axis=axis)
274 return cmap, levels, norm
275 # If winds on Beaufort Scale use custom colorbar and levels
276 if any("Beaufort_Scale" in name for name in varnames):
277 cmap, levels, norm = _custom_beaufort_scale(cube, axis=axis)
278 return cmap, levels, norm
279 # If probability is plotted use custom colorbar and levels
280 if any("probability_of_" in name for name in varnames):
281 cmap, levels, norm = _custom_colormap_probability(cube, axis=axis)
282 return cmap, levels, norm
283 # If aviation colour state use custom colorbar and levels
284 if any("aviation_colour_state" in name for name in varnames):
285 cmap, levels, norm = _custom_colormap_aviation_colour_state(cube)
286 return cmap, levels, norm
288 # If no valid colormap has been defined, use defaults and return.
289 if not cmap:
290 logging.warning("No colorbar definition exists for %s.", cube.name())
291 cmap, levels, norm = mpl.colormaps["viridis"], None, None
292 return cmap, levels, norm
294 # Test if pressure-level specific settings are provided for cube.
295 if pressure_level:
296 try:
297 var_colorbar = colorbar[varname_key]["pressure_levels"][pressure_level]
298 except KeyError:
299 logging.debug(
300 "%s has no colorbar definition for pressure level %s.",
301 varname,
302 pressure_level,
303 )
305 # Check for availability of x-axis or y-axis user-specific overrides
306 # for setting level bounds for line plot types and return just levels.
307 # Line plots do not need a colormap, and just use the data range.
308 if axis:
309 if axis == "x":
310 try:
311 vmin, vmax = var_colorbar["xmin"], var_colorbar["xmax"]
312 except KeyError:
313 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
314 if axis == "y":
315 try:
316 vmin, vmax = var_colorbar["ymin"], var_colorbar["ymax"]
317 except KeyError:
318 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
319 # Check if user-specified auto-scaling for this variable
320 if vmin == "auto" or vmax == "auto":
321 levels = None
322 else:
323 levels = [vmin, vmax]
324 return None, levels, None
325 # Get and use the colorbar levels for this variable if spatial or histogram.
326 else:
327 try:
328 levels = var_colorbar["levels"]
329 # Use discrete bins when levels are specified, rather
330 # than a smooth range.
331 norm = mpl.colors.BoundaryNorm(levels, ncolors=cmap.N)
332 logging.debug("Using levels for %s colorbar.", varname)
333 logging.info("Using levels: %s", levels)
334 except KeyError:
335 # Get the range for this variable.
336 vmin, vmax = var_colorbar["min"], var_colorbar["max"]
337 logging.debug("Using min and max for %s colorbar.", varname)
338 # Calculate levels from range.
339 levels = np.linspace(vmin, vmax, 101)
340 norm = None
342 # Overwrite cmap, levels and norm for specific variables that
343 # require custom colorbar_map as these can not be defined in the
344 # JSON file.
345 cmap, levels, norm = _custom_colourmap_precipitation(cube, cmap, levels, norm)
346 cmap, levels, norm = _custom_colourmap_visibility_in_air(
347 cube, cmap, levels, norm
348 )
349 cmap, levels, norm = _custom_colormap_celsius(cube, cmap, levels, norm)
350 return cmap, levels, norm
353def _setup_spatial_map(
354 cube: iris.cube.Cube,
355 figure,
356 cmap,
357 grid_size: int | None = None,
358 subplot: int | None = None,
359):
360 """Define map projections, extent and add coastlines and borderlines for spatial plots.
362 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
363 is specified, and map extent defined based on the input data.
365 Parameters
366 ----------
367 cube: Cube
368 2 dimensional (lat and lon) Cube of the data to plot.
369 figure:
370 Matplotlib Figure object holding all plot elements.
371 cmap:
372 Matplotlib colormap.
373 grid_size: int, optional
374 Size of grid for subplots if multiple spatial subplots in figure.
375 subplot: int, optional
376 Subplot index if multiple spatial subplots in figure.
378 Returns
379 -------
380 axes:
381 Matplotlib GeoAxes definition.
382 """
383 # Identify min/max plot bounds.
384 try:
385 lat_axis, lon_axis = get_cube_yxcoordname(cube)
386 x1 = np.min(cube.coord(lon_axis).points)
387 x2 = np.max(cube.coord(lon_axis).points)
388 y1 = np.min(cube.coord(lat_axis).points)
389 y2 = np.max(cube.coord(lat_axis).points)
391 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
392 if np.abs(x2 - x1) > 180.0:
393 x1 = x1 - 180.0
394 x2 = x2 - 180.0
395 logging.debug("Adjusting plot bounds to fit global extent.")
397 # Consider map projection orientation.
398 # Adapting orientation enables plotting across international dateline.
399 # Users can adapt the default central_longitude if alternative projections views.
400 if x2 > 180.0 or x1 < -180.0:
401 central_longitude = 180.0
402 else:
403 central_longitude = 0.0
405 # Define spatial map projection.
406 coord_system = cube.coord(lat_axis).coord_system
407 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
408 # Define rotated pole map projection for rotated pole inputs.
409 projection = ccrs.RotatedPole(
410 pole_longitude=coord_system.grid_north_pole_longitude,
411 pole_latitude=coord_system.grid_north_pole_latitude,
412 central_rotated_longitude=central_longitude,
413 )
414 crs = projection
415 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 415 ↛ 417line 415 didn't jump to line 417 because the condition on line 415 was never true
416 # Define Transverse Mercator projection for TM inputs.
417 projection = ccrs.TransverseMercator(
418 central_longitude=coord_system.longitude_of_central_meridian,
419 central_latitude=coord_system.latitude_of_projection_origin,
420 false_easting=coord_system.false_easting,
421 false_northing=coord_system.false_northing,
422 scale_factor=coord_system.scale_factor_at_central_meridian,
423 )
424 crs = projection
425 else:
426 # Define regular map projection for non-rotated pole inputs.
427 # Alternatives might include e.g. for global model outputs:
428 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
429 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
430 projection = ccrs.PlateCarree(central_longitude=central_longitude)
431 crs = ccrs.PlateCarree()
433 # Define axes for plot (or subplot) with required map projection.
434 if subplot is not None:
435 axes = figure.add_subplot(
436 grid_size, grid_size, subplot, projection=projection
437 )
438 else:
439 axes = figure.add_subplot(projection=projection)
441 # Add coastlines and borderlines if cube contains x and y map coordinates.
442 if cmap.name in ["viridis", "Greys"]:
443 coastcol = "magenta"
444 else:
445 coastcol = "black"
446 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
447 axes.coastlines(resolution="10m", color=coastcol)
448 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
450 # Add gridlines.
451 if subplot is None:
452 draw_labels = True
453 else:
454 draw_labels = False
455 gl = axes.gridlines(
456 alpha=0.3,
457 draw_labels=draw_labels,
458 dms=False,
459 x_inline=False,
460 y_inline=False,
461 )
462 gl.top_labels = False
463 gl.right_labels = False
465 # If is lat/lon spatial map, fix extent to keep plot tight.
466 # Specifying crs within set_extent helps ensure only data region is shown.
467 if isinstance(coord_system, iris.coord_systems.GeogCS):
468 axes.set_extent([x1, x2, y1, y2], crs=crs)
470 except ValueError:
471 # Skip if not both x and y map coordinates.
472 axes = figure.gca()
473 pass
475 return axes
478def _get_plot_resolution() -> int:
479 """Get resolution of rasterised plots in pixels per inch."""
480 return get_recipe_metadata().get("plot_resolution", 100)
483def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
484 """Return title and filename based on start and end points or bounds."""
485 if use_bounds and seq_coord.has_bounds():
486 vals = seq_coord.bounds.flatten()
487 else:
488 vals = seq_coord.points
489 start = seq_coord.units.title(vals[0])
490 end = seq_coord.units.title(vals[-1])
492 if start == end:
493 sequence_title = f"\n [{start}]"
494 sequence_fname = f"_{filename_slugify(start)}"
495 else:
496 sequence_title = f"\n [{start} to {end}]"
497 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
499 return sequence_title, sequence_fname
502def _set_title_and_filename(
503 seq_coord: iris.coords.Coord,
504 nplot: int,
505 recipe_title: str,
506 filename: str,
507):
508 """Set plot title and filename based on cube coordinate.
510 Parameters
511 ----------
512 sequence_coordinate: iris.coords.Coord
513 Coordinate about which to make a plot sequence.
514 nplot: int
515 Number of output plots to generate - controls title/naming.
516 recipe_title: str
517 Default plot title, potentially to update.
518 filename: str
519 Input plot filename, potentially to update.
521 Returns
522 -------
523 plot_title: str
524 Output formatted plot title string, based on plotted data.
525 plot_filename: str
526 Output formatted plot filename string.
527 """
528 ndim = seq_coord.ndim
529 npoints = np.size(seq_coord.points)
530 sequence_title = ""
531 sequence_fname = ""
533 # Case 1: Multiple dimension sequence input - list number of aggregated cases
534 # (e.g. aggregation histogram plots)
535 if ndim > 1:
536 ncase = np.shape(seq_coord)[0]
537 sequence_title = f"\n [{ncase} cases]"
538 sequence_fname = f"_{ncase}cases"
540 # Case 2: Single dimension input
541 else:
542 # Single sequence point
543 if npoints == 1:
544 if nplot > 1:
545 # Default labels for sequence inputs
546 sequence_value = seq_coord.units.title(seq_coord.points[0])
547 sequence_title = f"\n [{sequence_value}]"
548 sequence_fname = f"_{filename_slugify(sequence_value)}"
549 else:
550 # Aggregated attribute available where input collapsed over aggregation
551 try:
552 ncase = seq_coord.attributes["number_reference_times"]
553 sequence_title = f"\n [{ncase} cases]"
554 sequence_fname = f"_{ncase}cases"
555 except KeyError:
556 sequence_title, sequence_fname = _get_start_end_strings(
557 seq_coord, use_bounds=seq_coord.has_bounds()
558 )
559 # Multiple sequence (e.g. time) points
560 else:
561 sequence_title, sequence_fname = _get_start_end_strings(
562 seq_coord, use_bounds=False
563 )
565 # Set plot title and filename
566 plot_title = f"{recipe_title}{sequence_title}"
568 # Set plot filename, defaulting to user input if provided.
569 if filename is None:
570 filename = slugify(recipe_title)
571 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
572 else:
573 if nplot > 1:
574 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
575 else:
576 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
578 return plot_title, plot_filename
581def _plot_and_save_spatial_plot(
582 cube: iris.cube.Cube,
583 filename: str,
584 title: str,
585 method: Literal["contourf", "pcolormesh"],
586 overlay_cube: iris.cube.Cube | None = None,
587 contour_cube: iris.cube.Cube | None = None,
588 **kwargs,
589):
590 """Plot and save a spatial plot.
592 Parameters
593 ----------
594 cube: Cube
595 2 dimensional (lat and lon) Cube of the data to plot.
596 filename: str
597 Filename of the plot to write.
598 title: str
599 Plot title.
600 method: "contourf" | "pcolormesh"
601 The plotting method to use.
602 overlay_cube: Cube, optional
603 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
604 contour_cube: Cube, optional
605 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
606 """
607 # Setup plot details, size, resolution, etc.
608 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
610 # Specify the color bar
611 cmap, levels, norm = _colorbar_map_levels(cube)
613 # If overplotting, set required colorbars
614 if overlay_cube:
615 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
616 if contour_cube:
617 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
619 # Setup plot map projection, extent and coastlines and borderlines.
620 axes = _setup_spatial_map(cube, fig, cmap)
622 # Plot the field.
623 if method == "contourf":
624 # Filled contour plot of the field.
625 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
626 elif method == "pcolormesh":
627 try:
628 vmin = min(levels)
629 vmax = max(levels)
630 except TypeError:
631 vmin, vmax = None, None
632 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
633 # if levels are defined.
634 if norm is not None:
635 vmin = None
636 vmax = None
637 logging.debug("Plotting using defined levels.")
638 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
639 else:
640 raise ValueError(f"Unknown plotting method: {method}")
642 # Overplot overlay field, if required
643 if overlay_cube:
644 try:
645 over_vmin = min(over_levels)
646 over_vmax = max(over_levels)
647 except TypeError:
648 over_vmin, over_vmax = None, None
649 if over_norm is not None: 649 ↛ 650line 649 didn't jump to line 650 because the condition on line 649 was never true
650 over_vmin = None
651 over_vmax = None
652 overlay = iplt.pcolormesh(
653 overlay_cube,
654 cmap=over_cmap,
655 norm=over_norm,
656 alpha=0.8,
657 vmin=over_vmin,
658 vmax=over_vmax,
659 )
660 # Overplot contour field, if required, with contour labelling.
661 if contour_cube:
662 contour = iplt.contour(
663 contour_cube,
664 colors="darkgray",
665 levels=cntr_levels,
666 norm=cntr_norm,
667 alpha=0.5,
668 linestyles="--",
669 linewidths=1,
670 )
671 plt.clabel(contour)
673 # Check to see if transect, and if so, adjust y axis.
674 if is_transect(cube):
675 if "pressure" in [coord.name() for coord in cube.coords()]:
676 axes.invert_yaxis()
677 axes.set_yscale("log")
678 axes.set_ylim(1100, 100)
679 # If both model_level_number and level_height exists, iplt can construct
680 # plot as a function of height above orography (NOT sea level).
681 elif {"model_level_number", "level_height"}.issubset( 681 ↛ 686line 681 didn't jump to line 686 because the condition on line 681 was always true
682 {coord.name() for coord in cube.coords()}
683 ):
684 axes.set_yscale("log")
686 axes.set_title(
687 f"{title}\n"
688 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
689 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
690 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
691 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
692 fontsize=16,
693 )
695 # Inset code
696 axins = inset_axes(
697 axes,
698 width="20%",
699 height="20%",
700 loc="upper right",
701 axes_class=GeoAxes,
702 axes_kwargs={"map_projection": ccrs.PlateCarree()},
703 )
705 # Slightly transparent to reduce plot blocking.
706 axins.patch.set_alpha(0.4)
708 axins.coastlines(resolution="50m")
709 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
711 SLat, SLon, ELat, ELon = (
712 float(coord) for coord in cube.attributes["transect_coords"].split("_")
713 )
715 # Draw line between them
716 axins.plot(
717 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
718 )
720 # Plot points (note: lon, lat order for Cartopy)
721 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
722 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
724 lon_min, lon_max = sorted([SLon, ELon])
725 lat_min, lat_max = sorted([SLat, ELat])
727 # Midpoints
728 lon_mid = (lon_min + lon_max) / 2
729 lat_mid = (lat_min + lat_max) / 2
731 # Maximum half-range
732 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
733 if half_range == 0: # points identical → provide small default 733 ↛ 737line 733 didn't jump to line 737 because the condition on line 733 was always true
734 half_range = 1
736 # Set square extent
737 axins.set_extent(
738 [
739 lon_mid - half_range,
740 lon_mid + half_range,
741 lat_mid - half_range,
742 lat_mid + half_range,
743 ],
744 crs=ccrs.PlateCarree(),
745 )
747 # Ensure square aspect
748 axins.set_aspect("equal")
750 else:
751 # Add title.
752 axes.set_title(title, fontsize=16)
754 # Adjust padding if spatial plot or transect
755 if is_transect(cube):
756 yinfopad = -0.1
757 ycbarpad = 0.1
758 else:
759 yinfopad = 0.01
760 ycbarpad = 0.042
762 # Add watermark with min/max/mean. Currently not user togglable.
763 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
764 axes.annotate(
765 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
766 xy=(1, yinfopad),
767 xycoords="axes fraction",
768 xytext=(-5, 5),
769 textcoords="offset points",
770 ha="right",
771 va="bottom",
772 size=11,
773 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
774 )
776 # Add secondary colour bar for overlay_cube field if required.
777 if overlay_cube:
778 cbarB = fig.colorbar(
779 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
780 )
781 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
782 # add ticks and tick_labels for every levels if less than 20 levels exist
783 if over_levels is not None and len(over_levels) < 20: 783 ↛ 784line 783 didn't jump to line 784 because the condition on line 783 was never true
784 cbarB.set_ticks(over_levels)
785 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
786 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
787 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
788 logging.debug("Set secondary colorbar ticks and labels.")
790 # Add main colour bar.
791 cbar = fig.colorbar(
792 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
793 )
795 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
796 # add ticks and tick_labels for every levels if less than 20 levels exist
797 if levels is not None and len(levels) < 20:
798 cbar.set_ticks(levels)
799 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
800 if "rainfall" or "snowfall" or "visibility" in cube.name(): 800 ↛ 803line 800 didn't jump to line 803 because the condition on line 800 was always true
801 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
802 # Tick labels for rainfall rates from Nimrod radar data.
803 if "rainfall rate composite" in cube.name(): 803 ↛ 804line 803 didn't jump to line 804 because the condition on line 803 was never true
804 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
805 # Tick labels for rain accumulations from Nimrod radar data.
806 if "rain accumulation" in cube.name(): 806 ↛ 807line 806 didn't jump to line 807 because the condition on line 806 was never true
807 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
808 # Tick labels for model rainfall data.
809 if "surface_microphysical" in cube.name(): 809 ↛ 811line 809 didn't jump to line 811 because the condition on line 809 was always true
810 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
811 logging.debug("Set colorbar ticks and labels.")
813 # Save plot.
814 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
815 logging.info("Saved spatial plot to %s", filename)
816 plt.close(fig)
819def _plot_and_save_postage_stamp_spatial_plot(
820 cube: iris.cube.Cube,
821 filename: str,
822 stamp_coordinate: str,
823 title: str,
824 method: Literal["contourf", "pcolormesh"],
825 overlay_cube: iris.cube.Cube | None = None,
826 contour_cube: iris.cube.Cube | None = None,
827 **kwargs,
828):
829 """Plot postage stamp spatial plots from an ensemble.
831 Parameters
832 ----------
833 cube: Cube
834 Iris cube of data to be plotted. It must have the stamp coordinate.
835 filename: str
836 Filename of the plot to write.
837 stamp_coordinate: str
838 Coordinate that becomes different plots.
839 method: "contourf" | "pcolormesh"
840 The plotting method to use.
841 overlay_cube: Cube, optional
842 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
843 contour_cube: Cube, optional
844 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
846 Raises
847 ------
848 ValueError
849 If the cube doesn't have the right dimensions.
850 """
851 # Use the smallest square grid that will fit the members.
852 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
854 fig = plt.figure(figsize=(10, 10))
856 # Specify the color bar
857 cmap, levels, norm = _colorbar_map_levels(cube)
858 # If overplotting, set required colorbars
859 if overlay_cube: 859 ↛ 860line 859 didn't jump to line 860 because the condition on line 859 was never true
860 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
861 if contour_cube: 861 ↛ 862line 861 didn't jump to line 862 because the condition on line 861 was never true
862 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
864 # Make a subplot for each member.
865 for member, subplot in zip(
866 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
867 ):
868 # Setup subplot map projection, extent and coastlines and borderlines.
869 axes = _setup_spatial_map(
870 member, fig, cmap, grid_size=grid_size, subplot=subplot
871 )
872 if method == "contourf":
873 # Filled contour plot of the field.
874 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
875 elif method == "pcolormesh":
876 if levels is not None:
877 vmin = min(levels)
878 vmax = max(levels)
879 else:
880 raise TypeError("Unknown vmin and vmax range.")
881 vmin, vmax = None, None
882 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
883 # if levels are defined.
884 if norm is not None: 884 ↛ 885line 884 didn't jump to line 885 because the condition on line 884 was never true
885 vmin = None
886 vmax = None
887 # pcolormesh plot of the field.
888 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
889 else:
890 raise ValueError(f"Unknown plotting method: {method}")
892 # Overplot overlay field, if required
893 if overlay_cube: 893 ↛ 894line 893 didn't jump to line 894 because the condition on line 893 was never true
894 try:
895 over_vmin = min(over_levels)
896 over_vmax = max(over_levels)
897 except TypeError:
898 over_vmin, over_vmax = None, None
899 if over_norm is not None:
900 over_vmin = None
901 over_vmax = None
902 iplt.pcolormesh(
903 overlay_cube[member.coord(stamp_coordinate).points[0]],
904 cmap=over_cmap,
905 norm=over_norm,
906 alpha=0.6,
907 vmin=over_vmin,
908 vmax=over_vmax,
909 )
910 # Overplot contour field, if required
911 if contour_cube: 911 ↛ 912line 911 didn't jump to line 912 because the condition on line 911 was never true
912 iplt.contour(
913 contour_cube[member.coord(stamp_coordinate).points[0]],
914 colors="darkgray",
915 levels=cntr_levels,
916 norm=cntr_norm,
917 alpha=0.6,
918 linestyles="--",
919 linewidths=1,
920 )
921 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
922 axes.set_axis_off()
924 # Put the shared colorbar in its own axes.
925 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])
926 colorbar = fig.colorbar(
927 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
928 )
929 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
931 # Overall figure title.
932 fig.suptitle(title, fontsize=16)
934 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
935 logging.info("Saved contour postage stamp plot to %s", filename)
936 plt.close(fig)
939def _plot_and_save_line_series(
940 cubes: iris.cube.CubeList,
941 coords: list[iris.coords.Coord],
942 ensemble_coord: str,
943 filename: str,
944 title: str,
945 **kwargs,
946):
947 """Plot and save a 1D line series.
949 Parameters
950 ----------
951 cubes: Cube or CubeList
952 Cube or CubeList containing the cubes to plot on the y-axis.
953 coords: list[Coord]
954 Coordinates to plot on the x-axis, one per cube.
955 ensemble_coord: str
956 Ensemble coordinate in the cube.
957 filename: str
958 Filename of the plot to write.
959 title: str
960 Plot title.
961 """
962 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
964 model_colors_map = _get_model_colors_map(cubes)
966 # Store min/max ranges.
967 y_levels = []
969 # Check match-up across sequence coords gives consistent sizes
970 _validate_cubes_coords(cubes, coords)
972 for cube, coord in zip(cubes, coords, strict=True):
973 label = None
974 color = "black"
975 if model_colors_map:
976 label = cube.attributes.get("model_name")
977 color = model_colors_map.get(label)
978 for cube_slice in cube.slices_over(ensemble_coord):
979 # Label with (control) if part of an ensemble or not otherwise.
980 if cube_slice.coord(ensemble_coord).points == [0]:
981 iplt.plot(
982 coord,
983 cube_slice,
984 color=color,
985 marker="o",
986 ls="-",
987 lw=3,
988 label=f"{label} (control)"
989 if len(cube.coord(ensemble_coord).points) > 1
990 else label,
991 )
992 # Label with (perturbed) if part of an ensemble and not the control.
993 else:
994 iplt.plot(
995 coord,
996 cube_slice,
997 color=color,
998 ls="-",
999 lw=1.5,
1000 alpha=0.75,
1001 label=f"{label} (member)",
1002 )
1004 # Calculate the global min/max if multiple cubes are given.
1005 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1006 if levels is not None: 1006 ↛ 1007line 1006 didn't jump to line 1007 because the condition on line 1006 was never true
1007 y_levels.append(min(levels))
1008 y_levels.append(max(levels))
1010 # Get the current axes.
1011 ax = plt.gca()
1013 # Add some labels and tweak the style.
1014 # check if cubes[0] works for single cube if not CubeList
1015 if coords[0].name() == "time":
1016 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1017 else:
1018 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1019 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1020 ax.set_title(title, fontsize=16)
1022 ax.ticklabel_format(axis="y", useOffset=False)
1023 ax.tick_params(axis="x", labelrotation=15)
1024 ax.tick_params(axis="both", labelsize=12)
1026 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1027 if y_levels: 1027 ↛ 1028line 1027 didn't jump to line 1028 because the condition on line 1027 was never true
1028 ax.set_ylim(min(y_levels), max(y_levels))
1029 # Add zero line.
1030 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1031 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1032 logging.debug(
1033 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1034 )
1035 else:
1036 ax.autoscale()
1038 # Add gridlines
1039 ax.grid(linestyle="--", color="grey", linewidth=1)
1040 # Ientify unique labels for legend
1041 handles = list(
1042 {
1043 label: handle
1044 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1045 }.values()
1046 )
1047 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1049 # Save plot.
1050 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1051 logging.info("Saved line plot to %s", filename)
1052 plt.close(fig)
1055def _plot_and_save_vertical_line_series(
1056 cubes: iris.cube.CubeList,
1057 coords: list[iris.coords.Coord],
1058 ensemble_coord: str,
1059 filename: str,
1060 series_coordinate: str,
1061 title: str,
1062 vmin: float,
1063 vmax: float,
1064 **kwargs,
1065):
1066 """Plot and save a 1D line series in vertical.
1068 Parameters
1069 ----------
1070 cubes: CubeList
1071 1 dimensional Cube or CubeList of the data to plot on x-axis.
1072 coord: list[Coord]
1073 Coordinates to plot on the y-axis, one per cube.
1074 ensemble_coord: str
1075 Ensemble coordinate in the cube.
1076 filename: str
1077 Filename of the plot to write.
1078 series_coordinate: str
1079 Coordinate to use as vertical axis.
1080 title: str
1081 Plot title.
1082 vmin: float
1083 Minimum value for the x-axis.
1084 vmax: float
1085 Maximum value for the x-axis.
1086 """
1087 # plot the vertical pressure axis using log scale
1088 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1090 model_colors_map = _get_model_colors_map(cubes)
1092 # Check match-up across sequence coords gives consistent sizes
1093 _validate_cubes_coords(cubes, coords)
1095 for cube, coord in zip(cubes, coords, strict=True):
1096 label = None
1097 color = "black"
1098 if model_colors_map: 1098 ↛ 1099line 1098 didn't jump to line 1099 because the condition on line 1098 was never true
1099 label = cube.attributes.get("model_name")
1100 color = model_colors_map.get(label)
1102 for cube_slice in cube.slices_over(ensemble_coord):
1103 # If ensemble data given plot control member with (control)
1104 # unless single forecast.
1105 if cube_slice.coord(ensemble_coord).points == [0]:
1106 iplt.plot(
1107 cube_slice,
1108 coord,
1109 color=color,
1110 marker="o",
1111 ls="-",
1112 lw=3,
1113 label=f"{label} (control)"
1114 if len(cube.coord(ensemble_coord).points) > 1
1115 else label,
1116 )
1117 # If ensemble data given plot perturbed members with (perturbed).
1118 else:
1119 iplt.plot(
1120 cube_slice,
1121 coord,
1122 color=color,
1123 ls="-",
1124 lw=1.5,
1125 alpha=0.75,
1126 label=f"{label} (member)",
1127 )
1129 # Get the current axis
1130 ax = plt.gca()
1132 # Special handling for pressure level data.
1133 if series_coordinate == "pressure": 1133 ↛ 1155line 1133 didn't jump to line 1155 because the condition on line 1133 was always true
1134 # Invert y-axis and set to log scale.
1135 ax.invert_yaxis()
1136 ax.set_yscale("log")
1138 # Define y-ticks and labels for pressure log axis.
1139 y_tick_labels = [
1140 "1000",
1141 "850",
1142 "700",
1143 "500",
1144 "300",
1145 "200",
1146 "100",
1147 ]
1148 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1150 # Set y-axis limits and ticks.
1151 ax.set_ylim(1100, 100)
1153 # Test if series_coordinate is model level data. The UM data uses
1154 # model_level_number and lfric uses full_levels as coordinate.
1155 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1156 # Define y-ticks and labels for vertical axis.
1157 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1158 y_tick_labels = [str(int(i)) for i in y_ticks]
1159 ax.set_ylim(min(y_ticks), max(y_ticks))
1161 ax.set_yticks(y_ticks)
1162 ax.set_yticklabels(y_tick_labels)
1164 # Set x-axis limits.
1165 ax.set_xlim(vmin, vmax)
1166 # Mark y=0 if present in plot.
1167 if vmin < 0.0 and vmax > 0.0: 1167 ↛ 1168line 1167 didn't jump to line 1168 because the condition on line 1167 was never true
1168 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1170 # Add some labels and tweak the style.
1171 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1172 ax.set_xlabel(
1173 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1174 )
1175 ax.set_title(title, fontsize=16)
1176 ax.ticklabel_format(axis="x")
1177 ax.tick_params(axis="y")
1178 ax.tick_params(axis="both", labelsize=12)
1180 # Add gridlines
1181 ax.grid(linestyle="--", color="grey", linewidth=1)
1182 # Ientify unique labels for legend
1183 handles = list(
1184 {
1185 label: handle
1186 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1187 }.values()
1188 )
1189 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1191 # Save plot.
1192 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1193 logging.info("Saved line plot to %s", filename)
1194 plt.close(fig)
1197def _plot_and_save_scatter_plot(
1198 cube_x: iris.cube.Cube | iris.cube.CubeList,
1199 cube_y: iris.cube.Cube | iris.cube.CubeList,
1200 filename: str,
1201 title: str,
1202 one_to_one: bool,
1203 model_names: list[str] = None,
1204 **kwargs,
1205):
1206 """Plot and save a 2D scatter plot.
1208 Parameters
1209 ----------
1210 cube_x: Cube | CubeList
1211 1 dimensional Cube or CubeList of the data to plot on x-axis.
1212 cube_y: Cube | CubeList
1213 1 dimensional Cube or CubeList of the data to plot on y-axis.
1214 filename: str
1215 Filename of the plot to write.
1216 title: str
1217 Plot title.
1218 one_to_one: bool
1219 Whether a 1:1 line is plotted.
1220 """
1221 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1222 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1223 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1224 # over the pairs simultaneously.
1226 # Ensure cube_x and cube_y are iterable
1227 cube_x_iterable = iter_maybe(cube_x)
1228 cube_y_iterable = iter_maybe(cube_y)
1230 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1231 iplt.scatter(cube_x_iter, cube_y_iter)
1232 if one_to_one is True:
1233 plt.plot(
1234 [
1235 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1236 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1237 ],
1238 [
1239 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1240 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1241 ],
1242 "k",
1243 linestyle="--",
1244 )
1245 ax = plt.gca()
1247 # Add some labels and tweak the style.
1248 if model_names is None:
1249 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1250 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1251 else:
1252 # Add the model names, these should be order of base (x) and other (y).
1253 ax.set_xlabel(
1254 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1255 )
1256 ax.set_ylabel(
1257 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1258 )
1259 ax.set_title(title, fontsize=16)
1260 ax.ticklabel_format(axis="y", useOffset=False)
1261 ax.tick_params(axis="x", labelrotation=15)
1262 ax.tick_params(axis="both", labelsize=12)
1263 ax.autoscale()
1265 # Save plot.
1266 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1267 logging.info("Saved scatter plot to %s", filename)
1268 plt.close(fig)
1271def _plot_and_save_vector_plot(
1272 cube_u: iris.cube.Cube,
1273 cube_v: iris.cube.Cube,
1274 filename: str,
1275 title: str,
1276 method: Literal["contourf", "pcolormesh"],
1277 **kwargs,
1278):
1279 """Plot and save a 2D vector plot.
1281 Parameters
1282 ----------
1283 cube_u: Cube
1284 2 dimensional Cube of u component of the data.
1285 cube_v: Cube
1286 2 dimensional Cube of v component of the data.
1287 filename: str
1288 Filename of the plot to write.
1289 title: str
1290 Plot title.
1291 """
1292 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1294 # Create a cube containing the magnitude of the vector field.
1295 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1296 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1298 # Specify the color bar
1299 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1301 # Setup plot map projection, extent and coastlines and borderlines.
1302 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1304 if method == "contourf":
1305 # Filled contour plot of the field.
1306 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1307 elif method == "pcolormesh":
1308 try:
1309 vmin = min(levels)
1310 vmax = max(levels)
1311 except TypeError:
1312 vmin, vmax = None, None
1313 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1314 # if levels are defined.
1315 if norm is not None:
1316 vmin = None
1317 vmax = None
1318 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1319 else:
1320 raise ValueError(f"Unknown plotting method: {method}")
1322 # Check to see if transect, and if so, adjust y axis.
1323 if is_transect(cube_vec_mag):
1324 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1325 axes.invert_yaxis()
1326 axes.set_yscale("log")
1327 axes.set_ylim(1100, 100)
1328 # If both model_level_number and level_height exists, iplt can construct
1329 # plot as a function of height above orography (NOT sea level).
1330 elif {"model_level_number", "level_height"}.issubset(
1331 {coord.name() for coord in cube_vec_mag.coords()}
1332 ):
1333 axes.set_yscale("log")
1335 axes.set_title(
1336 f"{title}\n"
1337 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1338 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1339 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1340 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1341 fontsize=16,
1342 )
1344 else:
1345 # Add title.
1346 axes.set_title(title, fontsize=16)
1348 # Add watermark with min/max/mean. Currently not user togglable.
1349 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1350 axes.annotate(
1351 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}",
1352 xy=(1, -0.05),
1353 xycoords="axes fraction",
1354 xytext=(-5, 5),
1355 textcoords="offset points",
1356 ha="right",
1357 va="bottom",
1358 size=11,
1359 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1360 )
1362 # Add colour bar.
1363 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1364 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1365 # add ticks and tick_labels for every levels if less than 20 levels exist
1366 if levels is not None and len(levels) < 20:
1367 cbar.set_ticks(levels)
1368 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1370 # 30 barbs along the longest axis of the plot, or a barb per point for data
1371 # with less than 30 points.
1372 step = max(max(cube_u.shape) // 30, 1)
1373 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1375 # Save plot.
1376 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1377 logging.info("Saved vector plot to %s", filename)
1378 plt.close(fig)
1381def _plot_and_save_histogram_series(
1382 cubes: iris.cube.Cube | iris.cube.CubeList,
1383 filename: str,
1384 title: str,
1385 vmin: float,
1386 vmax: float,
1387 **kwargs,
1388):
1389 """Plot and save a histogram series.
1391 Parameters
1392 ----------
1393 cubes: Cube or CubeList
1394 2 dimensional Cube or CubeList of the data to plot as histogram.
1395 filename: str
1396 Filename of the plot to write.
1397 title: str
1398 Plot title.
1399 vmin: float
1400 minimum for colorbar
1401 vmax: float
1402 maximum for colorbar
1403 """
1404 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1405 ax = plt.gca()
1407 model_colors_map = _get_model_colors_map(cubes)
1409 # Set default that histograms will produce probability density function
1410 # at each bin (integral over range sums to 1).
1411 density = True
1413 for cube in iter_maybe(cubes):
1414 # Easier to check title (where var name originates)
1415 # than seeing if long names exist etc.
1416 # Exception case, where distribution better fits log scales/bins.
1417 if (
1418 ("surface_microphysical" in title)
1419 or ("rain accumulation" in title)
1420 or ("Nimrod_5min" in title)
1421 ):
1422 if "amount" in title: 1422 ↛ 1424line 1422 didn't jump to line 1424 because the condition on line 1422 was never true
1423 # Compute histogram following Klingaman et al. (2017): ASoP
1424 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1425 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1426 density = False
1427 else:
1428 bins = 10.0 ** (
1429 np.arange(-10, 27, 1) / 10.0
1430 ) # Suggestion from RMED toolbox.
1431 bins = np.insert(bins, 0, 0)
1432 ax.set_yscale("log")
1433 vmin = bins[1]
1434 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1435 ax.set_xscale("log")
1436 elif "lightning" in title:
1437 bins = [0, 1, 2, 3, 4, 5]
1438 else:
1439 bins = np.linspace(vmin, vmax, 51)
1440 logging.debug(
1441 "Plotting histogram with %s bins %s - %s.",
1442 np.size(bins),
1443 np.min(bins),
1444 np.max(bins),
1445 )
1447 # Reshape cube data into a single array to allow for a single histogram.
1448 # Otherwise we plot xdim histograms stacked.
1449 cube_data_1d = (cube.data).flatten()
1451 label = None
1452 color = "black"
1453 if model_colors_map: 1453 ↛ 1454line 1453 didn't jump to line 1454 because the condition on line 1453 was never true
1454 label = cube.attributes.get("model_name")
1455 color = model_colors_map[label]
1456 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1458 # Compute area under curve.
1459 if ( 1459 ↛ 1464line 1459 didn't jump to line 1464 because the condition on line 1459 was never true
1460 ("surface_microphysical" in title and "amount" in title)
1461 or ("rain_accumulation" in title)
1462 or ("Nimrod_5min" in title)
1463 ):
1464 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1465 x = x * bin_mean / x.sum()
1466 x = x[1:]
1467 y = y[1:]
1469 ax.plot(
1470 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1471 )
1473 # Add some labels and tweak the style.
1474 ax.set_title(title, fontsize=16)
1475 ax.set_xlabel(
1476 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1477 )
1478 ax.set_ylabel("Normalised probability density", fontsize=14)
1479 if ( 1479 ↛ 1484line 1479 didn't jump to line 1484 because the condition on line 1479 was never true
1480 ("surface_microphysical" in title and "amount" in title)
1481 or ("rain accumulation" in title)
1482 or ("Nimrod_5min" in title)
1483 ):
1484 ax.set_ylabel(
1485 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1486 )
1487 ax.set_xlim(vmin, vmax)
1488 ax.tick_params(axis="both", labelsize=12)
1490 # Overlay grid-lines onto histogram plot.
1491 ax.grid(linestyle="--", color="grey", linewidth=1)
1492 if model_colors_map: 1492 ↛ 1493line 1492 didn't jump to line 1493 because the condition on line 1492 was never true
1493 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1495 # Save plot.
1496 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1497 logging.info("Saved histogram plot to %s", filename)
1498 plt.close(fig)
1501def _plot_and_save_postage_stamp_histogram_series(
1502 cube: iris.cube.Cube,
1503 filename: str,
1504 title: str,
1505 stamp_coordinate: str,
1506 vmin: float,
1507 vmax: float,
1508 **kwargs,
1509):
1510 """Plot and save postage (ensemble members) stamps for a histogram series.
1512 Parameters
1513 ----------
1514 cube: Cube
1515 2 dimensional Cube of the data to plot as histogram.
1516 filename: str
1517 Filename of the plot to write.
1518 title: str
1519 Plot title.
1520 stamp_coordinate: str
1521 Coordinate that becomes different plots.
1522 vmin: float
1523 minimum for pdf x-axis
1524 vmax: float
1525 maximum for pdf x-axis
1526 """
1527 # Use the smallest square grid that will fit the members.
1528 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1530 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1531 # Make a subplot for each member.
1532 for member, subplot in zip(
1533 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1534 ):
1535 # Implicit interface is much easier here, due to needing to have the
1536 # cartopy GeoAxes generated.
1537 plt.subplot(grid_size, grid_size, subplot)
1538 # Reshape cube data into a single array to allow for a single histogram.
1539 # Otherwise we plot xdim histograms stacked.
1540 member_data_1d = (member.data).flatten()
1541 plt.hist(member_data_1d, density=True, stacked=True)
1542 ax = plt.gca()
1543 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1544 ax.set_xlim(vmin, vmax)
1546 # Overall figure title.
1547 fig.suptitle(title, fontsize=16)
1549 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1550 logging.info("Saved histogram postage stamp plot to %s", filename)
1551 plt.close(fig)
1554def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1555 cube: iris.cube.Cube,
1556 filename: str,
1557 title: str,
1558 stamp_coordinate: str,
1559 vmin: float,
1560 vmax: float,
1561 **kwargs,
1562):
1563 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1564 ax.set_title(title, fontsize=16)
1565 ax.set_xlim(vmin, vmax)
1566 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1567 ax.set_ylabel("normalised probability density", fontsize=14)
1568 # Loop over all slices along the stamp_coordinate
1569 for member in cube.slices_over(stamp_coordinate):
1570 # Flatten the member data to 1D
1571 member_data_1d = member.data.flatten()
1572 # Plot the histogram using plt.hist
1573 plt.hist(
1574 member_data_1d,
1575 density=True,
1576 stacked=True,
1577 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1578 )
1580 # Add a legend
1581 ax.legend(fontsize=16)
1583 # Save the figure to a file
1584 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1585 logging.info("Saved histogram postage stamp plot to %s", filename)
1587 # Close the figure
1588 plt.close(fig)
1591def _plot_and_save_scattermap_plot(
1592 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1593):
1594 """Plot and save a geographical scatter plot.
1596 Parameters
1597 ----------
1598 cube: Cube
1599 1 dimensional Cube of the data points with auxiliary latitude and
1600 longitude coordinates,
1601 filename: str
1602 Filename of the plot to write.
1603 title: str
1604 Plot title.
1605 projection: str
1606 Mapping projection to be used by cartopy.
1607 """
1608 # Setup plot details, size, resolution, etc.
1609 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1610 if projection is not None:
1611 # Apart from the default, the only projection we currently support is
1612 # a stereographic projection over the North Pole.
1613 if projection == "NP_Stereo":
1614 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1615 else:
1616 raise ValueError(f"Unknown projection: {projection}")
1617 else:
1618 axes = plt.axes(projection=ccrs.PlateCarree())
1620 # Scatter plot of the field. The marker size is chosen to give
1621 # symbols that decrease in size as the number of observations
1622 # increases, although the fraction of the figure covered by
1623 # symbols increases roughly as N^(1/2), disregarding overlaps,
1624 # and has been selected for the default figure size of (10, 10).
1625 # Should this be changed, the marker size should be adjusted in
1626 # proportion to the area of the figure.
1627 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1628 klon = None
1629 klat = None
1630 for kc in range(len(cube.aux_coords)):
1631 if cube.aux_coords[kc].standard_name == "latitude":
1632 klat = kc
1633 elif cube.aux_coords[kc].standard_name == "longitude":
1634 klon = kc
1635 scatter_map = iplt.scatter(
1636 cube.aux_coords[klon],
1637 cube.aux_coords[klat],
1638 c=cube.data[:],
1639 s=mrk_size,
1640 cmap="jet",
1641 edgecolors="k",
1642 )
1644 # Add coastlines and borderlines.
1645 try:
1646 axes.coastlines(resolution="10m")
1647 axes.add_feature(cfeature.BORDERS)
1648 except AttributeError:
1649 pass
1651 # Add title.
1652 axes.set_title(title, fontsize=16)
1654 # Add colour bar.
1655 cbar = fig.colorbar(scatter_map)
1656 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1658 # Save plot.
1659 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1660 logging.info("Saved geographical scatter plot to %s", filename)
1661 plt.close(fig)
1664def _plot_and_save_power_spectrum_series(
1665 cubes: iris.cube.Cube | iris.cube.CubeList,
1666 filename: str,
1667 title: str,
1668 **kwargs,
1669):
1670 """Plot and save a power spectrum series.
1672 Parameters
1673 ----------
1674 cubes: Cube or CubeList
1675 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1676 filename: str
1677 Filename of the plot to write.
1678 title: str
1679 Plot title.
1680 """
1681 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1682 ax = plt.gca()
1684 model_colors_map = _get_model_colors_map(cubes)
1686 for cube in iter_maybe(cubes):
1687 # Calculate power spectrum
1689 # Extract time coordinate and convert to datetime
1690 time_coord = cube.coord("time")
1691 time_points = time_coord.units.num2date(time_coord.points)
1693 # Choose one time point (e.g., the first one)
1694 target_time = time_points[0]
1696 # Bind target_time inside the lambda using a default argument
1697 time_constraint = iris.Constraint(
1698 time=lambda cell, target_time=target_time: cell.point == target_time
1699 )
1701 cube = cube.extract(time_constraint)
1703 if cube.ndim == 2:
1704 cube_3d = cube.data[np.newaxis, :, :]
1705 logging.debug("Adding in new axis for a 2 dimensional cube.")
1706 elif cube.ndim == 3: 1706 ↛ 1707line 1706 didn't jump to line 1707 because the condition on line 1706 was never true
1707 cube_3d = cube.data
1708 else:
1709 raise ValueError("Cube dimensions unsuitable for power spectra code")
1710 raise ValueError(
1711 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1712 )
1714 # Calculate spectra
1715 ps_array = _DCT_ps(cube_3d)
1717 ps_cube = iris.cube.Cube(
1718 ps_array,
1719 long_name="power_spectra",
1720 )
1722 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1724 # Create a frequency/wavelength array for coordinate
1725 ps_len = ps_cube.data.shape[1]
1726 freqs = np.arange(1, ps_len + 1)
1727 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1729 # Convert datetime to numeric time using original units
1730 numeric_time = time_coord.units.date2num(time_points)
1731 # Create a new DimCoord with numeric time
1732 new_time_coord = iris.coords.DimCoord(
1733 numeric_time, standard_name="time", units=time_coord.units
1734 )
1736 # Add time and frequency coordinate to spectra cube.
1737 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1738 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1740 # Extract data from the cube
1741 frequency = ps_cube.coord("frequency").points
1742 power_spectrum = ps_cube.data
1744 label = None
1745 color = "black"
1746 if model_colors_map: 1746 ↛ 1747line 1746 didn't jump to line 1747 because the condition on line 1746 was never true
1747 label = ps_cube.attributes.get("model_name")
1748 color = model_colors_map[label]
1749 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1751 # Add some labels and tweak the style.
1752 ax.set_title(title, fontsize=16)
1753 ax.set_xlabel("Wavenumber", fontsize=14)
1754 ax.set_ylabel("Power", fontsize=14)
1755 ax.tick_params(axis="both", labelsize=12)
1757 # Set log-log scale
1758 ax.set_xscale("log")
1759 ax.set_yscale("log")
1761 # Overlay grid-lines onto power spectrum plot.
1762 ax.grid(linestyle="--", color="grey", linewidth=1)
1763 if model_colors_map: 1763 ↛ 1764line 1763 didn't jump to line 1764 because the condition on line 1763 was never true
1764 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1766 # Save plot.
1767 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1768 logging.info("Saved power spectrum plot to %s", filename)
1769 plt.close(fig)
1772def _plot_and_save_postage_stamp_power_spectrum_series(
1773 cube: iris.cube.Cube,
1774 filename: str,
1775 title: str,
1776 stamp_coordinate: str,
1777 **kwargs,
1778):
1779 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1781 Parameters
1782 ----------
1783 cube: Cube
1784 2 dimensional Cube of the data to plot as power spectrum.
1785 filename: str
1786 Filename of the plot to write.
1787 title: str
1788 Plot title.
1789 stamp_coordinate: str
1790 Coordinate that becomes different plots.
1791 """
1792 # Use the smallest square grid that will fit the members.
1793 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1795 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1796 # Make a subplot for each member.
1797 for member, subplot in zip(
1798 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1799 ):
1800 # Implicit interface is much easier here, due to needing to have the
1801 # cartopy GeoAxes generated.
1802 plt.subplot(grid_size, grid_size, subplot)
1804 frequency = member.coord("frequency").points
1806 ax = plt.gca()
1807 ax.plot(frequency, member.data)
1808 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1810 # Overall figure title.
1811 fig.suptitle(title, fontsize=16)
1813 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1814 logging.info("Saved power spectra postage stamp plot to %s", filename)
1815 plt.close(fig)
1818def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1819 cube: iris.cube.Cube,
1820 filename: str,
1821 title: str,
1822 stamp_coordinate: str,
1823 **kwargs,
1824):
1825 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1826 ax.set_title(title, fontsize=16)
1827 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1828 ax.set_ylabel("Power", fontsize=14)
1829 # Loop over all slices along the stamp_coordinate
1830 for member in cube.slices_over(stamp_coordinate):
1831 frequency = member.coord("frequency").points
1832 ax.plot(
1833 frequency,
1834 member.data,
1835 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1836 )
1838 # Add a legend
1839 ax.legend(fontsize=16)
1841 # Save the figure to a file
1842 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1843 logging.info("Saved power spectra plot to %s", filename)
1845 # Close the figure
1846 plt.close(fig)
1849def _spatial_plot(
1850 method: Literal["contourf", "pcolormesh"],
1851 cube: iris.cube.Cube,
1852 filename: str | None,
1853 sequence_coordinate: str,
1854 stamp_coordinate: str,
1855 overlay_cube: iris.cube.Cube | None = None,
1856 contour_cube: iris.cube.Cube | None = None,
1857 **kwargs,
1858):
1859 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1861 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1862 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1863 is present then postage stamp plots will be produced.
1865 If an overlay_cube and/or contour_cube are specified, multiple variables can
1866 be overplotted on the same figure.
1868 Parameters
1869 ----------
1870 method: "contourf" | "pcolormesh"
1871 The plotting method to use.
1872 cube: Cube
1873 Iris cube of the data to plot. It should have two spatial dimensions,
1874 such as lat and lon, and may also have a another two dimension to be
1875 plotted sequentially and/or as postage stamp plots.
1876 filename: str | None
1877 Name of the plot to write, used as a prefix for plot sequences. If None
1878 uses the recipe name.
1879 sequence_coordinate: str
1880 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1881 This coordinate must exist in the cube.
1882 stamp_coordinate: str
1883 Coordinate about which to plot postage stamp plots. Defaults to
1884 ``"realization"``.
1885 overlay_cube: Cube | None, optional
1886 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1887 contour_cube: Cube | None, optional
1888 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1890 Raises
1891 ------
1892 ValueError
1893 If the cube doesn't have the right dimensions.
1894 TypeError
1895 If the cube isn't a single cube.
1896 """
1897 recipe_title = get_recipe_metadata().get("title", "Untitled")
1899 # Ensure we've got a single cube.
1900 cube = _check_single_cube(cube)
1902 # Make postage stamp plots if stamp_coordinate exists and has more than a
1903 # single point.
1904 plotting_func = _plot_and_save_spatial_plot
1905 try:
1906 if cube.coord(stamp_coordinate).shape[0] > 1:
1907 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1908 except iris.exceptions.CoordinateNotFoundError:
1909 pass
1911 # Produce a geographical scatter plot if the data have a
1912 # dimension called observation or model_obs_error
1913 if any( 1913 ↛ 1917line 1913 didn't jump to line 1917 because the condition on line 1913 was never true
1914 crd.var_name == "station" or crd.var_name == "model_obs_error"
1915 for crd in cube.coords()
1916 ):
1917 plotting_func = _plot_and_save_scattermap_plot
1919 # Must have a sequence coordinate.
1920 try:
1921 cube.coord(sequence_coordinate)
1922 except iris.exceptions.CoordinateNotFoundError as err:
1923 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1925 # Create a plot for each value of the sequence coordinate.
1926 plot_index = []
1927 nplot = np.size(cube.coord(sequence_coordinate).points)
1929 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1930 # Set plot titles and filename
1931 seq_coord = cube_slice.coord(sequence_coordinate)
1932 plot_title, plot_filename = _set_title_and_filename(
1933 seq_coord, nplot, recipe_title, filename
1934 )
1936 # Extract sequence slice for overlay_cube and contour_cube if required.
1937 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1938 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1940 # Do the actual plotting.
1941 plotting_func(
1942 cube_slice,
1943 filename=plot_filename,
1944 stamp_coordinate=stamp_coordinate,
1945 title=plot_title,
1946 method=method,
1947 overlay_cube=overlay_slice,
1948 contour_cube=contour_slice,
1949 **kwargs,
1950 )
1951 plot_index.append(plot_filename)
1953 # Add list of plots to plot metadata.
1954 complete_plot_index = _append_to_plot_index(plot_index)
1956 # Make a page to display the plots.
1957 _make_plot_html_page(complete_plot_index)
1960def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1961 """Get colourmap for mask.
1963 If "mask_for_" appears anywhere in the name of a cube this function will be called
1964 regardless of the name of the variable to ensure a consistent plot.
1966 Parameters
1967 ----------
1968 cube: Cube
1969 Cube of variable for which the colorbar information is desired.
1970 axis: "x", "y", optional
1971 Select the levels for just this axis of a line plot. The min and max
1972 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1973 setting a universal range is not desirable (e.g. temperature), users
1974 can set ymin/ymax values to "auto" in the colorbar definitions file.
1975 Where no additional xmin/xmax or ymin/ymax values are provided, the
1976 axis bounds default to use the vmin/vmax values provided.
1978 Returns
1979 -------
1980 cmap:
1981 Matplotlib colormap.
1982 levels:
1983 List of levels to use for plotting. For continuous plots the min and max
1984 should be taken as the range.
1985 norm:
1986 BoundaryNorm information.
1987 """
1988 if "difference" not in cube.long_name:
1989 if axis:
1990 levels = [0, 1]
1991 # Complete settings based on levels.
1992 return None, levels, None
1993 else:
1994 # Define the levels and colors.
1995 levels = [0, 1, 2]
1996 colors = ["white", "dodgerblue"]
1997 # Create a custom color map.
1998 cmap = mcolors.ListedColormap(colors)
1999 # Normalize the levels.
2000 norm = mcolors.BoundaryNorm(levels, cmap.N)
2001 logging.debug("Colourmap for %s.", cube.long_name)
2002 return cmap, levels, norm
2003 else:
2004 if axis:
2005 levels = [-1, 1]
2006 return None, levels, None
2007 else:
2008 # Search for if mask difference, set to +/- 0.5 as values plotted <
2009 # not <=.
2010 levels = [-2, -0.5, 0.5, 2]
2011 colors = ["goldenrod", "white", "teal"]
2012 cmap = mcolors.ListedColormap(colors)
2013 norm = mcolors.BoundaryNorm(levels, cmap.N)
2014 logging.debug("Colourmap for %s.", cube.long_name)
2015 return cmap, levels, norm
2018def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2019 """Get a custom colorbar for a cube in the Beaufort Scale.
2021 Specific variable ranges can be separately set in user-supplied definition
2022 for x- or y-axis limits, or indicate where automated range preferred.
2024 Parameters
2025 ----------
2026 cube: Cube
2027 Cube of variable with Beaufort Scale in name.
2028 axis: "x", "y", optional
2029 Select the levels for just this axis of a line plot. The min and max
2030 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2031 setting a universal range is not desirable (e.g. temperature), users
2032 can set ymin/ymax values to "auto" in the colorbar definitions file.
2033 Where no additional xmin/xmax or ymin/ymax values are provided, the
2034 axis bounds default to use the vmin/vmax values provided.
2036 Returns
2037 -------
2038 cmap:
2039 Matplotlib colormap.
2040 levels:
2041 List of levels to use for plotting. For continuous plots the min and max
2042 should be taken as the range.
2043 norm:
2044 BoundaryNorm information.
2045 """
2046 if "difference" not in cube.long_name:
2047 if axis:
2048 levels = [0, 12]
2049 return None, levels, None
2050 else:
2051 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
2052 colors = [
2053 "black",
2054 (0, 0, 0.6),
2055 "blue",
2056 "cyan",
2057 "green",
2058 "yellow",
2059 (1, 0.5, 0),
2060 "red",
2061 "pink",
2062 "magenta",
2063 "purple",
2064 "maroon",
2065 "white",
2066 ]
2067 cmap = mcolors.ListedColormap(colors)
2068 norm = mcolors.BoundaryNorm(levels, cmap.N)
2069 logging.info("change colormap for Beaufort Scale colorbar.")
2070 return cmap, levels, norm
2071 else:
2072 if axis:
2073 levels = [-4, 4]
2074 return None, levels, None
2075 else:
2076 levels = [
2077 -3.5,
2078 -2.5,
2079 -1.5,
2080 -0.5,
2081 0.5,
2082 1.5,
2083 2.5,
2084 3.5,
2085 ]
2086 cmap = plt.get_cmap("bwr", 8)
2087 norm = mcolors.BoundaryNorm(levels, cmap.N)
2088 return cmap, levels, norm
2091def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
2092 """Return altered colourmap for temperature with change in units to Celsius.
2094 If "Celsius" appears anywhere in the name of a cube this function will be called.
2096 Parameters
2097 ----------
2098 cube: Cube
2099 Cube of variable for which the colorbar information is desired.
2100 cmap: Matplotlib colormap.
2101 levels: List
2102 List of levels to use for plotting. For continuous plots the min and max
2103 should be taken as the range.
2104 norm: BoundaryNorm.
2106 Returns
2107 -------
2108 cmap: Matplotlib colormap.
2109 levels: List
2110 List of levels to use for plotting. For continuous plots the min and max
2111 should be taken as the range.
2112 norm: BoundaryNorm.
2113 """
2114 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2115 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
2116 levels = np.array(levels)
2117 levels -= 273
2118 levels = levels.tolist()
2119 else:
2120 # Do nothing keep the existing colourbar attributes
2121 levels = levels
2122 cmap = cmap
2123 norm = norm
2124 return cmap, levels, norm
2127def _custom_colormap_probability(
2128 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
2129):
2130 """Get a custom colorbar for a probability cube.
2132 Specific variable ranges can be separately set in user-supplied definition
2133 for x- or y-axis limits, or indicate where automated range preferred.
2135 Parameters
2136 ----------
2137 cube: Cube
2138 Cube of variable with probability in name.
2139 axis: "x", "y", optional
2140 Select the levels for just this axis of a line plot. The min and max
2141 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2142 setting a universal range is not desirable (e.g. temperature), users
2143 can set ymin/ymax values to "auto" in the colorbar definitions file.
2144 Where no additional xmin/xmax or ymin/ymax values are provided, the
2145 axis bounds default to use the vmin/vmax values provided.
2147 Returns
2148 -------
2149 cmap:
2150 Matplotlib colormap.
2151 levels:
2152 List of levels to use for plotting. For continuous plots the min and max
2153 should be taken as the range.
2154 norm:
2155 BoundaryNorm information.
2156 """
2157 if axis:
2158 levels = [0, 1]
2159 return None, levels, None
2160 else:
2161 cmap = mcolors.ListedColormap(
2162 [
2163 "#FFFFFF",
2164 "#636363",
2165 "#e1dada",
2166 "#B5CAFF",
2167 "#8FB3FF",
2168 "#7F97FF",
2169 "#ABCF63",
2170 "#E8F59E",
2171 "#FFFA14",
2172 "#FFD121",
2173 "#FFA30A",
2174 ]
2175 )
2176 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]
2177 norm = mcolors.BoundaryNorm(levels, cmap.N)
2178 return cmap, levels, norm
2181def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
2182 """Return a custom colourmap for the current recipe."""
2183 varnames_lower = [
2184 n.lower() for n in (cube.long_name, cube.standard_name, cube.var_name) if n
2185 ]
2187 is_rainfall_var = any(
2188 key in name
2189 for name in varnames_lower
2190 for key in (
2191 "surface_microphysical",
2192 "rainfall rate composite",
2193 "nimrod5min",
2194 "nimrod_5min",
2195 "rain_accumulation",
2196 "rain accumulation",
2197 )
2198 )
2200 if is_rainfall_var:
2201 print("varnames_lower ", varnames_lower)
2202 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
2203 colors = [
2204 "w",
2205 (0, 0, 0.6),
2206 "b",
2207 "c",
2208 "g",
2209 "y",
2210 (1, 0.5, 0),
2211 "r",
2212 "pink",
2213 "m",
2214 "purple",
2215 "maroon",
2216 "gray",
2217 ]
2218 # Create a custom colormap
2219 cmap = mcolors.ListedColormap(colors)
2220 # Normalize the levels
2221 norm = mcolors.BoundaryNorm(levels, cmap.N)
2222 logging.info("Using custom rainfall colourmap.")
2224 return cmap, levels, norm
2227def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
2228 """Return custom colourmap for aviation colour state.
2230 If "aviation_colour_state" appears anywhere in the name of a cube
2231 this function will be called.
2233 Parameters
2234 ----------
2235 cube: Cube
2236 Cube of variable for which the colorbar information is desired.
2238 Returns
2239 -------
2240 cmap: Matplotlib colormap.
2241 levels: List
2242 List of levels to use for plotting. For continuous plots the min and max
2243 should be taken as the range.
2244 norm: BoundaryNorm.
2245 """
2246 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
2247 colors = [
2248 "#87ceeb",
2249 "#ffffff",
2250 "#8ced69",
2251 "#ffff00",
2252 "#ffd700",
2253 "#ffa500",
2254 "#fe3620",
2255 ]
2256 # Create a custom colormap
2257 cmap = mcolors.ListedColormap(colors)
2258 # Normalise the levels
2259 norm = mcolors.BoundaryNorm(levels, cmap.N)
2260 return cmap, levels, norm
2263def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
2264 """Return a custom colourmap for the current recipe."""
2265 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2266 if (
2267 any("visibility_in_air" in name for name in varnames)
2268 and "difference" not in cube.long_name
2269 and "mask" not in cube.long_name
2270 ):
2271 # Define the levels and colors (in km)
2272 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]
2273 norm = mcolors.BoundaryNorm(levels, cmap.N)
2274 colours = [
2275 "#8f00d6",
2276 "#d10000",
2277 "#ff9700",
2278 "#ffff00",
2279 "#00007f",
2280 "#6c9ccd",
2281 "#aae8ff",
2282 "#37a648",
2283 "#8edc64",
2284 "#c5ffc5",
2285 "#dcdcdc",
2286 "#ffffff",
2287 ]
2288 # Create a custom colormap
2289 cmap = mcolors.ListedColormap(colours)
2290 # Normalize the levels
2291 norm = mcolors.BoundaryNorm(levels, cmap.N)
2292 logging.info("change colormap for visibility_in_air variable colorbar.")
2293 else:
2294 # do nothing and keep existing colorbar attributes
2295 cmap = cmap
2296 levels = levels
2297 norm = norm
2298 return cmap, levels, norm
2301def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
2302 """Return number of models based on cube attributes."""
2303 model_names = list(
2304 filter(
2305 lambda x: x is not None,
2306 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
2307 )
2308 )
2309 if not model_names:
2310 logging.debug("Missing model names. Will assume single model.")
2311 return 1
2312 else:
2313 return len(model_names)
2316def _validate_cube_shape(
2317 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2318) -> None:
2319 """Check all cubes have a model name."""
2320 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2320 ↛ 2321line 2320 didn't jump to line 2321 because the condition on line 2320 was never true
2321 raise ValueError(
2322 f"The number of model names ({num_models}) should equal the number "
2323 f"of cubes ({len(cube)})."
2324 )
2327def _validate_cubes_coords(
2328 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2329) -> None:
2330 """Check same number of cubes as sequence coordinate for zip functions."""
2331 if len(cubes) != len(coords): 2331 ↛ 2332line 2331 didn't jump to line 2332 because the condition on line 2331 was never true
2332 raise ValueError(
2333 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2334 f"of sequence coordinates ({len(coords)})."
2335 f"Check that number of time entries in input data are consistent if "
2336 f"performing time-averaging steps prior to plotting outputs."
2337 )
2340####################
2341# Public functions #
2342####################
2345def spatial_contour_plot(
2346 cube: iris.cube.Cube,
2347 filename: str = None,
2348 sequence_coordinate: str = "time",
2349 stamp_coordinate: str = "realization",
2350 **kwargs,
2351) -> iris.cube.Cube:
2352 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2354 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2355 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2356 is present then postage stamp plots will be produced.
2358 Parameters
2359 ----------
2360 cube: Cube
2361 Iris cube of the data to plot. It should have two spatial dimensions,
2362 such as lat and lon, and may also have a another two dimension to be
2363 plotted sequentially and/or as postage stamp plots.
2364 filename: str, optional
2365 Name of the plot to write, used as a prefix for plot sequences. Defaults
2366 to the recipe name.
2367 sequence_coordinate: str, optional
2368 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2369 This coordinate must exist in the cube.
2370 stamp_coordinate: str, optional
2371 Coordinate about which to plot postage stamp plots. Defaults to
2372 ``"realization"``.
2374 Returns
2375 -------
2376 Cube
2377 The original cube (so further operations can be applied).
2379 Raises
2380 ------
2381 ValueError
2382 If the cube doesn't have the right dimensions.
2383 TypeError
2384 If the cube isn't a single cube.
2385 """
2386 _spatial_plot(
2387 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2388 )
2389 return cube
2392def spatial_pcolormesh_plot(
2393 cube: iris.cube.Cube,
2394 filename: str = None,
2395 sequence_coordinate: str = "time",
2396 stamp_coordinate: str = "realization",
2397 **kwargs,
2398) -> iris.cube.Cube:
2399 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2401 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2402 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2403 is present then postage stamp plots will be produced.
2405 This function is significantly faster than ``spatial_contour_plot``,
2406 especially at high resolutions, and should be preferred unless contiguous
2407 contour areas are important.
2409 Parameters
2410 ----------
2411 cube: Cube
2412 Iris cube of the data to plot. It should have two spatial dimensions,
2413 such as lat and lon, and may also have a another two dimension to be
2414 plotted sequentially and/or as postage stamp plots.
2415 filename: str, optional
2416 Name of the plot to write, used as a prefix for plot sequences. Defaults
2417 to the recipe name.
2418 sequence_coordinate: str, optional
2419 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2420 This coordinate must exist in the cube.
2421 stamp_coordinate: str, optional
2422 Coordinate about which to plot postage stamp plots. Defaults to
2423 ``"realization"``.
2425 Returns
2426 -------
2427 Cube
2428 The original cube (so further operations can be applied).
2430 Raises
2431 ------
2432 ValueError
2433 If the cube doesn't have the right dimensions.
2434 TypeError
2435 If the cube isn't a single cube.
2436 """
2437 _spatial_plot(
2438 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2439 )
2440 return cube
2443def spatial_multi_pcolormesh_plot(
2444 cube: iris.cube.Cube,
2445 overlay_cube: iris.cube.Cube,
2446 contour_cube: iris.cube.Cube,
2447 filename: str = None,
2448 sequence_coordinate: str = "time",
2449 stamp_coordinate: str = "realization",
2450 **kwargs,
2451) -> iris.cube.Cube:
2452 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2454 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2455 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2456 is present then postage stamp plots will be produced.
2458 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2460 If specified, contours of a contour_cube can be overplotted on top of those.
2462 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2464 This function is significantly faster than ``spatial_contour_plot``,
2465 especially at high resolutions, and should be preferred unless contiguous
2466 contour areas are important.
2468 Parameters
2469 ----------
2470 cube: Cube
2471 Iris cube of the data to plot. It should have two spatial dimensions,
2472 such as lat and lon, and may also have a another two dimension to be
2473 plotted sequentially and/or as postage stamp plots.
2474 overlay_cube: Cube
2475 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2476 such as lat and lon, and may also have a another two dimension to be
2477 plotted sequentially and/or as postage stamp plots. This is likely to be a masked cube in order not to hide the underlying basis cube.
2478 contour_cube: Cube
2479 Iris cube of the data to plot as a contour overlay on top of basis cube and overlay_cube. It should have two spatial dimensions,
2480 such as lat and lon, and may also have a another two dimension to be
2481 plotted sequentially and/or as postage stamp plots.
2482 filename: str, optional
2483 Name of the plot to write, used as a prefix for plot sequences. Defaults
2484 to the recipe name.
2485 sequence_coordinate: str, optional
2486 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2487 This coordinate must exist in the cube.
2488 stamp_coordinate: str, optional
2489 Coordinate about which to plot postage stamp plots. Defaults to
2490 ``"realization"``.
2492 Returns
2493 -------
2494 Cube
2495 The original cube (so further operations can be applied).
2497 Raises
2498 ------
2499 ValueError
2500 If the cube doesn't have the right dimensions.
2501 TypeError
2502 If the cube isn't a single cube.
2503 """
2504 _spatial_plot(
2505 "pcolormesh",
2506 cube,
2507 filename,
2508 sequence_coordinate,
2509 stamp_coordinate,
2510 overlay_cube=overlay_cube,
2511 contour_cube=contour_cube,
2512 )
2513 return cube, overlay_cube, contour_cube
2516# TODO: Expand function to handle ensemble data.
2517# line_coordinate: str, optional
2518# Coordinate about which to plot multiple lines. Defaults to
2519# ``"realization"``.
2520def plot_line_series(
2521 cube: iris.cube.Cube | iris.cube.CubeList,
2522 filename: str = None,
2523 series_coordinate: str = "time",
2524 # line_coordinate: str = "realization",
2525 **kwargs,
2526) -> iris.cube.Cube | iris.cube.CubeList:
2527 """Plot a line plot for the specified coordinate.
2529 The Cube or CubeList must be 1D.
2531 Parameters
2532 ----------
2533 iris.cube | iris.cube.CubeList
2534 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2535 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2536 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2537 filename: str, optional
2538 Name of the plot to write, used as a prefix for plot sequences. Defaults
2539 to the recipe name.
2540 series_coordinate: str, optional
2541 Coordinate about which to make a series. Defaults to ``"time"``. This
2542 coordinate must exist in the cube.
2544 Returns
2545 -------
2546 iris.cube.Cube | iris.cube.CubeList
2547 The original Cube or CubeList (so further operations can be applied).
2548 Plotted data.
2550 Raises
2551 ------
2552 ValueError
2553 If the cubes don't have the right dimensions.
2554 TypeError
2555 If the cube isn't a Cube or CubeList.
2556 """
2557 # Ensure we have a name for the plot file.
2558 recipe_title = get_recipe_metadata().get("title", "Untitled")
2560 num_models = _get_num_models(cube)
2562 _validate_cube_shape(cube, num_models)
2564 # Iterate over all cubes and extract coordinate to plot.
2565 cubes = iter_maybe(cube)
2566 coords = []
2567 for cube in cubes:
2568 try:
2569 coords.append(cube.coord(series_coordinate))
2570 except iris.exceptions.CoordinateNotFoundError as err:
2571 raise ValueError(
2572 f"Cube must have a {series_coordinate} coordinate."
2573 ) from err
2574 if cube.ndim > 2 or not cube.coords("realization"):
2575 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2577 # Format the title and filename using plotted series coordinate
2578 nplot = 1
2579 seq_coord = coords[0]
2580 plot_title, plot_filename = _set_title_and_filename(
2581 seq_coord, nplot, recipe_title, filename
2582 )
2584 # Do the actual plotting.
2585 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, plot_title)
2587 # Add list of plots to plot metadata.
2588 plot_index = _append_to_plot_index([plot_filename])
2590 # Make a page to display the plots.
2591 _make_plot_html_page(plot_index)
2593 return cube
2596def plot_vertical_line_series(
2597 cubes: iris.cube.Cube | iris.cube.CubeList,
2598 filename: str = None,
2599 series_coordinate: str = "model_level_number",
2600 sequence_coordinate: str = "time",
2601 # line_coordinate: str = "realization",
2602 **kwargs,
2603) -> iris.cube.Cube | iris.cube.CubeList:
2604 """Plot a line plot against a type of vertical coordinate.
2606 The Cube or CubeList must be 1D.
2608 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2609 then a sequence of plots will be produced.
2611 Parameters
2612 ----------
2613 iris.cube | iris.cube.CubeList
2614 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2615 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2616 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2617 filename: str, optional
2618 Name of the plot to write, used as a prefix for plot sequences. Defaults
2619 to the recipe name.
2620 series_coordinate: str, optional
2621 Coordinate to plot on the y-axis. Can be ``pressure`` or
2622 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2623 for LFRic. Defaults to ``model_level_number``.
2624 This coordinate must exist in the cube.
2625 sequence_coordinate: str, optional
2626 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2627 This coordinate must exist in the cube.
2629 Returns
2630 -------
2631 iris.cube.Cube | iris.cube.CubeList
2632 The original Cube or CubeList (so further operations can be applied).
2633 Plotted data.
2635 Raises
2636 ------
2637 ValueError
2638 If the cubes doesn't have the right dimensions.
2639 TypeError
2640 If the cube isn't a Cube or CubeList.
2641 """
2642 # Ensure we have a name for the plot file.
2643 recipe_title = get_recipe_metadata().get("title", "Untitled")
2645 cubes = iter_maybe(cubes)
2646 # Initialise empty list to hold all data from all cubes in a CubeList
2647 all_data = []
2649 # Store min/max ranges for x range.
2650 x_levels = []
2652 num_models = _get_num_models(cubes)
2654 _validate_cube_shape(cubes, num_models)
2656 # Iterate over all cubes in cube or CubeList and plot.
2657 coords = []
2658 for cube in cubes:
2659 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2660 try:
2661 coords.append(cube.coord(series_coordinate))
2662 except iris.exceptions.CoordinateNotFoundError as err:
2663 raise ValueError(
2664 f"Cube must have a {series_coordinate} coordinate."
2665 ) from err
2667 try:
2668 if cube.ndim > 1 or not cube.coords("realization"): 2668 ↛ 2676line 2668 didn't jump to line 2676 because the condition on line 2668 was always true
2669 cube.coord(sequence_coordinate)
2670 except iris.exceptions.CoordinateNotFoundError as err:
2671 raise ValueError(
2672 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2673 ) from err
2675 # Get minimum and maximum from levels information.
2676 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2677 if levels is not None: 2677 ↛ 2681line 2677 didn't jump to line 2681 because the condition on line 2677 was always true
2678 x_levels.append(min(levels))
2679 x_levels.append(max(levels))
2680 else:
2681 all_data.append(cube.data)
2683 if len(x_levels) == 0: 2683 ↛ 2685line 2683 didn't jump to line 2685 because the condition on line 2683 was never true
2684 # Combine all data into a single NumPy array
2685 combined_data = np.concatenate(all_data)
2687 # Set the lower and upper limit for the x-axis to ensure all plots have
2688 # same range. This needs to read the whole cube over the range of the
2689 # sequence and if applicable postage stamp coordinate.
2690 vmin = np.floor(combined_data.min())
2691 vmax = np.ceil(combined_data.max())
2692 else:
2693 vmin = min(x_levels)
2694 vmax = max(x_levels)
2696 # Matching the slices (matching by seq coord point; it may happen that
2697 # evaluated models do not cover the same seq coord range, hence matching
2698 # necessary)
2699 def filter_cube_iterables(cube_iterables) -> bool:
2700 return len(cube_iterables) == len(coords)
2702 cube_iterables = filter(
2703 filter_cube_iterables,
2704 (
2705 iris.cube.CubeList(
2706 s
2707 for s in itertools.chain.from_iterable(
2708 cb.slices_over(sequence_coordinate) for cb in cubes
2709 )
2710 if s.coord(sequence_coordinate).points[0] == point
2711 )
2712 for point in sorted(
2713 set(
2714 itertools.chain.from_iterable(
2715 cb.coord(sequence_coordinate).points for cb in cubes
2716 )
2717 )
2718 )
2719 ),
2720 )
2722 # Create a plot for each value of the sequence coordinate.
2723 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2724 # similar sequence values. Passing a CubeList into the internal plotting function
2725 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2726 # or an iris.cube.CubeList.
2727 plot_index = []
2728 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2729 for cubes_slice in cube_iterables:
2730 # Format the coordinate value in a unit appropriate way.
2731 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2732 plot_title, plot_filename = _set_title_and_filename(
2733 seq_coord, nplot, recipe_title, filename
2734 )
2736 # Do the actual plotting.
2737 _plot_and_save_vertical_line_series(
2738 cubes_slice,
2739 coords,
2740 "realization",
2741 plot_filename,
2742 series_coordinate,
2743 title=plot_title,
2744 vmin=vmin,
2745 vmax=vmax,
2746 )
2747 plot_index.append(plot_filename)
2749 # Add list of plots to plot metadata.
2750 complete_plot_index = _append_to_plot_index(plot_index)
2752 # Make a page to display the plots.
2753 _make_plot_html_page(complete_plot_index)
2755 return cubes
2758def qq_plot(
2759 cubes: iris.cube.CubeList,
2760 coordinates: list[str],
2761 percentiles: list[float],
2762 model_names: list[str],
2763 filename: str = None,
2764 one_to_one: bool = True,
2765 **kwargs,
2766) -> iris.cube.CubeList:
2767 """Plot a Quantile-Quantile plot between two models for common time points.
2769 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2770 collapsed within the operator over all specified coordinates such as
2771 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2772 ensemble members to ensure a 1D cube (array).
2774 Parameters
2775 ----------
2776 cubes: iris.cube.CubeList
2777 Two cubes of the same variable with different models.
2778 coordinate: list[str]
2779 The list of coordinates to collapse over. This list should be
2780 every coordinate within the cube to result in a 1D cube around
2781 the percentile coordinate.
2782 percent: list[float]
2783 A list of percentiles to appear in the plot.
2784 model_names: list[str]
2785 A list of model names to appear on the axis of the plot.
2786 filename: str, optional
2787 Filename of the plot to write.
2788 one_to_one: bool, optional
2789 If True a 1:1 line is plotted; if False it is not. Default is True.
2791 Raises
2792 ------
2793 ValueError
2794 When the cubes are not compatible.
2796 Notes
2797 -----
2798 The quantile-quantile plot is a variant on the scatter plot representing
2799 two datasets by their quantiles (percentiles) for common time points.
2800 This plot does not use a theoretical distribution to compare against, but
2801 compares percentiles of two datasets. This plot does
2802 not use all raw data points, but plots the selected percentiles (quantiles) of
2803 each variable instead for the two datasets, thereby normalising the data for a
2804 direct comparison between the selected percentiles of the two dataset distributions.
2806 Quantile-quantile plots are valuable for comparing against
2807 observations and other models. Identical percentiles between the variables
2808 will lie on the one-to-one line implying the values correspond well to each
2809 other. Where there is a deviation from the one-to-one line a range of
2810 possibilities exist depending on how and where the data is shifted (e.g.,
2811 Wilks 2011 [Wilks2011]_).
2813 For distributions above the one-to-one line the distribution is left-skewed;
2814 below is right-skewed. A distinct break implies a bimodal distribution, and
2815 closer values/values further apart at the tails imply poor representation of
2816 the extremes.
2818 References
2819 ----------
2820 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2821 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2822 """
2823 # Check cubes using same functionality as the difference operator.
2824 if len(cubes) != 2:
2825 raise ValueError("cubes should contain exactly 2 cubes.")
2826 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2827 other: Cube = cubes.extract_cube(
2828 iris.Constraint(
2829 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2830 )
2831 )
2833 # Get spatial coord names.
2834 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2835 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2837 # Ensure cubes to compare are on common differencing grid.
2838 # This is triggered if either
2839 # i) latitude and longitude shapes are not the same. Note grid points
2840 # are not compared directly as these can differ through rounding
2841 # errors.
2842 # ii) or variables are known to often sit on different grid staggering
2843 # in different models (e.g. cell center vs cell edge), as is the case
2844 # for UM and LFRic comparisons.
2845 # In future greater choice of regridding method might be applied depending
2846 # on variable type. Linear regridding can in general be appropriate for smooth
2847 # variables. Care should be taken with interpretation of differences
2848 # given this dependency on regridding.
2849 if (
2850 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2851 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2852 ) or (
2853 base.long_name
2854 in [
2855 "eastward_wind_at_10m",
2856 "northward_wind_at_10m",
2857 "northward_wind_at_cell_centres",
2858 "eastward_wind_at_cell_centres",
2859 "zonal_wind_at_pressure_levels",
2860 "meridional_wind_at_pressure_levels",
2861 "potential_vorticity_at_pressure_levels",
2862 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2863 ]
2864 ):
2865 logging.debug(
2866 "Linear regridding base cube to other grid to compute differences"
2867 )
2868 base = regrid_onto_cube(base, other, method="Linear")
2870 # Extract just common time points.
2871 base, other = _extract_common_time_points(base, other)
2873 # Equalise attributes so we can merge.
2874 fully_equalise_attributes([base, other])
2875 logging.debug("Base: %s\nOther: %s", base, other)
2877 # Collapse cubes.
2878 base = collapse(
2879 base,
2880 coordinate=coordinates,
2881 method="PERCENTILE",
2882 additional_percent=percentiles,
2883 )
2884 other = collapse(
2885 other,
2886 coordinate=coordinates,
2887 method="PERCENTILE",
2888 additional_percent=percentiles,
2889 )
2891 # Ensure we have a name for the plot file.
2892 recipe_title = get_recipe_metadata().get("title", "Untitled")
2893 title = f"{recipe_title}"
2895 if filename is None:
2896 filename = slugify(recipe_title)
2898 # Add file extension.
2899 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2901 # Do the actual plotting on a scatter plot
2902 _plot_and_save_scatter_plot(
2903 base, other, plot_filename, title, one_to_one, model_names
2904 )
2906 # Add list of plots to plot metadata.
2907 plot_index = _append_to_plot_index([plot_filename])
2909 # Make a page to display the plots.
2910 _make_plot_html_page(plot_index)
2912 return iris.cube.CubeList([base, other])
2915def scatter_plot(
2916 cube_x: iris.cube.Cube | iris.cube.CubeList,
2917 cube_y: iris.cube.Cube | iris.cube.CubeList,
2918 filename: str = None,
2919 one_to_one: bool = True,
2920 **kwargs,
2921) -> iris.cube.CubeList:
2922 """Plot a scatter plot between two variables.
2924 Both cubes must be 1D.
2926 Parameters
2927 ----------
2928 cube_x: Cube | CubeList
2929 1 dimensional Cube of the data to plot on y-axis.
2930 cube_y: Cube | CubeList
2931 1 dimensional Cube of the data to plot on x-axis.
2932 filename: str, optional
2933 Filename of the plot to write.
2934 one_to_one: bool, optional
2935 If True a 1:1 line is plotted; if False it is not. Default is True.
2937 Returns
2938 -------
2939 cubes: CubeList
2940 CubeList of the original x and y cubes for further processing.
2942 Raises
2943 ------
2944 ValueError
2945 If the cube doesn't have the right dimensions and cubes not the same
2946 size.
2947 TypeError
2948 If the cube isn't a single cube.
2950 Notes
2951 -----
2952 Scatter plots are used for determining if there is a relationship between
2953 two variables. Positive relations have a slope going from bottom left to top
2954 right; Negative relations have a slope going from top left to bottom right.
2955 """
2956 # Iterate over all cubes in cube or CubeList and plot.
2957 for cube_iter in iter_maybe(cube_x):
2958 # Check cubes are correct shape.
2959 cube_iter = _check_single_cube(cube_iter)
2960 if cube_iter.ndim > 1:
2961 raise ValueError("cube_x must be 1D.")
2963 # Iterate over all cubes in cube or CubeList and plot.
2964 for cube_iter in iter_maybe(cube_y):
2965 # Check cubes are correct shape.
2966 cube_iter = _check_single_cube(cube_iter)
2967 if cube_iter.ndim > 1:
2968 raise ValueError("cube_y must be 1D.")
2970 # Ensure we have a name for the plot file.
2971 recipe_title = get_recipe_metadata().get("title", "Untitled")
2972 title = f"{recipe_title}"
2974 if filename is None:
2975 filename = slugify(recipe_title)
2977 # Add file extension.
2978 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2980 # Do the actual plotting.
2981 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2983 # Add list of plots to plot metadata.
2984 plot_index = _append_to_plot_index([plot_filename])
2986 # Make a page to display the plots.
2987 _make_plot_html_page(plot_index)
2989 return iris.cube.CubeList([cube_x, cube_y])
2992def vector_plot(
2993 cube_u: iris.cube.Cube,
2994 cube_v: iris.cube.Cube,
2995 filename: str = None,
2996 sequence_coordinate: str = "time",
2997 **kwargs,
2998) -> iris.cube.CubeList:
2999 """Plot a vector plot based on the input u and v components."""
3000 recipe_title = get_recipe_metadata().get("title", "Untitled")
3002 # Cubes must have a matching sequence coordinate.
3003 try:
3004 # Check that the u and v cubes have the same sequence coordinate.
3005 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3005 ↛ anywhereline 3005 didn't jump anywhere: it always raised an exception.
3006 raise ValueError("Coordinates do not match.")
3007 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3008 raise ValueError(
3009 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3010 ) from err
3012 # Create a plot for each value of the sequence coordinate.
3013 plot_index = []
3014 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3015 for cube_u_slice, cube_v_slice in zip(
3016 cube_u.slices_over(sequence_coordinate),
3017 cube_v.slices_over(sequence_coordinate),
3018 strict=True,
3019 ):
3020 # Format the coordinate value in a unit appropriate way.
3021 seq_coord = cube_u_slice.coord(sequence_coordinate)
3022 plot_title, plot_filename = _set_title_and_filename(
3023 seq_coord, nplot, recipe_title, filename
3024 )
3026 # Do the actual plotting.
3027 _plot_and_save_vector_plot(
3028 cube_u_slice,
3029 cube_v_slice,
3030 filename=plot_filename,
3031 title=plot_title,
3032 method="contourf",
3033 )
3034 plot_index.append(plot_filename)
3036 # Add list of plots to plot metadata.
3037 complete_plot_index = _append_to_plot_index(plot_index)
3039 # Make a page to display the plots.
3040 _make_plot_html_page(complete_plot_index)
3042 return iris.cube.CubeList([cube_u, cube_v])
3045def plot_histogram_series(
3046 cubes: iris.cube.Cube | iris.cube.CubeList,
3047 filename: str = None,
3048 sequence_coordinate: str = "time",
3049 stamp_coordinate: str = "realization",
3050 single_plot: bool = False,
3051 **kwargs,
3052) -> iris.cube.Cube | iris.cube.CubeList:
3053 """Plot a histogram plot for each vertical level provided.
3055 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3056 is present then a sequence of plots will be produced using the time slider
3057 functionality to scroll through histograms against time. If a
3058 stamp_coordinate is present then postage stamp plots will be produced. If
3059 stamp_coordinate and single_plot is True, all postage stamp plots will be
3060 plotted in a single plot instead of separate postage stamp plots.
3062 Parameters
3063 ----------
3064 cubes: Cube | iris.cube.CubeList
3065 Iris cube or CubeList of the data to plot. It should have a single dimension other
3066 than the stamp coordinate.
3067 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3068 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3069 filename: str, optional
3070 Name of the plot to write, used as a prefix for plot sequences. Defaults
3071 to the recipe name.
3072 sequence_coordinate: str, optional
3073 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3074 This coordinate must exist in the cube and will be used for the time
3075 slider.
3076 stamp_coordinate: str, optional
3077 Coordinate about which to plot postage stamp plots. Defaults to
3078 ``"realization"``.
3079 single_plot: bool, optional
3080 If True, all postage stamp plots will be plotted in a single plot. If
3081 False, each postage stamp plot will be plotted separately. Is only valid
3082 if stamp_coordinate exists and has more than a single point.
3084 Returns
3085 -------
3086 iris.cube.Cube | iris.cube.CubeList
3087 The original Cube or CubeList (so further operations can be applied).
3088 Plotted data.
3090 Raises
3091 ------
3092 ValueError
3093 If the cube doesn't have the right dimensions.
3094 TypeError
3095 If the cube isn't a Cube or CubeList.
3096 """
3097 recipe_title = get_recipe_metadata().get("title", "Untitled")
3099 cubes = iter_maybe(cubes)
3101 # Internal plotting function.
3102 plotting_func = _plot_and_save_histogram_series
3104 num_models = _get_num_models(cubes)
3106 _validate_cube_shape(cubes, num_models)
3108 # If several histograms are plotted with time as sequence_coordinate for the
3109 # time slider option.
3110 for cube in cubes:
3111 try:
3112 cube.coord(sequence_coordinate)
3113 except iris.exceptions.CoordinateNotFoundError as err:
3114 raise ValueError(
3115 f"Cube must have a {sequence_coordinate} coordinate."
3116 ) from err
3118 # Get minimum and maximum from levels information.
3119 levels = None
3120 for cube in cubes: 3120 ↛ 3136line 3120 didn't jump to line 3136 because the loop on line 3120 didn't complete
3121 # First check if user-specified "auto" range variable.
3122 # This maintains the value of levels as None, so proceed.
3123 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3124 if levels is None:
3125 break
3126 # If levels is changed, recheck to use the vmin,vmax or
3127 # levels-based ranges for histogram plots.
3128 _, levels, _ = _colorbar_map_levels(cube)
3129 logging.debug("levels: %s", levels)
3130 if levels is not None: 3130 ↛ 3120line 3130 didn't jump to line 3120 because the condition on line 3130 was always true
3131 vmin = min(levels)
3132 vmax = max(levels)
3133 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
3134 break
3136 if levels is None:
3137 vmin = min(cb.data.min() for cb in cubes)
3138 vmax = max(cb.data.max() for cb in cubes)
3140 # Make postage stamp plots if stamp_coordinate exists and has more than a
3141 # single point. If single_plot is True:
3142 # -- all postage stamp plots will be plotted in a single plot instead of
3143 # separate postage stamp plots.
3144 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3145 # produced per single model only
3146 if num_models == 1: 3146 ↛ 3159line 3146 didn't jump to line 3159 because the condition on line 3146 was always true
3147 if ( 3147 ↛ 3151line 3147 didn't jump to line 3151 because the condition on line 3147 was never true
3148 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3149 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3150 ):
3151 if single_plot:
3152 plotting_func = (
3153 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3154 )
3155 else:
3156 plotting_func = _plot_and_save_postage_stamp_histogram_series
3157 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3158 else:
3159 all_points = sorted(
3160 set(
3161 itertools.chain.from_iterable(
3162 cb.coord(sequence_coordinate).points for cb in cubes
3163 )
3164 )
3165 )
3166 all_slices = list(
3167 itertools.chain.from_iterable(
3168 cb.slices_over(sequence_coordinate) for cb in cubes
3169 )
3170 )
3171 # Matched slices (matched by seq coord point; it may happen that
3172 # evaluated models do not cover the same seq coord range, hence matching
3173 # necessary)
3174 cube_iterables = [
3175 iris.cube.CubeList(
3176 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3177 )
3178 for point in all_points
3179 ]
3181 plot_index = []
3182 nplot = np.size(cube.coord(sequence_coordinate).points)
3183 # Create a plot for each value of the sequence coordinate. Allowing for
3184 # multiple cubes in a CubeList to be plotted in the same plot for similar
3185 # sequence values. Passing a CubeList into the internal plotting function
3186 # for similar values of the sequence coordinate. cube_slice can be an
3187 # iris.cube.Cube or an iris.cube.CubeList.
3188 for cube_slice in cube_iterables:
3189 single_cube = cube_slice
3190 if isinstance(cube_slice, iris.cube.CubeList): 3190 ↛ 3191line 3190 didn't jump to line 3191 because the condition on line 3190 was never true
3191 single_cube = cube_slice[0]
3193 # Set plot titles and filename, based on sequence coordinate
3194 seq_coord = single_cube.coord(sequence_coordinate)
3195 # Use time coordinate in title and filename if single histogram output.
3196 if sequence_coordinate == "realization" and nplot == 1: 3196 ↛ 3197line 3196 didn't jump to line 3197 because the condition on line 3196 was never true
3197 seq_coord = single_cube.coord("time")
3198 plot_title, plot_filename = _set_title_and_filename(
3199 seq_coord, nplot, recipe_title, filename
3200 )
3202 # Do the actual plotting.
3203 plotting_func(
3204 cube_slice,
3205 filename=plot_filename,
3206 stamp_coordinate=stamp_coordinate,
3207 title=plot_title,
3208 vmin=vmin,
3209 vmax=vmax,
3210 )
3211 plot_index.append(plot_filename)
3213 # Add list of plots to plot metadata.
3214 complete_plot_index = _append_to_plot_index(plot_index)
3216 # Make a page to display the plots.
3217 _make_plot_html_page(complete_plot_index)
3219 return cubes
3222def plot_power_spectrum_series(
3223 cubes: iris.cube.Cube | iris.cube.CubeList,
3224 filename: str = None,
3225 sequence_coordinate: str = "time",
3226 stamp_coordinate: str = "realization",
3227 single_plot: bool = False,
3228 **kwargs,
3229) -> iris.cube.Cube | iris.cube.CubeList:
3230 """Plot a power spectrum plot for each vertical level provided.
3232 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
3233 is present then a sequence of plots will be produced using the time slider
3234 functionality to scroll through power spectrum against time. If a
3235 stamp_coordinate is present then postage stamp plots will be produced. If
3236 stamp_coordinate and single_plot is True, all postage stamp plots will be
3237 plotted in a single plot instead of separate postage stamp plots.
3239 Parameters
3240 ----------
3241 cubes: Cube | iris.cube.CubeList
3242 Iris cube or CubeList of the data to plot. It should have a single dimension other
3243 than the stamp coordinate.
3244 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3245 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3246 filename: str, optional
3247 Name of the plot to write, used as a prefix for plot sequences. Defaults
3248 to the recipe name.
3249 sequence_coordinate: str, optional
3250 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3251 This coordinate must exist in the cube and will be used for the time
3252 slider.
3253 stamp_coordinate: str, optional
3254 Coordinate about which to plot postage stamp plots. Defaults to
3255 ``"realization"``.
3256 single_plot: bool, optional
3257 If True, all postage stamp plots will be plotted in a single plot. If
3258 False, each postage stamp plot will be plotted separately. Is only valid
3259 if stamp_coordinate exists and has more than a single point.
3261 Returns
3262 -------
3263 iris.cube.Cube | iris.cube.CubeList
3264 The original Cube or CubeList (so further operations can be applied).
3265 Plotted data.
3267 Raises
3268 ------
3269 ValueError
3270 If the cube doesn't have the right dimensions.
3271 TypeError
3272 If the cube isn't a Cube or CubeList.
3273 """
3274 recipe_title = get_recipe_metadata().get("title", "Untitled")
3276 cubes = iter_maybe(cubes)
3278 # Internal plotting function.
3279 plotting_func = _plot_and_save_power_spectrum_series
3281 num_models = _get_num_models(cubes)
3283 _validate_cube_shape(cubes, num_models)
3285 # If several power spectra are plotted with time as sequence_coordinate for the
3286 # time slider option.
3287 for cube in cubes:
3288 try:
3289 cube.coord(sequence_coordinate)
3290 except iris.exceptions.CoordinateNotFoundError as err:
3291 raise ValueError(
3292 f"Cube must have a {sequence_coordinate} coordinate."
3293 ) from err
3295 # Make postage stamp plots if stamp_coordinate exists and has more than a
3296 # single point. If single_plot is True:
3297 # -- all postage stamp plots will be plotted in a single plot instead of
3298 # separate postage stamp plots.
3299 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3300 # produced per single model only
3301 if num_models == 1: 3301 ↛ 3314line 3301 didn't jump to line 3314 because the condition on line 3301 was always true
3302 if ( 3302 ↛ 3306line 3302 didn't jump to line 3306 because the condition on line 3302 was never true
3303 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3304 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3305 ):
3306 if single_plot:
3307 plotting_func = (
3308 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3309 )
3310 else:
3311 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3312 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3313 else:
3314 all_points = sorted(
3315 set(
3316 itertools.chain.from_iterable(
3317 cb.coord(sequence_coordinate).points for cb in cubes
3318 )
3319 )
3320 )
3321 all_slices = list(
3322 itertools.chain.from_iterable(
3323 cb.slices_over(sequence_coordinate) for cb in cubes
3324 )
3325 )
3326 # Matched slices (matched by seq coord point; it may happen that
3327 # evaluated models do not cover the same seq coord range, hence matching
3328 # necessary)
3329 cube_iterables = [
3330 iris.cube.CubeList(
3331 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3332 )
3333 for point in all_points
3334 ]
3336 plot_index = []
3337 nplot = np.size(cube.coord(sequence_coordinate).points)
3338 # Create a plot for each value of the sequence coordinate. Allowing for
3339 # multiple cubes in a CubeList to be plotted in the same plot for similar
3340 # sequence values. Passing a CubeList into the internal plotting function
3341 # for similar values of the sequence coordinate. cube_slice can be an
3342 # iris.cube.Cube or an iris.cube.CubeList.
3343 for cube_slice in cube_iterables:
3344 single_cube = cube_slice
3345 if isinstance(cube_slice, iris.cube.CubeList): 3345 ↛ 3346line 3345 didn't jump to line 3346 because the condition on line 3345 was never true
3346 single_cube = cube_slice[0]
3348 # Set plot title and filenames based on sequence values
3349 seq_coord = single_cube.coord(sequence_coordinate)
3350 plot_title, plot_filename = _set_title_and_filename(
3351 seq_coord, nplot, recipe_title, filename
3352 )
3354 # Do the actual plotting.
3355 plotting_func(
3356 cube_slice,
3357 filename=plot_filename,
3358 stamp_coordinate=stamp_coordinate,
3359 title=plot_title,
3360 )
3361 plot_index.append(plot_filename)
3363 # Add list of plots to plot metadata.
3364 complete_plot_index = _append_to_plot_index(plot_index)
3366 # Make a page to display the plots.
3367 _make_plot_html_page(complete_plot_index)
3369 return cubes
3372def _DCT_ps(y_3d):
3373 """Calculate power spectra for regional domains.
3375 Parameters
3376 ----------
3377 y_3d: 3D array
3378 3 dimensional array to calculate spectrum for.
3379 (2D field data with 3rd dimension of time)
3381 Returns: ps_array
3382 Array of power spectra values calculated for input field (for each time)
3384 Method for regional domains:
3385 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3386 as described in Denis et al 2002 [Denis_etal_2002]_.
3388 References
3389 ----------
3390 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3391 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3392 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3393 Monthly Weather Review, Vol. 130, 1812-1828
3394 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3395 """
3396 Nt, Ny, Nx = y_3d.shape
3398 # Max coefficient
3399 Nmin = min(Nx - 1, Ny - 1)
3401 # Create alpha matrix (of wavenumbers)
3402 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3404 # Prepare output array
3405 ps_array = np.zeros((Nt, Nmin))
3407 # Loop over time to get spectrum for each time.
3408 for t in range(Nt):
3409 y_2d = y_3d[t]
3411 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3412 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3413 # cosine basis functions at different spatial frequencies.
3415 # normalise spectrum to allow comparison between models.
3416 fkk = fft.dctn(y_2d, norm="ortho")
3418 # Normalise fkk
3419 fkk = fkk / np.sqrt(Ny * Nx)
3421 # calculate variance of spectral coefficient
3422 sigma_2 = fkk**2 / Nx / Ny
3424 # Group ellipses of alphas into the same wavenumber k/Nmin
3425 for k in range(1, Nmin + 1):
3426 alpha = k / Nmin
3427 alpha_p1 = (k + 1) / Nmin
3429 # Sum up elements matching k
3430 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3431 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3433 return ps_array
3436def _create_alpha_matrix(Ny, Nx):
3437 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3439 Parameters
3440 ----------
3441 Ny, Nx:
3442 Dimensions of the 2D field for which the power spectra is calculated. Used to
3443 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3444 single-scale parameter.
3446 Returns: alpha_matrix
3447 normalisation of 2D wavenumber axes, transforming the spectral domain into
3448 an elliptic coordinate system.
3450 """
3451 # Create x_indices: each row is [1, 2, ..., Nx]
3452 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3454 # Create y_indices: each column is [1, 2, ..., Ny]
3455 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3457 # Compute alpha_matrix
3458 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3460 return alpha_matrix