Coverage for src / CSET / operators / plot.py: 83%
1056 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-28 09:58 +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: 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, optional
378 Size of grid 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, grid_size, 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 if subplot is None:
462 draw_labels = True
463 else:
464 draw_labels = False
465 gl = axes.gridlines(
466 alpha=0.3,
467 draw_labels=draw_labels,
468 dms=False,
469 x_inline=False,
470 y_inline=False,
471 )
472 gl.top_labels = False
473 gl.right_labels = False
475 # If is lat/lon spatial map, fix extent to keep plot tight.
476 # Specifying crs within set_extent helps ensure only data region is shown.
477 if isinstance(coord_system, iris.coord_systems.GeogCS):
478 axes.set_extent([x1, x2, y1, y2], crs=crs)
480 except ValueError:
481 # Skip if not both x and y map coordinates.
482 axes = figure.gca()
483 pass
485 return axes
488def _get_plot_resolution() -> int:
489 """Get resolution of rasterised plots in pixels per inch."""
490 return get_recipe_metadata().get("plot_resolution", 100)
493def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
494 """Return title and filename based on start and end points or bounds."""
495 if use_bounds and seq_coord.has_bounds():
496 vals = seq_coord.bounds.flatten()
497 else:
498 vals = seq_coord.points
499 start = seq_coord.units.title(vals[0])
500 end = seq_coord.units.title(vals[-1])
502 if start == end:
503 sequence_title = f"\n [{start}]"
504 sequence_fname = f"_{filename_slugify(start)}"
505 else:
506 sequence_title = f"\n [{start} to {end}]"
507 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
509 # Do not include time if coord set to zero.
510 if (
511 seq_coord.units == "hours since 0001-01-01 00:00:00"
512 and vals[0] == 0
513 and vals[-1] == 0
514 ):
515 sequence_title = ""
516 sequence_fname = ""
518 return sequence_title, sequence_fname
521def _set_title_and_filename(
522 seq_coord: iris.coords.Coord,
523 nplot: int,
524 recipe_title: str,
525 filename: str,
526):
527 """Set plot title and filename based on cube coordinate.
529 Parameters
530 ----------
531 sequence_coordinate: iris.coords.Coord
532 Coordinate about which to make a plot sequence.
533 nplot: int
534 Number of output plots to generate - controls title/naming.
535 recipe_title: str
536 Default plot title, potentially to update.
537 filename: str
538 Input plot filename, potentially to update.
540 Returns
541 -------
542 plot_title: str
543 Output formatted plot title string, based on plotted data.
544 plot_filename: str
545 Output formatted plot filename string.
546 """
547 ndim = seq_coord.ndim
548 npoints = np.size(seq_coord.points)
549 sequence_title = ""
550 sequence_fname = ""
552 # Case 1: Multiple dimension sequence input - list number of aggregated cases
553 # (e.g. aggregation histogram plots)
554 if ndim > 1:
555 ncase = np.shape(seq_coord)[0]
556 sequence_title = f"\n [{ncase} cases]"
557 sequence_fname = f"_{ncase}cases"
559 # Case 2: Single dimension input
560 else:
561 # Single sequence point
562 if npoints == 1:
563 if nplot > 1:
564 # Default labels for sequence inputs
565 sequence_value = seq_coord.units.title(seq_coord.points[0])
566 sequence_title = f"\n [{sequence_value}]"
567 sequence_fname = f"_{filename_slugify(sequence_value)}"
568 else:
569 # Aggregated attribute available where input collapsed over aggregation
570 try:
571 ncase = seq_coord.attributes["number_reference_times"]
572 sequence_title = f"\n [{ncase} cases]"
573 sequence_fname = f"_{ncase}cases"
574 except KeyError:
575 sequence_title, sequence_fname = _get_start_end_strings(
576 seq_coord, use_bounds=seq_coord.has_bounds()
577 )
578 # Multiple sequence (e.g. time) points
579 else:
580 sequence_title, sequence_fname = _get_start_end_strings(
581 seq_coord, use_bounds=False
582 )
584 # Set plot title and filename
585 plot_title = f"{recipe_title}{sequence_title}"
587 # Set plot filename, defaulting to user input if provided.
588 if filename is None:
589 filename = slugify(recipe_title)
590 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
591 else:
592 if nplot > 1:
593 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
594 else:
595 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
597 return plot_title, plot_filename
600def _plot_and_save_spatial_plot(
601 cube: iris.cube.Cube,
602 filename: str,
603 title: str,
604 method: Literal["contourf", "pcolormesh"],
605 overlay_cube: iris.cube.Cube | None = None,
606 contour_cube: iris.cube.Cube | None = None,
607 **kwargs,
608):
609 """Plot and save a spatial plot.
611 Parameters
612 ----------
613 cube: Cube
614 2 dimensional (lat and lon) Cube of the data to plot.
615 filename: str
616 Filename of the plot to write.
617 title: str
618 Plot title.
619 method: "contourf" | "pcolormesh"
620 The plotting method to use.
621 overlay_cube: Cube, optional
622 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
623 contour_cube: Cube, optional
624 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
625 """
626 # Setup plot details, size, resolution, etc.
627 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
629 # Specify the color bar
630 cmap, levels, norm = _colorbar_map_levels(cube)
632 # If overplotting, set required colorbars
633 if overlay_cube:
634 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
635 if contour_cube:
636 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
638 # Setup plot map projection, extent and coastlines and borderlines.
639 axes = _setup_spatial_map(cube, fig, cmap)
641 # Plot the field.
642 if method == "contourf":
643 # Filled contour plot of the field.
644 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
645 elif method == "pcolormesh":
646 try:
647 vmin = min(levels)
648 vmax = max(levels)
649 except TypeError:
650 vmin, vmax = None, None
651 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
652 # if levels are defined.
653 if norm is not None:
654 vmin = None
655 vmax = None
656 logging.debug("Plotting using defined levels.")
657 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
658 else:
659 raise ValueError(f"Unknown plotting method: {method}")
661 # Overplot overlay field, if required
662 if overlay_cube:
663 try:
664 over_vmin = min(over_levels)
665 over_vmax = max(over_levels)
666 except TypeError:
667 over_vmin, over_vmax = None, None
668 if over_norm is not None: 668 ↛ 669line 668 didn't jump to line 669 because the condition on line 668 was never true
669 over_vmin = None
670 over_vmax = None
671 overlay = iplt.pcolormesh(
672 overlay_cube,
673 cmap=over_cmap,
674 norm=over_norm,
675 alpha=0.8,
676 vmin=over_vmin,
677 vmax=over_vmax,
678 )
679 # Overplot contour field, if required, with contour labelling.
680 if contour_cube:
681 contour = iplt.contour(
682 contour_cube,
683 colors="darkgray",
684 levels=cntr_levels,
685 norm=cntr_norm,
686 alpha=0.5,
687 linestyles="--",
688 linewidths=1,
689 )
690 plt.clabel(contour)
692 # Check to see if transect, and if so, adjust y axis.
693 if is_transect(cube):
694 if "pressure" in [coord.name() for coord in cube.coords()]:
695 axes.invert_yaxis()
696 axes.set_yscale("log")
697 axes.set_ylim(1100, 100)
698 # If both model_level_number and level_height exists, iplt can construct
699 # plot as a function of height above orography (NOT sea level).
700 elif {"model_level_number", "level_height"}.issubset( 700 ↛ 705line 700 didn't jump to line 705 because the condition on line 700 was always true
701 {coord.name() for coord in cube.coords()}
702 ):
703 axes.set_yscale("log")
705 axes.set_title(
706 f"{title}\n"
707 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
708 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
709 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
710 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
711 fontsize=16,
712 )
714 # Inset code
715 axins = inset_axes(
716 axes,
717 width="20%",
718 height="20%",
719 loc="upper right",
720 axes_class=GeoAxes,
721 axes_kwargs={"map_projection": ccrs.PlateCarree()},
722 )
724 # Slightly transparent to reduce plot blocking.
725 axins.patch.set_alpha(0.4)
727 axins.coastlines(resolution="50m")
728 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
730 SLat, SLon, ELat, ELon = (
731 float(coord) for coord in cube.attributes["transect_coords"].split("_")
732 )
734 # Draw line between them
735 axins.plot(
736 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
737 )
739 # Plot points (note: lon, lat order for Cartopy)
740 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
741 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
743 lon_min, lon_max = sorted([SLon, ELon])
744 lat_min, lat_max = sorted([SLat, ELat])
746 # Midpoints
747 lon_mid = (lon_min + lon_max) / 2
748 lat_mid = (lat_min + lat_max) / 2
750 # Maximum half-range
751 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
752 if half_range == 0: # points identical → provide small default 752 ↛ 756line 752 didn't jump to line 756 because the condition on line 752 was always true
753 half_range = 1
755 # Set square extent
756 axins.set_extent(
757 [
758 lon_mid - half_range,
759 lon_mid + half_range,
760 lat_mid - half_range,
761 lat_mid + half_range,
762 ],
763 crs=ccrs.PlateCarree(),
764 )
766 # Ensure square aspect
767 axins.set_aspect("equal")
769 else:
770 # Add title.
771 axes.set_title(title, fontsize=16)
773 # Adjust padding if spatial plot or transect
774 if is_transect(cube):
775 yinfopad = -0.1
776 ycbarpad = 0.1
777 else:
778 yinfopad = 0.01
779 ycbarpad = 0.042
781 # Add watermark with min/max/mean. Currently not user togglable.
782 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
783 axes.annotate(
784 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
785 xy=(0.025, yinfopad),
786 xycoords="axes fraction",
787 xytext=(-5, 5),
788 textcoords="offset points",
789 ha="left",
790 va="bottom",
791 size=11,
792 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
793 )
795 # Add secondary colour bar for overlay_cube field if required.
796 if overlay_cube:
797 cbarB = fig.colorbar(
798 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
799 )
800 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
801 # add ticks and tick_labels for every levels if less than 20 levels exist
802 if over_levels is not None and len(over_levels) < 20: 802 ↛ 803line 802 didn't jump to line 803 because the condition on line 802 was never true
803 cbarB.set_ticks(over_levels)
804 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
805 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
806 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
807 logging.debug("Set secondary colorbar ticks and labels.")
809 # Add main colour bar.
810 cbar = fig.colorbar(
811 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
812 )
814 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
815 # add ticks and tick_labels for every levels if less than 20 levels exist
816 if levels is not None and len(levels) < 20:
817 cbar.set_ticks(levels)
818 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
819 if "rainfall" or "snowfall" or "visibility" in cube.name(): 819 ↛ 821line 819 didn't jump to line 821 because the condition on line 819 was always true
820 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
821 logging.debug("Set colorbar ticks and labels.")
823 # Save plot.
824 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
825 logging.info("Saved spatial plot to %s", filename)
826 plt.close(fig)
829def _plot_and_save_postage_stamp_spatial_plot(
830 cube: iris.cube.Cube,
831 filename: str,
832 stamp_coordinate: str,
833 title: str,
834 method: Literal["contourf", "pcolormesh"],
835 overlay_cube: iris.cube.Cube | None = None,
836 contour_cube: iris.cube.Cube | None = None,
837 **kwargs,
838):
839 """Plot postage stamp spatial plots from an ensemble.
841 Parameters
842 ----------
843 cube: Cube
844 Iris cube of data to be plotted. It must have the stamp coordinate.
845 filename: str
846 Filename of the plot to write.
847 stamp_coordinate: str
848 Coordinate that becomes different plots.
849 method: "contourf" | "pcolormesh"
850 The plotting method to use.
851 overlay_cube: Cube, optional
852 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
853 contour_cube: Cube, optional
854 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
856 Raises
857 ------
858 ValueError
859 If the cube doesn't have the right dimensions.
860 """
861 # Use the smallest square grid that will fit the members.
862 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
864 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
866 # Specify the color bar
867 cmap, levels, norm = _colorbar_map_levels(cube)
868 # If overplotting, set required colorbars
869 if overlay_cube: 869 ↛ 870line 869 didn't jump to line 870 because the condition on line 869 was never true
870 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
871 if contour_cube: 871 ↛ 872line 871 didn't jump to line 872 because the condition on line 871 was never true
872 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
874 # Make a subplot for each member.
875 for member, subplot in zip(
876 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
877 ):
878 # Setup subplot map projection, extent and coastlines and borderlines.
879 axes = _setup_spatial_map(
880 member, fig, cmap, grid_size=grid_size, subplot=subplot
881 )
882 if method == "contourf":
883 # Filled contour plot of the field.
884 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
885 elif method == "pcolormesh":
886 if levels is not None:
887 vmin = min(levels)
888 vmax = max(levels)
889 else:
890 raise TypeError("Unknown vmin and vmax range.")
891 vmin, vmax = None, None
892 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
893 # if levels are defined.
894 if norm is not None: 894 ↛ 895line 894 didn't jump to line 895 because the condition on line 894 was never true
895 vmin = None
896 vmax = None
897 # pcolormesh plot of the field.
898 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
899 else:
900 raise ValueError(f"Unknown plotting method: {method}")
902 # Overplot overlay field, if required
903 if overlay_cube: 903 ↛ 904line 903 didn't jump to line 904 because the condition on line 903 was never true
904 try:
905 over_vmin = min(over_levels)
906 over_vmax = max(over_levels)
907 except TypeError:
908 over_vmin, over_vmax = None, None
909 if over_norm is not None:
910 over_vmin = None
911 over_vmax = None
912 iplt.pcolormesh(
913 overlay_cube[member.coord(stamp_coordinate).points[0]],
914 cmap=over_cmap,
915 norm=over_norm,
916 alpha=0.6,
917 vmin=over_vmin,
918 vmax=over_vmax,
919 )
920 # Overplot contour field, if required
921 if contour_cube: 921 ↛ 922line 921 didn't jump to line 922 because the condition on line 921 was never true
922 iplt.contour(
923 contour_cube[member.coord(stamp_coordinate).points[0]],
924 colors="darkgray",
925 levels=cntr_levels,
926 norm=cntr_norm,
927 alpha=0.6,
928 linestyles="--",
929 linewidths=1,
930 )
931 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
933 # Put the shared colorbar in its own axes.
934 colorbar_axes = fig.add_axes([0.15, 0.07, 0.7, 0.03])
935 colorbar = fig.colorbar(
936 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
937 )
938 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
940 # Overall figure title.
941 fig.suptitle(title, fontsize=16)
943 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
944 logging.info("Saved contour postage stamp plot to %s", filename)
945 plt.close(fig)
948def _plot_and_save_line_series(
949 cubes: iris.cube.CubeList,
950 coords: list[iris.coords.Coord],
951 ensemble_coord: str,
952 filename: str,
953 title: str,
954 **kwargs,
955):
956 """Plot and save a 1D line series.
958 Parameters
959 ----------
960 cubes: Cube or CubeList
961 Cube or CubeList containing the cubes to plot on the y-axis.
962 coords: list[Coord]
963 Coordinates to plot on the x-axis, one per cube.
964 ensemble_coord: str
965 Ensemble coordinate in the cube.
966 filename: str
967 Filename of the plot to write.
968 title: str
969 Plot title.
970 """
971 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
973 model_colors_map = _get_model_colors_map(cubes)
975 # Store min/max ranges.
976 y_levels = []
978 # Check match-up across sequence coords gives consistent sizes
979 _validate_cubes_coords(cubes, coords)
981 for cube, coord in zip(cubes, coords, strict=True):
982 label = None
983 color = "black"
984 if model_colors_map:
985 label = cube.attributes.get("model_name")
986 color = model_colors_map.get(label)
987 for cube_slice in cube.slices_over(ensemble_coord):
988 # Label with (control) if part of an ensemble or not otherwise.
989 if cube_slice.coord(ensemble_coord).points == [0]:
990 iplt.plot(
991 coord,
992 cube_slice,
993 color=color,
994 marker="o",
995 ls="-",
996 lw=3,
997 label=f"{label} (control)"
998 if len(cube.coord(ensemble_coord).points) > 1
999 else label,
1000 )
1001 # Label with (perturbed) if part of an ensemble and not the control.
1002 else:
1003 iplt.plot(
1004 coord,
1005 cube_slice,
1006 color=color,
1007 ls="-",
1008 lw=1.5,
1009 alpha=0.75,
1010 label=f"{label} (member)",
1011 )
1013 # Calculate the global min/max if multiple cubes are given.
1014 _, levels, _ = _colorbar_map_levels(cube, axis="y")
1015 if levels is not None: 1015 ↛ 1016line 1015 didn't jump to line 1016 because the condition on line 1015 was never true
1016 y_levels.append(min(levels))
1017 y_levels.append(max(levels))
1019 # Get the current axes.
1020 ax = plt.gca()
1022 # Add some labels and tweak the style.
1023 # check if cubes[0] works for single cube if not CubeList
1024 if coords[0].name() == "time":
1025 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1026 else:
1027 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1028 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1029 ax.set_title(title, fontsize=16)
1031 ax.ticklabel_format(axis="y", useOffset=False)
1032 ax.tick_params(axis="x", labelrotation=15)
1033 ax.tick_params(axis="both", labelsize=12)
1035 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1036 if y_levels: 1036 ↛ 1037line 1036 didn't jump to line 1037 because the condition on line 1036 was never true
1037 ax.set_ylim(min(y_levels), max(y_levels))
1038 # Add zero line.
1039 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1040 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1041 logging.debug(
1042 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1043 )
1044 else:
1045 ax.autoscale()
1047 # Add gridlines
1048 ax.grid(linestyle="--", color="grey", linewidth=1)
1049 # Ientify unique labels for legend
1050 handles = list(
1051 {
1052 label: handle
1053 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1054 }.values()
1055 )
1056 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1058 # Save plot.
1059 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1060 logging.info("Saved line plot to %s", filename)
1061 plt.close(fig)
1064def _plot_and_save_vertical_line_series(
1065 cubes: iris.cube.CubeList,
1066 coords: list[iris.coords.Coord],
1067 ensemble_coord: str,
1068 filename: str,
1069 series_coordinate: str,
1070 title: str,
1071 vmin: float,
1072 vmax: float,
1073 **kwargs,
1074):
1075 """Plot and save a 1D line series in vertical.
1077 Parameters
1078 ----------
1079 cubes: CubeList
1080 1 dimensional Cube or CubeList of the data to plot on x-axis.
1081 coord: list[Coord]
1082 Coordinates to plot on the y-axis, one per cube.
1083 ensemble_coord: str
1084 Ensemble coordinate in the cube.
1085 filename: str
1086 Filename of the plot to write.
1087 series_coordinate: str
1088 Coordinate to use as vertical axis.
1089 title: str
1090 Plot title.
1091 vmin: float
1092 Minimum value for the x-axis.
1093 vmax: float
1094 Maximum value for the x-axis.
1095 """
1096 # plot the vertical pressure axis using log scale
1097 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1099 model_colors_map = _get_model_colors_map(cubes)
1101 # Check match-up across sequence coords gives consistent sizes
1102 _validate_cubes_coords(cubes, coords)
1104 for cube, coord in zip(cubes, coords, strict=True):
1105 label = None
1106 color = "black"
1107 if model_colors_map: 1107 ↛ 1108line 1107 didn't jump to line 1108 because the condition on line 1107 was never true
1108 label = cube.attributes.get("model_name")
1109 color = model_colors_map.get(label)
1111 for cube_slice in cube.slices_over(ensemble_coord):
1112 # If ensemble data given plot control member with (control)
1113 # unless single forecast.
1114 if cube_slice.coord(ensemble_coord).points == [0]:
1115 iplt.plot(
1116 cube_slice,
1117 coord,
1118 color=color,
1119 marker="o",
1120 ls="-",
1121 lw=3,
1122 label=f"{label} (control)"
1123 if len(cube.coord(ensemble_coord).points) > 1
1124 else label,
1125 )
1126 # If ensemble data given plot perturbed members with (perturbed).
1127 else:
1128 iplt.plot(
1129 cube_slice,
1130 coord,
1131 color=color,
1132 ls="-",
1133 lw=1.5,
1134 alpha=0.75,
1135 label=f"{label} (member)",
1136 )
1138 # Get the current axis
1139 ax = plt.gca()
1141 # Special handling for pressure level data.
1142 if series_coordinate == "pressure": 1142 ↛ 1164line 1142 didn't jump to line 1164 because the condition on line 1142 was always true
1143 # Invert y-axis and set to log scale.
1144 ax.invert_yaxis()
1145 ax.set_yscale("log")
1147 # Define y-ticks and labels for pressure log axis.
1148 y_tick_labels = [
1149 "1000",
1150 "850",
1151 "700",
1152 "500",
1153 "300",
1154 "200",
1155 "100",
1156 ]
1157 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1159 # Set y-axis limits and ticks.
1160 ax.set_ylim(1100, 100)
1162 # Test if series_coordinate is model level data. The UM data uses
1163 # model_level_number and lfric uses full_levels as coordinate.
1164 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1165 # Define y-ticks and labels for vertical axis.
1166 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1167 y_tick_labels = [str(int(i)) for i in y_ticks]
1168 ax.set_ylim(min(y_ticks), max(y_ticks))
1170 ax.set_yticks(y_ticks)
1171 ax.set_yticklabels(y_tick_labels)
1173 # Set x-axis limits.
1174 ax.set_xlim(vmin, vmax)
1175 # Mark y=0 if present in plot.
1176 if vmin < 0.0 and vmax > 0.0: 1176 ↛ 1177line 1176 didn't jump to line 1177 because the condition on line 1176 was never true
1177 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1179 # Add some labels and tweak the style.
1180 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1181 ax.set_xlabel(
1182 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1183 )
1184 ax.set_title(title, fontsize=16)
1185 ax.ticklabel_format(axis="x")
1186 ax.tick_params(axis="y")
1187 ax.tick_params(axis="both", labelsize=12)
1189 # Add gridlines
1190 ax.grid(linestyle="--", color="grey", linewidth=1)
1191 # Ientify unique labels for legend
1192 handles = list(
1193 {
1194 label: handle
1195 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1196 }.values()
1197 )
1198 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1200 # Save plot.
1201 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1202 logging.info("Saved line plot to %s", filename)
1203 plt.close(fig)
1206def _plot_and_save_scatter_plot(
1207 cube_x: iris.cube.Cube | iris.cube.CubeList,
1208 cube_y: iris.cube.Cube | iris.cube.CubeList,
1209 filename: str,
1210 title: str,
1211 one_to_one: bool,
1212 model_names: list[str] = None,
1213 **kwargs,
1214):
1215 """Plot and save a 2D scatter plot.
1217 Parameters
1218 ----------
1219 cube_x: Cube | CubeList
1220 1 dimensional Cube or CubeList of the data to plot on x-axis.
1221 cube_y: Cube | CubeList
1222 1 dimensional Cube or CubeList of the data to plot on y-axis.
1223 filename: str
1224 Filename of the plot to write.
1225 title: str
1226 Plot title.
1227 one_to_one: bool
1228 Whether a 1:1 line is plotted.
1229 """
1230 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1231 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1232 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1233 # over the pairs simultaneously.
1235 # Ensure cube_x and cube_y are iterable
1236 cube_x_iterable = iter_maybe(cube_x)
1237 cube_y_iterable = iter_maybe(cube_y)
1239 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1240 iplt.scatter(cube_x_iter, cube_y_iter)
1241 if one_to_one is True:
1242 plt.plot(
1243 [
1244 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1245 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1246 ],
1247 [
1248 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1249 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1250 ],
1251 "k",
1252 linestyle="--",
1253 )
1254 ax = plt.gca()
1256 # Add some labels and tweak the style.
1257 if model_names is None:
1258 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1259 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1260 else:
1261 # Add the model names, these should be order of base (x) and other (y).
1262 ax.set_xlabel(
1263 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1264 )
1265 ax.set_ylabel(
1266 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1267 )
1268 ax.set_title(title, fontsize=16)
1269 ax.ticklabel_format(axis="y", useOffset=False)
1270 ax.tick_params(axis="x", labelrotation=15)
1271 ax.tick_params(axis="both", labelsize=12)
1272 ax.autoscale()
1274 # Save plot.
1275 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1276 logging.info("Saved scatter plot to %s", filename)
1277 plt.close(fig)
1280def _plot_and_save_vector_plot(
1281 cube_u: iris.cube.Cube,
1282 cube_v: iris.cube.Cube,
1283 filename: str,
1284 title: str,
1285 method: Literal["contourf", "pcolormesh"],
1286 **kwargs,
1287):
1288 """Plot and save a 2D vector plot.
1290 Parameters
1291 ----------
1292 cube_u: Cube
1293 2 dimensional Cube of u component of the data.
1294 cube_v: Cube
1295 2 dimensional Cube of v component of the data.
1296 filename: str
1297 Filename of the plot to write.
1298 title: str
1299 Plot title.
1300 """
1301 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1303 # Create a cube containing the magnitude of the vector field.
1304 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1305 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1307 # Specify the color bar
1308 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1310 # Setup plot map projection, extent and coastlines and borderlines.
1311 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1313 if method == "contourf":
1314 # Filled contour plot of the field.
1315 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1316 elif method == "pcolormesh":
1317 try:
1318 vmin = min(levels)
1319 vmax = max(levels)
1320 except TypeError:
1321 vmin, vmax = None, None
1322 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1323 # if levels are defined.
1324 if norm is not None:
1325 vmin = None
1326 vmax = None
1327 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1328 else:
1329 raise ValueError(f"Unknown plotting method: {method}")
1331 # Check to see if transect, and if so, adjust y axis.
1332 if is_transect(cube_vec_mag):
1333 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1334 axes.invert_yaxis()
1335 axes.set_yscale("log")
1336 axes.set_ylim(1100, 100)
1337 # If both model_level_number and level_height exists, iplt can construct
1338 # plot as a function of height above orography (NOT sea level).
1339 elif {"model_level_number", "level_height"}.issubset(
1340 {coord.name() for coord in cube_vec_mag.coords()}
1341 ):
1342 axes.set_yscale("log")
1344 axes.set_title(
1345 f"{title}\n"
1346 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1347 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1348 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1349 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1350 fontsize=16,
1351 )
1353 else:
1354 # Add title.
1355 axes.set_title(title, fontsize=16)
1357 # Add watermark with min/max/mean. Currently not user togglable.
1358 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1359 axes.annotate(
1360 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}",
1361 xy=(0.05, -0.05),
1362 xycoords="axes fraction",
1363 xytext=(-5, 5),
1364 textcoords="offset points",
1365 ha="right",
1366 va="bottom",
1367 size=11,
1368 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1369 )
1371 # Add colour bar.
1372 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1373 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1374 # add ticks and tick_labels for every levels if less than 20 levels exist
1375 if levels is not None and len(levels) < 20:
1376 cbar.set_ticks(levels)
1377 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1379 # 30 barbs along the longest axis of the plot, or a barb per point for data
1380 # with less than 30 points.
1381 step = max(max(cube_u.shape) // 30, 1)
1382 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1384 # Save plot.
1385 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1386 logging.info("Saved vector plot to %s", filename)
1387 plt.close(fig)
1390def _plot_and_save_histogram_series(
1391 cubes: iris.cube.Cube | iris.cube.CubeList,
1392 filename: str,
1393 title: str,
1394 vmin: float,
1395 vmax: float,
1396 **kwargs,
1397):
1398 """Plot and save a histogram series.
1400 Parameters
1401 ----------
1402 cubes: Cube or CubeList
1403 2 dimensional Cube or CubeList of the data to plot as histogram.
1404 filename: str
1405 Filename of the plot to write.
1406 title: str
1407 Plot title.
1408 vmin: float
1409 minimum for colorbar
1410 vmax: float
1411 maximum for colorbar
1412 """
1413 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1414 ax = plt.gca()
1416 model_colors_map = _get_model_colors_map(cubes)
1418 # Set default that histograms will produce probability density function
1419 # at each bin (integral over range sums to 1).
1420 density = True
1422 for cube in iter_maybe(cubes):
1423 # Easier to check title (where var name originates)
1424 # than seeing if long names exist etc.
1425 # Exception case, where distribution better fits log scales/bins.
1426 if "surface_microphysical" in title:
1427 if "amount" in title: 1427 ↛ 1429line 1427 didn't jump to line 1429 because the condition on line 1427 was never true
1428 # Compute histogram following Klingaman et al. (2017): ASoP
1429 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1430 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1431 density = False
1432 else:
1433 bins = 10.0 ** (
1434 np.arange(-10, 27, 1) / 10.0
1435 ) # Suggestion from RMED toolbox.
1436 bins = np.insert(bins, 0, 0)
1437 ax.set_yscale("log")
1438 vmin = bins[1]
1439 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1440 ax.set_xscale("log")
1441 elif "lightning" in title:
1442 bins = [0, 1, 2, 3, 4, 5]
1443 else:
1444 bins = np.linspace(vmin, vmax, 51)
1445 logging.debug(
1446 "Plotting histogram with %s bins %s - %s.",
1447 np.size(bins),
1448 np.min(bins),
1449 np.max(bins),
1450 )
1452 # Reshape cube data into a single array to allow for a single histogram.
1453 # Otherwise we plot xdim histograms stacked.
1454 cube_data_1d = (cube.data).flatten()
1456 label = None
1457 color = "black"
1458 if model_colors_map: 1458 ↛ 1459line 1458 didn't jump to line 1459 because the condition on line 1458 was never true
1459 label = cube.attributes.get("model_name")
1460 color = model_colors_map[label]
1461 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1463 # Compute area under curve.
1464 if "surface_microphysical" in title and "amount" in title: 1464 ↛ 1465line 1464 didn't jump to line 1465 because the condition on line 1464 was never true
1465 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1466 x = x * bin_mean / x.sum()
1467 x = x[1:]
1468 y = y[1:]
1470 ax.plot(
1471 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1472 )
1474 # Add some labels and tweak the style.
1475 ax.set_title(title, fontsize=16)
1476 ax.set_xlabel(
1477 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1478 )
1479 ax.set_ylabel("Normalised probability density", fontsize=14)
1480 if "surface_microphysical" in title and "amount" in title: 1480 ↛ 1481line 1480 didn't jump to line 1481 because the condition on line 1480 was never true
1481 ax.set_ylabel(
1482 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1483 )
1484 ax.set_xlim(vmin, vmax)
1485 ax.tick_params(axis="both", labelsize=12)
1487 # Overlay grid-lines onto histogram plot.
1488 ax.grid(linestyle="--", color="grey", linewidth=1)
1489 if model_colors_map: 1489 ↛ 1490line 1489 didn't jump to line 1490 because the condition on line 1489 was never true
1490 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1492 # Save plot.
1493 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1494 logging.info("Saved histogram plot to %s", filename)
1495 plt.close(fig)
1498def _plot_and_save_postage_stamp_histogram_series(
1499 cube: iris.cube.Cube,
1500 filename: str,
1501 title: str,
1502 stamp_coordinate: str,
1503 vmin: float,
1504 vmax: float,
1505 **kwargs,
1506):
1507 """Plot and save postage (ensemble members) stamps for a histogram series.
1509 Parameters
1510 ----------
1511 cube: Cube
1512 2 dimensional Cube of the data to plot as histogram.
1513 filename: str
1514 Filename of the plot to write.
1515 title: str
1516 Plot title.
1517 stamp_coordinate: str
1518 Coordinate that becomes different plots.
1519 vmin: float
1520 minimum for pdf x-axis
1521 vmax: float
1522 maximum for pdf x-axis
1523 """
1524 # Use the smallest square grid that will fit the members.
1525 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1527 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1528 # Make a subplot for each member.
1529 for member, subplot in zip(
1530 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1531 ):
1532 # Implicit interface is much easier here, due to needing to have the
1533 # cartopy GeoAxes generated.
1534 plt.subplot(grid_size, grid_size, subplot)
1535 # Reshape cube data into a single array to allow for a single histogram.
1536 # Otherwise we plot xdim histograms stacked.
1537 member_data_1d = (member.data).flatten()
1538 plt.hist(member_data_1d, density=True, stacked=True)
1539 axes = plt.gca()
1540 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1541 axes.set_xlim(vmin, vmax)
1543 # Overall figure title.
1544 fig.suptitle(title, fontsize=16)
1546 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1547 logging.info("Saved histogram postage stamp plot to %s", filename)
1548 plt.close(fig)
1551def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1552 cube: iris.cube.Cube,
1553 filename: str,
1554 title: str,
1555 stamp_coordinate: str,
1556 vmin: float,
1557 vmax: float,
1558 **kwargs,
1559):
1560 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1561 ax.set_title(title, fontsize=16)
1562 ax.set_xlim(vmin, vmax)
1563 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1564 ax.set_ylabel("normalised probability density", fontsize=14)
1565 # Loop over all slices along the stamp_coordinate
1566 for member in cube.slices_over(stamp_coordinate):
1567 # Flatten the member data to 1D
1568 member_data_1d = member.data.flatten()
1569 # Plot the histogram using plt.hist
1570 plt.hist(
1571 member_data_1d,
1572 density=True,
1573 stacked=True,
1574 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1575 )
1577 # Add a legend
1578 ax.legend(fontsize=16)
1580 # Save the figure to a file
1581 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1582 logging.info("Saved histogram postage stamp plot to %s", filename)
1584 # Close the figure
1585 plt.close(fig)
1588def _plot_and_save_scattermap_plot(
1589 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1590):
1591 """Plot and save a geographical scatter plot.
1593 Parameters
1594 ----------
1595 cube: Cube
1596 1 dimensional Cube of the data points with auxiliary latitude and
1597 longitude coordinates,
1598 filename: str
1599 Filename of the plot to write.
1600 title: str
1601 Plot title.
1602 projection: str
1603 Mapping projection to be used by cartopy.
1604 """
1605 # Setup plot details, size, resolution, etc.
1606 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1607 if projection is not None:
1608 # Apart from the default, the only projection we currently support is
1609 # a stereographic projection over the North Pole.
1610 if projection == "NP_Stereo":
1611 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1612 else:
1613 raise ValueError(f"Unknown projection: {projection}")
1614 else:
1615 axes = plt.axes(projection=ccrs.PlateCarree())
1617 # Scatter plot of the field. The marker size is chosen to give
1618 # symbols that decrease in size as the number of observations
1619 # increases, although the fraction of the figure covered by
1620 # symbols increases roughly as N^(1/2), disregarding overlaps,
1621 # and has been selected for the default figure size of (10, 10).
1622 # Should this be changed, the marker size should be adjusted in
1623 # proportion to the area of the figure.
1624 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1625 klon = None
1626 klat = None
1627 for kc in range(len(cube.aux_coords)):
1628 if cube.aux_coords[kc].standard_name == "latitude":
1629 klat = kc
1630 elif cube.aux_coords[kc].standard_name == "longitude":
1631 klon = kc
1632 scatter_map = iplt.scatter(
1633 cube.aux_coords[klon],
1634 cube.aux_coords[klat],
1635 c=cube.data[:],
1636 s=mrk_size,
1637 cmap="jet",
1638 edgecolors="k",
1639 )
1641 # Add coastlines and borderlines.
1642 try:
1643 axes.coastlines(resolution="10m")
1644 axes.add_feature(cfeature.BORDERS)
1645 except AttributeError:
1646 pass
1648 # Add title.
1649 axes.set_title(title, fontsize=16)
1651 # Add colour bar.
1652 cbar = fig.colorbar(scatter_map)
1653 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1655 # Save plot.
1656 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1657 logging.info("Saved geographical scatter plot to %s", filename)
1658 plt.close(fig)
1661def _plot_and_save_power_spectrum_series(
1662 cubes: iris.cube.Cube | iris.cube.CubeList,
1663 filename: str,
1664 title: str,
1665 **kwargs,
1666):
1667 """Plot and save a power spectrum series.
1669 Parameters
1670 ----------
1671 cubes: Cube or CubeList
1672 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1673 filename: str
1674 Filename of the plot to write.
1675 title: str
1676 Plot title.
1677 """
1678 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1679 ax = plt.gca()
1681 model_colors_map = _get_model_colors_map(cubes)
1683 for cube in iter_maybe(cubes):
1684 # Calculate power spectrum
1686 # Extract time coordinate and convert to datetime
1687 time_coord = cube.coord("time")
1688 time_points = time_coord.units.num2date(time_coord.points)
1690 # Choose one time point (e.g., the first one)
1691 target_time = time_points[0]
1693 # Bind target_time inside the lambda using a default argument
1694 time_constraint = iris.Constraint(
1695 time=lambda cell, target_time=target_time: cell.point == target_time
1696 )
1698 cube = cube.extract(time_constraint)
1700 if cube.ndim == 2:
1701 cube_3d = cube.data[np.newaxis, :, :]
1702 logging.debug("Adding in new axis for a 2 dimensional cube.")
1703 elif cube.ndim == 3: 1703 ↛ 1704line 1703 didn't jump to line 1704 because the condition on line 1703 was never true
1704 cube_3d = cube.data
1705 else:
1706 raise ValueError("Cube dimensions unsuitable for power spectra code")
1707 raise ValueError(
1708 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1709 )
1711 # Calculate spectra
1712 ps_array = _DCT_ps(cube_3d)
1714 ps_cube = iris.cube.Cube(
1715 ps_array,
1716 long_name="power_spectra",
1717 )
1719 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1721 # Create a frequency/wavelength array for coordinate
1722 ps_len = ps_cube.data.shape[1]
1723 freqs = np.arange(1, ps_len + 1)
1724 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1726 # Convert datetime to numeric time using original units
1727 numeric_time = time_coord.units.date2num(time_points)
1728 # Create a new DimCoord with numeric time
1729 new_time_coord = iris.coords.DimCoord(
1730 numeric_time, standard_name="time", units=time_coord.units
1731 )
1733 # Add time and frequency coordinate to spectra cube.
1734 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1735 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1737 # Extract data from the cube
1738 frequency = ps_cube.coord("frequency").points
1739 power_spectrum = ps_cube.data
1741 label = None
1742 color = "black"
1743 if model_colors_map: 1743 ↛ 1744line 1743 didn't jump to line 1744 because the condition on line 1743 was never true
1744 label = ps_cube.attributes.get("model_name")
1745 color = model_colors_map[label]
1746 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1748 # Add some labels and tweak the style.
1749 ax.set_title(title, fontsize=16)
1750 ax.set_xlabel("Wavenumber", fontsize=14)
1751 ax.set_ylabel("Power", fontsize=14)
1752 ax.tick_params(axis="both", labelsize=12)
1754 # Set log-log scale
1755 ax.set_xscale("log")
1756 ax.set_yscale("log")
1758 # Overlay grid-lines onto power spectrum plot.
1759 ax.grid(linestyle="--", color="grey", linewidth=1)
1760 if model_colors_map: 1760 ↛ 1761line 1760 didn't jump to line 1761 because the condition on line 1760 was never true
1761 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1763 # Save plot.
1764 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1765 logging.info("Saved power spectrum plot to %s", filename)
1766 plt.close(fig)
1769def _plot_and_save_postage_stamp_power_spectrum_series(
1770 cube: iris.cube.Cube,
1771 filename: str,
1772 title: str,
1773 stamp_coordinate: str,
1774 **kwargs,
1775):
1776 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1778 Parameters
1779 ----------
1780 cube: Cube
1781 2 dimensional Cube of the data to plot as power spectrum.
1782 filename: str
1783 Filename of the plot to write.
1784 title: str
1785 Plot title.
1786 stamp_coordinate: str
1787 Coordinate that becomes different plots.
1788 """
1789 # Use the smallest square grid that will fit the members.
1790 grid_size = int(math.ceil(math.sqrt(len(cube.coord(stamp_coordinate).points))))
1792 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1793 # Make a subplot for each member.
1794 for member, subplot in zip(
1795 cube.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
1796 ):
1797 # Implicit interface is much easier here, due to needing to have the
1798 # cartopy GeoAxes generated.
1799 plt.subplot(grid_size, grid_size, subplot)
1801 frequency = member.coord("frequency").points
1803 axes = plt.gca()
1804 axes.plot(frequency, member.data)
1805 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1807 # Overall figure title.
1808 fig.suptitle(title, fontsize=16)
1810 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1811 logging.info("Saved power spectra postage stamp plot to %s", filename)
1812 plt.close(fig)
1815def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1816 cube: iris.cube.Cube,
1817 filename: str,
1818 title: str,
1819 stamp_coordinate: str,
1820 **kwargs,
1821):
1822 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1823 ax.set_title(title, fontsize=16)
1824 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1825 ax.set_ylabel("Power", fontsize=14)
1826 # Loop over all slices along the stamp_coordinate
1827 for member in cube.slices_over(stamp_coordinate):
1828 frequency = member.coord("frequency").points
1829 ax.plot(
1830 frequency,
1831 member.data,
1832 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1833 )
1835 # Add a legend
1836 ax.legend(fontsize=16)
1838 # Save the figure to a file
1839 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1840 logging.info("Saved power spectra plot to %s", filename)
1842 # Close the figure
1843 plt.close(fig)
1846def _spatial_plot(
1847 method: Literal["contourf", "pcolormesh"],
1848 cube: iris.cube.Cube,
1849 filename: str | None,
1850 sequence_coordinate: str,
1851 stamp_coordinate: str,
1852 overlay_cube: iris.cube.Cube | None = None,
1853 contour_cube: iris.cube.Cube | None = None,
1854 **kwargs,
1855):
1856 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1858 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1859 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1860 is present then postage stamp plots will be produced.
1862 If an overlay_cube and/or contour_cube are specified, multiple variables can
1863 be overplotted on the same figure.
1865 Parameters
1866 ----------
1867 method: "contourf" | "pcolormesh"
1868 The plotting method to use.
1869 cube: Cube
1870 Iris cube of the data to plot. It should have two spatial dimensions,
1871 such as lat and lon, and may also have a another two dimension to be
1872 plotted sequentially and/or as postage stamp plots.
1873 filename: str | None
1874 Name of the plot to write, used as a prefix for plot sequences. If None
1875 uses the recipe name.
1876 sequence_coordinate: str
1877 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1878 This coordinate must exist in the cube.
1879 stamp_coordinate: str
1880 Coordinate about which to plot postage stamp plots. Defaults to
1881 ``"realization"``.
1882 overlay_cube: Cube | None, optional
1883 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1884 contour_cube: Cube | None, optional
1885 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1887 Raises
1888 ------
1889 ValueError
1890 If the cube doesn't have the right dimensions.
1891 TypeError
1892 If the cube isn't a single cube.
1893 """
1894 recipe_title = get_recipe_metadata().get("title", "Untitled")
1896 # Ensure we've got a single cube.
1897 cube = _check_single_cube(cube)
1899 # Test if valid stamp coordinate in cube dimensions.
1900 stamp_coordinate = check_stamp_coordinate(cube)
1902 # Make postage stamp plots if stamp_coordinate exists and has more than a
1903 # single point.
1904 plotting_func = _plot_and_save_spatial_plot
1905 try:
1906 if cube.coord(stamp_coordinate).shape[0] > 1:
1907 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1908 except iris.exceptions.CoordinateNotFoundError:
1909 pass
1911 # Produce a geographical scatter plot if the data have a
1912 # dimension called observation or model_obs_error
1913 if any( 1913 ↛ 1917line 1913 didn't jump to line 1917 because the condition on line 1913 was never true
1914 crd.var_name == "station" or crd.var_name == "model_obs_error"
1915 for crd in cube.coords()
1916 ):
1917 plotting_func = _plot_and_save_scattermap_plot
1919 # Must have a sequence coordinate.
1920 try:
1921 cube.coord(sequence_coordinate)
1922 except iris.exceptions.CoordinateNotFoundError as err:
1923 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1925 # Create a plot for each value of the sequence coordinate.
1926 plot_index = []
1927 nplot = np.size(cube.coord(sequence_coordinate).points)
1929 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1930 # Set plot titles and filename
1931 seq_coord = cube_slice.coord(sequence_coordinate)
1932 plot_title, plot_filename = _set_title_and_filename(
1933 seq_coord, nplot, recipe_title, filename
1934 )
1936 # Extract sequence slice for overlay_cube and contour_cube if required.
1937 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1938 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1940 # Do the actual plotting.
1941 plotting_func(
1942 cube_slice,
1943 filename=plot_filename,
1944 stamp_coordinate=stamp_coordinate,
1945 title=plot_title,
1946 method=method,
1947 overlay_cube=overlay_slice,
1948 contour_cube=contour_slice,
1949 **kwargs,
1950 )
1951 plot_index.append(plot_filename)
1953 # Add list of plots to plot metadata.
1954 complete_plot_index = _append_to_plot_index(plot_index)
1956 # Make a page to display the plots.
1957 _make_plot_html_page(complete_plot_index)
1960def _custom_colormap_mask(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
1961 """Get colourmap for mask.
1963 If "mask_for_" appears anywhere in the name of a cube this function will be called
1964 regardless of the name of the variable to ensure a consistent plot.
1966 Parameters
1967 ----------
1968 cube: Cube
1969 Cube of variable for which the colorbar information is desired.
1970 axis: "x", "y", optional
1971 Select the levels for just this axis of a line plot. The min and max
1972 can be set by xmin/xmax or ymin/ymax respectively. For variables where
1973 setting a universal range is not desirable (e.g. temperature), users
1974 can set ymin/ymax values to "auto" in the colorbar definitions file.
1975 Where no additional xmin/xmax or ymin/ymax values are provided, the
1976 axis bounds default to use the vmin/vmax values provided.
1978 Returns
1979 -------
1980 cmap:
1981 Matplotlib colormap.
1982 levels:
1983 List of levels to use for plotting. For continuous plots the min and max
1984 should be taken as the range.
1985 norm:
1986 BoundaryNorm information.
1987 """
1988 if "difference" not in cube.long_name:
1989 if axis:
1990 levels = [0, 1]
1991 # Complete settings based on levels.
1992 return None, levels, None
1993 else:
1994 # Define the levels and colors.
1995 levels = [0, 1, 2]
1996 colors = ["white", "dodgerblue"]
1997 # Create a custom color map.
1998 cmap = mcolors.ListedColormap(colors)
1999 # Normalize the levels.
2000 norm = mcolors.BoundaryNorm(levels, cmap.N)
2001 logging.debug("Colourmap for %s.", cube.long_name)
2002 return cmap, levels, norm
2003 else:
2004 if axis:
2005 levels = [-1, 1]
2006 return None, levels, None
2007 else:
2008 # Search for if mask difference, set to +/- 0.5 as values plotted <
2009 # not <=.
2010 levels = [-2, -0.5, 0.5, 2]
2011 colors = ["goldenrod", "white", "teal"]
2012 cmap = mcolors.ListedColormap(colors)
2013 norm = mcolors.BoundaryNorm(levels, cmap.N)
2014 logging.debug("Colourmap for %s.", cube.long_name)
2015 return cmap, levels, norm
2018def _custom_beaufort_scale(cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None):
2019 """Get a custom colorbar for a cube in the Beaufort Scale.
2021 Specific variable ranges can be separately set in user-supplied definition
2022 for x- or y-axis limits, or indicate where automated range preferred.
2024 Parameters
2025 ----------
2026 cube: Cube
2027 Cube of variable with Beaufort Scale in name.
2028 axis: "x", "y", optional
2029 Select the levels for just this axis of a line plot. The min and max
2030 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2031 setting a universal range is not desirable (e.g. temperature), users
2032 can set ymin/ymax values to "auto" in the colorbar definitions file.
2033 Where no additional xmin/xmax or ymin/ymax values are provided, the
2034 axis bounds default to use the vmin/vmax values provided.
2036 Returns
2037 -------
2038 cmap:
2039 Matplotlib colormap.
2040 levels:
2041 List of levels to use for plotting. For continuous plots the min and max
2042 should be taken as the range.
2043 norm:
2044 BoundaryNorm information.
2045 """
2046 if "difference" not in cube.long_name:
2047 if axis:
2048 levels = [0, 12]
2049 return None, levels, None
2050 else:
2051 levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
2052 colors = [
2053 "black",
2054 (0, 0, 0.6),
2055 "blue",
2056 "cyan",
2057 "green",
2058 "yellow",
2059 (1, 0.5, 0),
2060 "red",
2061 "pink",
2062 "magenta",
2063 "purple",
2064 "maroon",
2065 "white",
2066 ]
2067 cmap = mcolors.ListedColormap(colors)
2068 norm = mcolors.BoundaryNorm(levels, cmap.N)
2069 logging.info("change colormap for Beaufort Scale colorbar.")
2070 return cmap, levels, norm
2071 else:
2072 if axis:
2073 levels = [-4, 4]
2074 return None, levels, None
2075 else:
2076 levels = [
2077 -3.5,
2078 -2.5,
2079 -1.5,
2080 -0.5,
2081 0.5,
2082 1.5,
2083 2.5,
2084 3.5,
2085 ]
2086 cmap = plt.get_cmap("bwr", 8)
2087 norm = mcolors.BoundaryNorm(levels, cmap.N)
2088 return cmap, levels, norm
2091def _custom_colormap_celsius(cube: iris.cube.Cube, cmap, levels, norm):
2092 """Return altered colourmap for temperature with change in units to Celsius.
2094 If "Celsius" appears anywhere in the name of a cube this function will be called.
2096 Parameters
2097 ----------
2098 cube: Cube
2099 Cube of variable for which the colorbar information is desired.
2100 cmap: Matplotlib colormap.
2101 levels: List
2102 List of levels to use for plotting. For continuous plots the min and max
2103 should be taken as the range.
2104 norm: BoundaryNorm.
2106 Returns
2107 -------
2108 cmap: Matplotlib colormap.
2109 levels: List
2110 List of levels to use for plotting. For continuous plots the min and max
2111 should be taken as the range.
2112 norm: BoundaryNorm.
2113 """
2114 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2115 if any("temperature" in name for name in varnames) and "Celsius" == cube.units:
2116 levels = np.array(levels)
2117 levels -= 273
2118 levels = levels.tolist()
2119 else:
2120 # Do nothing keep the existing colourbar attributes
2121 levels = levels
2122 cmap = cmap
2123 norm = norm
2124 return cmap, levels, norm
2127def _custom_colormap_probability(
2128 cube: iris.cube.Cube, axis: Literal["x", "y"] | None = None
2129):
2130 """Get a custom colorbar for a probability cube.
2132 Specific variable ranges can be separately set in user-supplied definition
2133 for x- or y-axis limits, or indicate where automated range preferred.
2135 Parameters
2136 ----------
2137 cube: Cube
2138 Cube of variable with probability in name.
2139 axis: "x", "y", optional
2140 Select the levels for just this axis of a line plot. The min and max
2141 can be set by xmin/xmax or ymin/ymax respectively. For variables where
2142 setting a universal range is not desirable (e.g. temperature), users
2143 can set ymin/ymax values to "auto" in the colorbar definitions file.
2144 Where no additional xmin/xmax or ymin/ymax values are provided, the
2145 axis bounds default to use the vmin/vmax values provided.
2147 Returns
2148 -------
2149 cmap:
2150 Matplotlib colormap.
2151 levels:
2152 List of levels to use for plotting. For continuous plots the min and max
2153 should be taken as the range.
2154 norm:
2155 BoundaryNorm information.
2156 """
2157 if axis:
2158 levels = [0, 1]
2159 return None, levels, None
2160 else:
2161 cmap = mcolors.ListedColormap(
2162 [
2163 "#FFFFFF",
2164 "#636363",
2165 "#e1dada",
2166 "#B5CAFF",
2167 "#8FB3FF",
2168 "#7F97FF",
2169 "#ABCF63",
2170 "#E8F59E",
2171 "#FFFA14",
2172 "#FFD121",
2173 "#FFA30A",
2174 ]
2175 )
2176 levels = [0.0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
2177 norm = mcolors.BoundaryNorm(levels, cmap.N)
2178 return cmap, levels, norm
2181def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
2182 """Return a custom colourmap for the current recipe."""
2183 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2184 if (
2185 any("surface_microphysical" in name for name in varnames)
2186 and "difference" not in cube.long_name
2187 and "mask" not in cube.long_name
2188 ):
2189 # Define the levels and colors
2190 levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
2191 colors = [
2192 "w",
2193 (0, 0, 0.6),
2194 "b",
2195 "c",
2196 "g",
2197 "y",
2198 (1, 0.5, 0),
2199 "r",
2200 "pink",
2201 "m",
2202 "purple",
2203 "maroon",
2204 "gray",
2205 ]
2206 # Create a custom colormap
2207 cmap = mcolors.ListedColormap(colors)
2208 # Normalize the levels
2209 norm = mcolors.BoundaryNorm(levels, cmap.N)
2210 logging.info("change colormap for surface_microphysical variable colorbar.")
2211 else:
2212 # do nothing and keep existing colorbar attributes
2213 cmap = cmap
2214 levels = levels
2215 norm = norm
2216 return cmap, levels, norm
2219def _custom_colormap_aviation_colour_state(cube: iris.cube.Cube):
2220 """Return custom colourmap for aviation colour state.
2222 If "aviation_colour_state" appears anywhere in the name of a cube
2223 this function will be called.
2225 Parameters
2226 ----------
2227 cube: Cube
2228 Cube of variable for which the colorbar information is desired.
2230 Returns
2231 -------
2232 cmap: Matplotlib colormap.
2233 levels: List
2234 List of levels to use for plotting. For continuous plots the min and max
2235 should be taken as the range.
2236 norm: BoundaryNorm.
2237 """
2238 levels = [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]
2239 colors = [
2240 "#87ceeb",
2241 "#ffffff",
2242 "#8ced69",
2243 "#ffff00",
2244 "#ffd700",
2245 "#ffa500",
2246 "#fe3620",
2247 ]
2248 # Create a custom colormap
2249 cmap = mcolors.ListedColormap(colors)
2250 # Normalise the levels
2251 norm = mcolors.BoundaryNorm(levels, cmap.N)
2252 return cmap, levels, norm
2255def _custom_colourmap_visibility_in_air(cube: iris.cube.Cube, cmap, levels, norm):
2256 """Return a custom colourmap for the current recipe."""
2257 varnames = filter(None, [cube.long_name, cube.standard_name, cube.var_name])
2258 if (
2259 any("visibility_in_air" in name for name in varnames)
2260 and "difference" not in cube.long_name
2261 and "mask" not in cube.long_name
2262 ):
2263 # Define the levels and colors (in km)
2264 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]
2265 norm = mcolors.BoundaryNorm(levels, cmap.N)
2266 colours = [
2267 "#8f00d6",
2268 "#d10000",
2269 "#ff9700",
2270 "#ffff00",
2271 "#00007f",
2272 "#6c9ccd",
2273 "#aae8ff",
2274 "#37a648",
2275 "#8edc64",
2276 "#c5ffc5",
2277 "#dcdcdc",
2278 "#ffffff",
2279 ]
2280 # Create a custom colormap
2281 cmap = mcolors.ListedColormap(colours)
2282 # Normalize the levels
2283 norm = mcolors.BoundaryNorm(levels, cmap.N)
2284 logging.info("change colormap for visibility_in_air variable colorbar.")
2285 else:
2286 # do nothing and keep existing colorbar attributes
2287 cmap = cmap
2288 levels = levels
2289 norm = norm
2290 return cmap, levels, norm
2293def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
2294 """Return number of models based on cube attributes."""
2295 model_names = list(
2296 filter(
2297 lambda x: x is not None,
2298 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
2299 )
2300 )
2301 if not model_names:
2302 logging.debug("Missing model names. Will assume single model.")
2303 return 1
2304 else:
2305 return len(model_names)
2308def _validate_cube_shape(
2309 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
2310) -> None:
2311 """Check all cubes have a model name."""
2312 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 2312 ↛ 2313line 2312 didn't jump to line 2313 because the condition on line 2312 was never true
2313 raise ValueError(
2314 f"The number of model names ({num_models}) should equal the number "
2315 f"of cubes ({len(cube)})."
2316 )
2319def _validate_cubes_coords(
2320 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
2321) -> None:
2322 """Check same number of cubes as sequence coordinate for zip functions."""
2323 if len(cubes) != len(coords): 2323 ↛ 2324line 2323 didn't jump to line 2324 because the condition on line 2323 was never true
2324 raise ValueError(
2325 f"The number of CubeList entries ({len(cubes)}) should equal the number "
2326 f"of sequence coordinates ({len(coords)})."
2327 f"Check that number of time entries in input data are consistent if "
2328 f"performing time-averaging steps prior to plotting outputs."
2329 )
2332####################
2333# Public functions #
2334####################
2337def spatial_contour_plot(
2338 cube: iris.cube.Cube,
2339 filename: str = None,
2340 sequence_coordinate: str = "time",
2341 stamp_coordinate: str = "realization",
2342 **kwargs,
2343) -> iris.cube.Cube:
2344 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2346 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2347 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2348 is present then postage stamp plots will be produced.
2350 Parameters
2351 ----------
2352 cube: Cube
2353 Iris cube of the data to plot. It should have two spatial dimensions,
2354 such as lat and lon, and may also have a another two dimension to be
2355 plotted sequentially and/or as postage stamp plots.
2356 filename: str, optional
2357 Name of the plot to write, used as a prefix for plot sequences. Defaults
2358 to the recipe name.
2359 sequence_coordinate: str, optional
2360 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2361 This coordinate must exist in the cube.
2362 stamp_coordinate: str, optional
2363 Coordinate about which to plot postage stamp plots. Defaults to
2364 ``"realization"``.
2366 Returns
2367 -------
2368 Cube
2369 The original cube (so further operations can be applied).
2371 Raises
2372 ------
2373 ValueError
2374 If the cube doesn't have the right dimensions.
2375 TypeError
2376 If the cube isn't a single cube.
2377 """
2378 _spatial_plot(
2379 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2380 )
2381 return cube
2384def spatial_pcolormesh_plot(
2385 cube: iris.cube.Cube,
2386 filename: str = None,
2387 sequence_coordinate: str = "time",
2388 stamp_coordinate: str = "realization",
2389 **kwargs,
2390) -> iris.cube.Cube:
2391 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2393 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2394 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2395 is present then postage stamp plots will be produced.
2397 This function is significantly faster than ``spatial_contour_plot``,
2398 especially at high resolutions, and should be preferred unless contiguous
2399 contour areas are important.
2401 Parameters
2402 ----------
2403 cube: Cube
2404 Iris cube of the data to plot. It should have two spatial dimensions,
2405 such as lat and lon, and may also have a another two dimension to be
2406 plotted sequentially and/or as postage stamp plots.
2407 filename: str, optional
2408 Name of the plot to write, used as a prefix for plot sequences. Defaults
2409 to the recipe name.
2410 sequence_coordinate: str, optional
2411 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2412 This coordinate must exist in the cube.
2413 stamp_coordinate: str, optional
2414 Coordinate about which to plot postage stamp plots. Defaults to
2415 ``"realization"``.
2417 Returns
2418 -------
2419 Cube
2420 The original cube (so further operations can be applied).
2422 Raises
2423 ------
2424 ValueError
2425 If the cube doesn't have the right dimensions.
2426 TypeError
2427 If the cube isn't a single cube.
2428 """
2429 _spatial_plot(
2430 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2431 )
2432 return cube
2435def spatial_multi_pcolormesh_plot(
2436 cube: iris.cube.Cube,
2437 overlay_cube: iris.cube.Cube,
2438 contour_cube: iris.cube.Cube,
2439 filename: str = None,
2440 sequence_coordinate: str = "time",
2441 stamp_coordinate: str = "realization",
2442 **kwargs,
2443) -> iris.cube.Cube:
2444 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2446 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2447 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2448 is present then postage stamp plots will be produced.
2450 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2452 If specified, contours of a contour_cube can be overplotted on top of those.
2454 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2456 This function is significantly faster than ``spatial_contour_plot``,
2457 especially at high resolutions, and should be preferred unless contiguous
2458 contour areas are important.
2460 Parameters
2461 ----------
2462 cube: Cube
2463 Iris cube of the data to plot. It should have two spatial dimensions,
2464 such as lat and lon, and may also have a another two dimension to be
2465 plotted sequentially and/or as postage stamp plots.
2466 overlay_cube: Cube
2467 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2468 such as lat and lon, and may also have a another two dimension to be
2469 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.
2470 contour_cube: Cube
2471 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,
2472 such as lat and lon, and may also have a another two dimension to be
2473 plotted sequentially and/or as postage stamp plots.
2474 filename: str, optional
2475 Name of the plot to write, used as a prefix for plot sequences. Defaults
2476 to the recipe name.
2477 sequence_coordinate: str, optional
2478 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2479 This coordinate must exist in the cube.
2480 stamp_coordinate: str, optional
2481 Coordinate about which to plot postage stamp plots. Defaults to
2482 ``"realization"``.
2484 Returns
2485 -------
2486 Cube
2487 The original cube (so further operations can be applied).
2489 Raises
2490 ------
2491 ValueError
2492 If the cube doesn't have the right dimensions.
2493 TypeError
2494 If the cube isn't a single cube.
2495 """
2496 _spatial_plot(
2497 "pcolormesh",
2498 cube,
2499 filename,
2500 sequence_coordinate,
2501 stamp_coordinate,
2502 overlay_cube=overlay_cube,
2503 contour_cube=contour_cube,
2504 )
2505 return cube, overlay_cube, contour_cube
2508# TODO: Expand function to handle ensemble data.
2509# line_coordinate: str, optional
2510# Coordinate about which to plot multiple lines. Defaults to
2511# ``"realization"``.
2512def plot_line_series(
2513 cube: iris.cube.Cube | iris.cube.CubeList,
2514 filename: str = None,
2515 series_coordinate: str = "time",
2516 # line_coordinate: str = "realization",
2517 **kwargs,
2518) -> iris.cube.Cube | iris.cube.CubeList:
2519 """Plot a line plot for the specified coordinate.
2521 The Cube or CubeList must be 1D.
2523 Parameters
2524 ----------
2525 iris.cube | iris.cube.CubeList
2526 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2527 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2528 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2529 filename: str, optional
2530 Name of the plot to write, used as a prefix for plot sequences. Defaults
2531 to the recipe name.
2532 series_coordinate: str, optional
2533 Coordinate about which to make a series. Defaults to ``"time"``. This
2534 coordinate must exist in the cube.
2536 Returns
2537 -------
2538 iris.cube.Cube | iris.cube.CubeList
2539 The original Cube or CubeList (so further operations can be applied).
2540 plotted data.
2542 Raises
2543 ------
2544 ValueError
2545 If the cubes don't have the right dimensions.
2546 TypeError
2547 If the cube isn't a Cube or CubeList.
2548 """
2549 # Ensure we have a name for the plot file.
2550 recipe_title = get_recipe_metadata().get("title", "Untitled")
2552 num_models = _get_num_models(cube)
2554 _validate_cube_shape(cube, num_models)
2556 # Iterate over all cubes and extract coordinate to plot.
2557 cubes = iter_maybe(cube)
2558 coords = []
2559 for cube in cubes:
2560 try:
2561 coords.append(cube.coord(series_coordinate))
2562 except iris.exceptions.CoordinateNotFoundError as err:
2563 raise ValueError(
2564 f"Cube must have a {series_coordinate} coordinate."
2565 ) from err
2566 if cube.ndim > 2 or not cube.coords("realization"):
2567 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2569 # Format the title and filename using plotted series coordinate
2570 nplot = 1
2571 seq_coord = coords[0]
2572 plot_title, plot_filename = _set_title_and_filename(
2573 seq_coord, nplot, recipe_title, filename
2574 )
2576 # Do the actual plotting.
2577 _plot_and_save_line_series(cubes, coords, "realization", plot_filename, plot_title)
2579 # Add list of plots to plot metadata.
2580 plot_index = _append_to_plot_index([plot_filename])
2582 # Make a page to display the plots.
2583 _make_plot_html_page(plot_index)
2585 return cube
2588def plot_vertical_line_series(
2589 cubes: iris.cube.Cube | iris.cube.CubeList,
2590 filename: str = None,
2591 series_coordinate: str = "model_level_number",
2592 sequence_coordinate: str = "time",
2593 # line_coordinate: str = "realization",
2594 **kwargs,
2595) -> iris.cube.Cube | iris.cube.CubeList:
2596 """Plot a line plot against a type of vertical coordinate.
2598 The Cube or CubeList must be 1D.
2600 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2601 then a sequence of plots will be produced.
2603 Parameters
2604 ----------
2605 iris.cube | iris.cube.CubeList
2606 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2607 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2608 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2609 filename: str, optional
2610 Name of the plot to write, used as a prefix for plot sequences. Defaults
2611 to the recipe name.
2612 series_coordinate: str, optional
2613 Coordinate to plot on the y-axis. Can be ``pressure`` or
2614 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2615 for LFRic. Defaults to ``model_level_number``.
2616 This coordinate must exist in the cube.
2617 sequence_coordinate: str, optional
2618 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2619 This coordinate must exist in the cube.
2621 Returns
2622 -------
2623 iris.cube.Cube | iris.cube.CubeList
2624 The original Cube or CubeList (so further operations can be applied).
2625 Plotted data.
2627 Raises
2628 ------
2629 ValueError
2630 If the cubes doesn't have the right dimensions.
2631 TypeError
2632 If the cube isn't a Cube or CubeList.
2633 """
2634 # Ensure we have a name for the plot file.
2635 recipe_title = get_recipe_metadata().get("title", "Untitled")
2637 cubes = iter_maybe(cubes)
2638 # Initialise empty list to hold all data from all cubes in a CubeList
2639 all_data = []
2641 # Store min/max ranges for x range.
2642 x_levels = []
2644 num_models = _get_num_models(cubes)
2646 _validate_cube_shape(cubes, num_models)
2648 # Iterate over all cubes in cube or CubeList and plot.
2649 coords = []
2650 for cube in cubes:
2651 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2652 try:
2653 coords.append(cube.coord(series_coordinate))
2654 except iris.exceptions.CoordinateNotFoundError as err:
2655 raise ValueError(
2656 f"Cube must have a {series_coordinate} coordinate."
2657 ) from err
2659 try:
2660 if cube.ndim > 1 or not cube.coords("realization"): 2660 ↛ 2668line 2660 didn't jump to line 2668 because the condition on line 2660 was always true
2661 cube.coord(sequence_coordinate)
2662 except iris.exceptions.CoordinateNotFoundError as err:
2663 raise ValueError(
2664 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2665 ) from err
2667 # Get minimum and maximum from levels information.
2668 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2669 if levels is not None: 2669 ↛ 2673line 2669 didn't jump to line 2673 because the condition on line 2669 was always true
2670 x_levels.append(min(levels))
2671 x_levels.append(max(levels))
2672 else:
2673 all_data.append(cube.data)
2675 if len(x_levels) == 0: 2675 ↛ 2677line 2675 didn't jump to line 2677 because the condition on line 2675 was never true
2676 # Combine all data into a single NumPy array
2677 combined_data = np.concatenate(all_data)
2679 # Set the lower and upper limit for the x-axis to ensure all plots have
2680 # same range. This needs to read the whole cube over the range of the
2681 # sequence and if applicable postage stamp coordinate.
2682 vmin = np.floor(combined_data.min())
2683 vmax = np.ceil(combined_data.max())
2684 else:
2685 vmin = min(x_levels)
2686 vmax = max(x_levels)
2688 # Matching the slices (matching by seq coord point; it may happen that
2689 # evaluated models do not cover the same seq coord range, hence matching
2690 # necessary)
2691 def filter_cube_iterables(cube_iterables) -> bool:
2692 return len(cube_iterables) == len(coords)
2694 cube_iterables = filter(
2695 filter_cube_iterables,
2696 (
2697 iris.cube.CubeList(
2698 s
2699 for s in itertools.chain.from_iterable(
2700 cb.slices_over(sequence_coordinate) for cb in cubes
2701 )
2702 if s.coord(sequence_coordinate).points[0] == point
2703 )
2704 for point in sorted(
2705 set(
2706 itertools.chain.from_iterable(
2707 cb.coord(sequence_coordinate).points for cb in cubes
2708 )
2709 )
2710 )
2711 ),
2712 )
2714 # Create a plot for each value of the sequence coordinate.
2715 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2716 # similar sequence values. Passing a CubeList into the internal plotting function
2717 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2718 # or an iris.cube.CubeList.
2719 plot_index = []
2720 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2721 for cubes_slice in cube_iterables:
2722 # Format the coordinate value in a unit appropriate way.
2723 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2724 plot_title, plot_filename = _set_title_and_filename(
2725 seq_coord, nplot, recipe_title, filename
2726 )
2728 # Do the actual plotting.
2729 _plot_and_save_vertical_line_series(
2730 cubes_slice,
2731 coords,
2732 "realization",
2733 plot_filename,
2734 series_coordinate,
2735 title=plot_title,
2736 vmin=vmin,
2737 vmax=vmax,
2738 )
2739 plot_index.append(plot_filename)
2741 # Add list of plots to plot metadata.
2742 complete_plot_index = _append_to_plot_index(plot_index)
2744 # Make a page to display the plots.
2745 _make_plot_html_page(complete_plot_index)
2747 return cubes
2750def qq_plot(
2751 cubes: iris.cube.CubeList,
2752 coordinates: list[str],
2753 percentiles: list[float],
2754 model_names: list[str],
2755 filename: str = None,
2756 one_to_one: bool = True,
2757 **kwargs,
2758) -> iris.cube.CubeList:
2759 """Plot a Quantile-Quantile plot between two models for common time points.
2761 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2762 collapsed within the operator over all specified coordinates such as
2763 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2764 ensemble members to ensure a 1D cube (array).
2766 Parameters
2767 ----------
2768 cubes: iris.cube.CubeList
2769 Two cubes of the same variable with different models.
2770 coordinate: list[str]
2771 The list of coordinates to collapse over. This list should be
2772 every coordinate within the cube to result in a 1D cube around
2773 the percentile coordinate.
2774 percent: list[float]
2775 A list of percentiles to appear in the plot.
2776 model_names: list[str]
2777 A list of model names to appear on the axis of the plot.
2778 filename: str, optional
2779 Filename of the plot to write.
2780 one_to_one: bool, optional
2781 If True a 1:1 line is plotted; if False it is not. Default is True.
2783 Raises
2784 ------
2785 ValueError
2786 When the cubes are not compatible.
2788 Notes
2789 -----
2790 The quantile-quantile plot is a variant on the scatter plot representing
2791 two datasets by their quantiles (percentiles) for common time points.
2792 This plot does not use a theoretical distribution to compare against, but
2793 compares percentiles of two datasets. This plot does
2794 not use all raw data points, but plots the selected percentiles (quantiles) of
2795 each variable instead for the two datasets, thereby normalising the data for a
2796 direct comparison between the selected percentiles of the two dataset distributions.
2798 Quantile-quantile plots are valuable for comparing against
2799 observations and other models. Identical percentiles between the variables
2800 will lie on the one-to-one line implying the values correspond well to each
2801 other. Where there is a deviation from the one-to-one line a range of
2802 possibilities exist depending on how and where the data is shifted (e.g.,
2803 Wilks 2011 [Wilks2011]_).
2805 For distributions above the one-to-one line the distribution is left-skewed;
2806 below is right-skewed. A distinct break implies a bimodal distribution, and
2807 closer values/values further apart at the tails imply poor representation of
2808 the extremes.
2810 References
2811 ----------
2812 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2813 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2814 """
2815 # Check cubes using same functionality as the difference operator.
2816 if len(cubes) != 2:
2817 raise ValueError("cubes should contain exactly 2 cubes.")
2818 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2819 other: Cube = cubes.extract_cube(
2820 iris.Constraint(
2821 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2822 )
2823 )
2825 # Get spatial coord names.
2826 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2827 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2829 # Ensure cubes to compare are on common differencing grid.
2830 # This is triggered if either
2831 # i) latitude and longitude shapes are not the same. Note grid points
2832 # are not compared directly as these can differ through rounding
2833 # errors.
2834 # ii) or variables are known to often sit on different grid staggering
2835 # in different models (e.g. cell center vs cell edge), as is the case
2836 # for UM and LFRic comparisons.
2837 # In future greater choice of regridding method might be applied depending
2838 # on variable type. Linear regridding can in general be appropriate for smooth
2839 # variables. Care should be taken with interpretation of differences
2840 # given this dependency on regridding.
2841 if (
2842 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2843 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2844 ) or (
2845 base.long_name
2846 in [
2847 "eastward_wind_at_10m",
2848 "northward_wind_at_10m",
2849 "northward_wind_at_cell_centres",
2850 "eastward_wind_at_cell_centres",
2851 "zonal_wind_at_pressure_levels",
2852 "meridional_wind_at_pressure_levels",
2853 "potential_vorticity_at_pressure_levels",
2854 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2855 ]
2856 ):
2857 logging.debug(
2858 "Linear regridding base cube to other grid to compute differences"
2859 )
2860 base = regrid_onto_cube(base, other, method="Linear")
2862 # Extract just common time points.
2863 base, other = _extract_common_time_points(base, other)
2865 # Equalise attributes so we can merge.
2866 fully_equalise_attributes([base, other])
2867 logging.debug("Base: %s\nOther: %s", base, other)
2869 # Collapse cubes.
2870 base = collapse(
2871 base,
2872 coordinate=coordinates,
2873 method="PERCENTILE",
2874 additional_percent=percentiles,
2875 )
2876 other = collapse(
2877 other,
2878 coordinate=coordinates,
2879 method="PERCENTILE",
2880 additional_percent=percentiles,
2881 )
2883 # Ensure we have a name for the plot file.
2884 recipe_title = get_recipe_metadata().get("title", "Untitled")
2885 title = f"{recipe_title}"
2887 if filename is None:
2888 filename = slugify(recipe_title)
2890 # Add file extension.
2891 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2893 # Do the actual plotting on a scatter plot
2894 _plot_and_save_scatter_plot(
2895 base, other, plot_filename, title, one_to_one, model_names
2896 )
2898 # Add list of plots to plot metadata.
2899 plot_index = _append_to_plot_index([plot_filename])
2901 # Make a page to display the plots.
2902 _make_plot_html_page(plot_index)
2904 return iris.cube.CubeList([base, other])
2907def scatter_plot(
2908 cube_x: iris.cube.Cube | iris.cube.CubeList,
2909 cube_y: iris.cube.Cube | iris.cube.CubeList,
2910 filename: str = None,
2911 one_to_one: bool = True,
2912 **kwargs,
2913) -> iris.cube.CubeList:
2914 """Plot a scatter plot between two variables.
2916 Both cubes must be 1D.
2918 Parameters
2919 ----------
2920 cube_x: Cube | CubeList
2921 1 dimensional Cube of the data to plot on y-axis.
2922 cube_y: Cube | CubeList
2923 1 dimensional Cube of the data to plot on x-axis.
2924 filename: str, optional
2925 Filename of the plot to write.
2926 one_to_one: bool, optional
2927 If True a 1:1 line is plotted; if False it is not. Default is True.
2929 Returns
2930 -------
2931 cubes: CubeList
2932 CubeList of the original x and y cubes for further processing.
2934 Raises
2935 ------
2936 ValueError
2937 If the cube doesn't have the right dimensions and cubes not the same
2938 size.
2939 TypeError
2940 If the cube isn't a single cube.
2942 Notes
2943 -----
2944 Scatter plots are used for determining if there is a relationship between
2945 two variables. Positive relations have a slope going from bottom left to top
2946 right; Negative relations have a slope going from top left to bottom right.
2947 """
2948 # Iterate over all cubes in cube or CubeList and plot.
2949 for cube_iter in iter_maybe(cube_x):
2950 # Check cubes are correct shape.
2951 cube_iter = _check_single_cube(cube_iter)
2952 if cube_iter.ndim > 1:
2953 raise ValueError("cube_x must be 1D.")
2955 # Iterate over all cubes in cube or CubeList and plot.
2956 for cube_iter in iter_maybe(cube_y):
2957 # Check cubes are correct shape.
2958 cube_iter = _check_single_cube(cube_iter)
2959 if cube_iter.ndim > 1:
2960 raise ValueError("cube_y must be 1D.")
2962 # Ensure we have a name for the plot file.
2963 recipe_title = get_recipe_metadata().get("title", "Untitled")
2964 title = f"{recipe_title}"
2966 if filename is None:
2967 filename = slugify(recipe_title)
2969 # Add file extension.
2970 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2972 # Do the actual plotting.
2973 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2975 # Add list of plots to plot metadata.
2976 plot_index = _append_to_plot_index([plot_filename])
2978 # Make a page to display the plots.
2979 _make_plot_html_page(plot_index)
2981 return iris.cube.CubeList([cube_x, cube_y])
2984def vector_plot(
2985 cube_u: iris.cube.Cube,
2986 cube_v: iris.cube.Cube,
2987 filename: str = None,
2988 sequence_coordinate: str = "time",
2989 **kwargs,
2990) -> iris.cube.CubeList:
2991 """Plot a vector plot based on the input u and v components."""
2992 recipe_title = get_recipe_metadata().get("title", "Untitled")
2994 # Cubes must have a matching sequence coordinate.
2995 try:
2996 # Check that the u and v cubes have the same sequence coordinate.
2997 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2997 ↛ anywhereline 2997 didn't jump anywhere: it always raised an exception.
2998 raise ValueError("Coordinates do not match.")
2999 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3000 raise ValueError(
3001 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3002 ) from err
3004 # Create a plot for each value of the sequence coordinate.
3005 plot_index = []
3006 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3007 for cube_u_slice, cube_v_slice in zip(
3008 cube_u.slices_over(sequence_coordinate),
3009 cube_v.slices_over(sequence_coordinate),
3010 strict=True,
3011 ):
3012 # Format the coordinate value in a unit appropriate way.
3013 seq_coord = cube_u_slice.coord(sequence_coordinate)
3014 plot_title, plot_filename = _set_title_and_filename(
3015 seq_coord, nplot, recipe_title, filename
3016 )
3018 # Do the actual plotting.
3019 _plot_and_save_vector_plot(
3020 cube_u_slice,
3021 cube_v_slice,
3022 filename=plot_filename,
3023 title=plot_title,
3024 method="contourf",
3025 )
3026 plot_index.append(plot_filename)
3028 # Add list of plots to plot metadata.
3029 complete_plot_index = _append_to_plot_index(plot_index)
3031 # Make a page to display the plots.
3032 _make_plot_html_page(complete_plot_index)
3034 return iris.cube.CubeList([cube_u, cube_v])
3037def plot_histogram_series(
3038 cubes: iris.cube.Cube | iris.cube.CubeList,
3039 filename: str = None,
3040 sequence_coordinate: str = "time",
3041 stamp_coordinate: str = "realization",
3042 single_plot: bool = False,
3043 **kwargs,
3044) -> iris.cube.Cube | iris.cube.CubeList:
3045 """Plot a histogram plot for each vertical level provided.
3047 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3048 is present then a sequence of plots will be produced using the time slider
3049 functionality to scroll through histograms against time. If a
3050 stamp_coordinate is present then postage stamp plots will be produced. If
3051 stamp_coordinate and single_plot is True, all postage stamp plots will be
3052 plotted in a single plot instead of separate postage stamp plots.
3054 Parameters
3055 ----------
3056 cubes: Cube | iris.cube.CubeList
3057 Iris cube or CubeList of the data to plot. It should have a single dimension other
3058 than the stamp coordinate.
3059 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3060 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3061 filename: str, optional
3062 Name of the plot to write, used as a prefix for plot sequences. Defaults
3063 to the recipe name.
3064 sequence_coordinate: str, optional
3065 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3066 This coordinate must exist in the cube and will be used for the time
3067 slider.
3068 stamp_coordinate: str, optional
3069 Coordinate about which to plot postage stamp plots. Defaults to
3070 ``"realization"``.
3071 single_plot: bool, optional
3072 If True, all postage stamp plots will be plotted in a single plot. If
3073 False, each postage stamp plot will be plotted separately. Is only valid
3074 if stamp_coordinate exists and has more than a single point.
3076 Returns
3077 -------
3078 iris.cube.Cube | iris.cube.CubeList
3079 The original Cube or CubeList (so further operations can be applied).
3080 Plotted data.
3082 Raises
3083 ------
3084 ValueError
3085 If the cube doesn't have the right dimensions.
3086 TypeError
3087 If the cube isn't a Cube or CubeList.
3088 """
3089 recipe_title = get_recipe_metadata().get("title", "Untitled")
3091 cubes = iter_maybe(cubes)
3093 # Internal plotting function.
3094 plotting_func = _plot_and_save_histogram_series
3096 num_models = _get_num_models(cubes)
3098 _validate_cube_shape(cubes, num_models)
3100 # If several histograms are plotted with time as sequence_coordinate for the
3101 # time slider option.
3102 for cube in cubes:
3103 try:
3104 cube.coord(sequence_coordinate)
3105 except iris.exceptions.CoordinateNotFoundError as err:
3106 raise ValueError(
3107 f"Cube must have a {sequence_coordinate} coordinate."
3108 ) from err
3110 # Get minimum and maximum from levels information.
3111 levels = None
3112 for cube in cubes: 3112 ↛ 3128line 3112 didn't jump to line 3128 because the loop on line 3112 didn't complete
3113 # First check if user-specified "auto" range variable.
3114 # This maintains the value of levels as None, so proceed.
3115 _, levels, _ = _colorbar_map_levels(cube, axis="y")
3116 if levels is None:
3117 break
3118 # If levels is changed, recheck to use the vmin,vmax or
3119 # levels-based ranges for histogram plots.
3120 _, levels, _ = _colorbar_map_levels(cube)
3121 logging.debug("levels: %s", levels)
3122 if levels is not None: 3122 ↛ 3112line 3122 didn't jump to line 3112 because the condition on line 3122 was always true
3123 vmin = min(levels)
3124 vmax = max(levels)
3125 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
3126 break
3128 if levels is None:
3129 vmin = min(cb.data.min() for cb in cubes)
3130 vmax = max(cb.data.max() for cb in cubes)
3132 # Make postage stamp plots if stamp_coordinate exists and has more than a
3133 # single point. If single_plot is True:
3134 # -- all postage stamp plots will be plotted in a single plot instead of
3135 # separate postage stamp plots.
3136 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3137 # produced per single model only
3138 if num_models == 1: 3138 ↛ 3151line 3138 didn't jump to line 3151 because the condition on line 3138 was always true
3139 if ( 3139 ↛ 3143line 3139 didn't jump to line 3143 because the condition on line 3139 was never true
3140 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3141 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3142 ):
3143 if single_plot:
3144 plotting_func = (
3145 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3146 )
3147 else:
3148 plotting_func = _plot_and_save_postage_stamp_histogram_series
3149 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3150 else:
3151 all_points = sorted(
3152 set(
3153 itertools.chain.from_iterable(
3154 cb.coord(sequence_coordinate).points for cb in cubes
3155 )
3156 )
3157 )
3158 all_slices = list(
3159 itertools.chain.from_iterable(
3160 cb.slices_over(sequence_coordinate) for cb in cubes
3161 )
3162 )
3163 # Matched slices (matched by seq coord point; it may happen that
3164 # evaluated models do not cover the same seq coord range, hence matching
3165 # necessary)
3166 cube_iterables = [
3167 iris.cube.CubeList(
3168 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3169 )
3170 for point in all_points
3171 ]
3173 plot_index = []
3174 nplot = np.size(cube.coord(sequence_coordinate).points)
3175 # Create a plot for each value of the sequence coordinate. Allowing for
3176 # multiple cubes in a CubeList to be plotted in the same plot for similar
3177 # sequence values. Passing a CubeList into the internal plotting function
3178 # for similar values of the sequence coordinate. cube_slice can be an
3179 # iris.cube.Cube or an iris.cube.CubeList.
3180 for cube_slice in cube_iterables:
3181 single_cube = cube_slice
3182 if isinstance(cube_slice, iris.cube.CubeList): 3182 ↛ 3183line 3182 didn't jump to line 3183 because the condition on line 3182 was never true
3183 single_cube = cube_slice[0]
3185 # Set plot titles and filename, based on sequence coordinate
3186 seq_coord = single_cube.coord(sequence_coordinate)
3187 # Use time coordinate in title and filename if single histogram output.
3188 if sequence_coordinate == "realization" and nplot == 1: 3188 ↛ 3189line 3188 didn't jump to line 3189 because the condition on line 3188 was never true
3189 seq_coord = single_cube.coord("time")
3190 plot_title, plot_filename = _set_title_and_filename(
3191 seq_coord, nplot, recipe_title, filename
3192 )
3194 # Do the actual plotting.
3195 plotting_func(
3196 cube_slice,
3197 filename=plot_filename,
3198 stamp_coordinate=stamp_coordinate,
3199 title=plot_title,
3200 vmin=vmin,
3201 vmax=vmax,
3202 )
3203 plot_index.append(plot_filename)
3205 # Add list of plots to plot metadata.
3206 complete_plot_index = _append_to_plot_index(plot_index)
3208 # Make a page to display the plots.
3209 _make_plot_html_page(complete_plot_index)
3211 return cubes
3214def plot_power_spectrum_series(
3215 cubes: iris.cube.Cube | iris.cube.CubeList,
3216 filename: str = None,
3217 sequence_coordinate: str = "time",
3218 stamp_coordinate: str = "realization",
3219 single_plot: bool = False,
3220 **kwargs,
3221) -> iris.cube.Cube | iris.cube.CubeList:
3222 """Plot a power spectrum plot for each vertical level provided.
3224 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
3225 is present then a sequence of plots will be produced using the time slider
3226 functionality to scroll through power spectrum against time. If a
3227 stamp_coordinate is present then postage stamp plots will be produced. If
3228 stamp_coordinate and single_plot is True, all postage stamp plots will be
3229 plotted in a single plot instead of separate postage stamp plots.
3231 Parameters
3232 ----------
3233 cubes: Cube | iris.cube.CubeList
3234 Iris cube or CubeList of the data to plot. It should have a single dimension other
3235 than the stamp coordinate.
3236 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3237 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3238 filename: str, optional
3239 Name of the plot to write, used as a prefix for plot sequences. Defaults
3240 to the recipe name.
3241 sequence_coordinate: str, optional
3242 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3243 This coordinate must exist in the cube and will be used for the time
3244 slider.
3245 stamp_coordinate: str, optional
3246 Coordinate about which to plot postage stamp plots. Defaults to
3247 ``"realization"``.
3248 single_plot: bool, optional
3249 If True, all postage stamp plots will be plotted in a single plot. If
3250 False, each postage stamp plot will be plotted separately. Is only valid
3251 if stamp_coordinate exists and has more than a single point.
3253 Returns
3254 -------
3255 iris.cube.Cube | iris.cube.CubeList
3256 The original Cube or CubeList (so further operations can be applied).
3257 Plotted data.
3259 Raises
3260 ------
3261 ValueError
3262 If the cube doesn't have the right dimensions.
3263 TypeError
3264 If the cube isn't a Cube or CubeList.
3265 """
3266 recipe_title = get_recipe_metadata().get("title", "Untitled")
3268 cubes = iter_maybe(cubes)
3270 # Internal plotting function.
3271 plotting_func = _plot_and_save_power_spectrum_series
3273 num_models = _get_num_models(cubes)
3275 _validate_cube_shape(cubes, num_models)
3277 # If several power spectra are plotted with time as sequence_coordinate for the
3278 # time slider option.
3279 for cube in cubes:
3280 try:
3281 cube.coord(sequence_coordinate)
3282 except iris.exceptions.CoordinateNotFoundError as err:
3283 raise ValueError(
3284 f"Cube must have a {sequence_coordinate} coordinate."
3285 ) from err
3287 # Make postage stamp plots if stamp_coordinate exists and has more than a
3288 # single point. If single_plot is True:
3289 # -- all postage stamp plots will be plotted in a single plot instead of
3290 # separate postage stamp plots.
3291 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3292 # produced per single model only
3293 if num_models == 1: 3293 ↛ 3306line 3293 didn't jump to line 3306 because the condition on line 3293 was always true
3294 if ( 3294 ↛ 3298line 3294 didn't jump to line 3298 because the condition on line 3294 was never true
3295 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3296 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3297 ):
3298 if single_plot:
3299 plotting_func = (
3300 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3301 )
3302 else:
3303 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3304 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3305 else:
3306 all_points = sorted(
3307 set(
3308 itertools.chain.from_iterable(
3309 cb.coord(sequence_coordinate).points for cb in cubes
3310 )
3311 )
3312 )
3313 all_slices = list(
3314 itertools.chain.from_iterable(
3315 cb.slices_over(sequence_coordinate) for cb in cubes
3316 )
3317 )
3318 # Matched slices (matched by seq coord point; it may happen that
3319 # evaluated models do not cover the same seq coord range, hence matching
3320 # necessary)
3321 cube_iterables = [
3322 iris.cube.CubeList(
3323 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
3324 )
3325 for point in all_points
3326 ]
3328 plot_index = []
3329 nplot = np.size(cube.coord(sequence_coordinate).points)
3330 # Create a plot for each value of the sequence coordinate. Allowing for
3331 # multiple cubes in a CubeList to be plotted in the same plot for similar
3332 # sequence values. Passing a CubeList into the internal plotting function
3333 # for similar values of the sequence coordinate. cube_slice can be an
3334 # iris.cube.Cube or an iris.cube.CubeList.
3335 for cube_slice in cube_iterables:
3336 single_cube = cube_slice
3337 if isinstance(cube_slice, iris.cube.CubeList): 3337 ↛ 3338line 3337 didn't jump to line 3338 because the condition on line 3337 was never true
3338 single_cube = cube_slice[0]
3340 # Set plot title and filenames based on sequence values
3341 seq_coord = single_cube.coord(sequence_coordinate)
3342 plot_title, plot_filename = _set_title_and_filename(
3343 seq_coord, nplot, recipe_title, filename
3344 )
3346 # Do the actual plotting.
3347 plotting_func(
3348 cube_slice,
3349 filename=plot_filename,
3350 stamp_coordinate=stamp_coordinate,
3351 title=plot_title,
3352 )
3353 plot_index.append(plot_filename)
3355 # Add list of plots to plot metadata.
3356 complete_plot_index = _append_to_plot_index(plot_index)
3358 # Make a page to display the plots.
3359 _make_plot_html_page(complete_plot_index)
3361 return cubes
3364def _DCT_ps(y_3d):
3365 """Calculate power spectra for regional domains.
3367 Parameters
3368 ----------
3369 y_3d: 3D array
3370 3 dimensional array to calculate spectrum for.
3371 (2D field data with 3rd dimension of time)
3373 Returns: ps_array
3374 Array of power spectra values calculated for input field (for each time)
3376 Method for regional domains:
3377 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3378 as described in Denis et al 2002 [Denis_etal_2002]_.
3380 References
3381 ----------
3382 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3383 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3384 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3385 Monthly Weather Review, Vol. 130, 1812-1828
3386 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3387 """
3388 Nt, Ny, Nx = y_3d.shape
3390 # Max coefficient
3391 Nmin = min(Nx - 1, Ny - 1)
3393 # Create alpha matrix (of wavenumbers)
3394 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3396 # Prepare output array
3397 ps_array = np.zeros((Nt, Nmin))
3399 # Loop over time to get spectrum for each time.
3400 for t in range(Nt):
3401 y_2d = y_3d[t]
3403 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3404 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3405 # cosine basis functions at different spatial frequencies.
3407 # normalise spectrum to allow comparison between models.
3408 fkk = fft.dctn(y_2d, norm="ortho")
3410 # Normalise fkk
3411 fkk = fkk / np.sqrt(Ny * Nx)
3413 # calculate variance of spectral coefficient
3414 sigma_2 = fkk**2 / Nx / Ny
3416 # Group ellipses of alphas into the same wavenumber k/Nmin
3417 for k in range(1, Nmin + 1):
3418 alpha = k / Nmin
3419 alpha_p1 = (k + 1) / Nmin
3421 # Sum up elements matching k
3422 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3423 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3425 return ps_array
3428def _create_alpha_matrix(Ny, Nx):
3429 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3431 Parameters
3432 ----------
3433 Ny, Nx:
3434 Dimensions of the 2D field for which the power spectra is calculated. Used to
3435 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3436 single-scale parameter.
3438 Returns: alpha_matrix
3439 normalisation of 2D wavenumber axes, transforming the spectral domain into
3440 an elliptic coordinate system.
3442 """
3443 # Create x_indices: each row is [1, 2, ..., Nx]
3444 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3446 # Create y_indices: each column is [1, 2, ..., Ny]
3447 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3449 # Compute alpha_matrix
3450 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3452 return alpha_matrix