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