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