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