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