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