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