Coverage for src/CSET/operators/plot.py: 83%
1103 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-07 15:12 +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 importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24import sys
25from typing import Literal
27import cartopy.crs as ccrs
28import cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.pyplot as plt
36import numpy as np
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._colormaps import (
50 colorbar_map_levels,
51 get_model_colors_map,
52)
53from CSET.operators._utils import (
54 check_sequence_coordinate,
55 check_single_cube,
56 check_stamp_coordinate,
57 fully_equalise_attributes,
58 get_cube_yxcoordname,
59 get_num_models,
60 is_transect,
61 slice_over_maybe,
62 validate_cube_shape,
63 validate_cubes_coords,
64)
65from CSET.operators.collapse import collapse
66from CSET.operators.misc import _extract_common_time_points
67from CSET.operators.regrid import regrid_onto_cube
69logger = logging.getLogger(__name__)
71# Use a non-interactive plotting backend.
72mpl.use("agg")
75############################
76# Private helper functions #
77############################
80def in_sphinx_gallery():
81 """Test if running plot code in sphinx-gallery context."""
82 return "sphinx_gallery" in sys.modules
85def _append_to_plot_index(plot_index: list) -> list:
86 """Add plots into the plot index, returning the complete plot index."""
87 with open("meta.json", "r+t", encoding="UTF-8") as fp:
88 fcntl.flock(fp, fcntl.LOCK_EX)
89 fp.seek(0)
90 meta = json.load(fp)
91 complete_plot_index = meta.get("plots", [])
92 complete_plot_index = complete_plot_index + plot_index
93 meta["plots"] = complete_plot_index
94 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
95 os.getenv("DO_CASE_AGGREGATION")
96 ):
97 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
98 fp.seek(0)
99 fp.truncate()
100 json.dump(meta, fp, indent=2)
101 return complete_plot_index
104def _make_plot_html_page(plots: list):
105 """Create a HTML page to display a plot image."""
106 # Debug check that plots actually contains some strings.
107 assert isinstance(plots[0], str)
109 # Load HTML template file.
110 operator_files = importlib.resources.files()
111 template_file = operator_files.joinpath("_plot_page_template.html")
113 # Get some metadata.
114 meta = get_recipe_metadata()
115 title = meta.get("title", "Untitled")
116 description = MarkdownIt().render(meta.get("description", "*No description.*"))
118 # Prepare template variables.
119 variables = {
120 "title": title,
121 "description": description,
122 "initial_plot": plots[0],
123 "plots": plots,
124 "title_slug": slugify(title),
125 }
127 # Render template.
128 html = render_file(template_file, **variables)
130 # Save completed HTML.
131 with open("index.html", "wt", encoding="UTF-8") as fp:
132 fp.write(html)
135def _save_close_figure(figure, plot_type: str, filename: str):
136 """Save generated plot figure file and close figure.
138 If running documentation gallery generation, avoid saving to file.
140 Parameters
141 ----------
142 figure:
143 Matplotlib Figure object holding all plot elements.
144 plot_type: str
145 String identifier for plot type for logging information.
146 filename: str
147 Filename for saved figure.
148 """
149 if not in_sphinx_gallery():
150 figure.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
151 logger.info("Saved %s plot to %s", plot_type, filename)
152 plt.close(figure)
155def _setup_spatial_map(
156 cube: iris.cube.Cube,
157 figure,
158 cmap,
159 grid_size: tuple[int, int] | None = None,
160 subplot: int | None = None,
161):
162 """Define map projections, extent and add coastlines and borderlines for spatial plots.
164 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
165 is specified, and map extent defined based on the input data.
167 Parameters
168 ----------
169 cube: Cube
170 2 dimensional (lat and lon) Cube of the data to plot.
171 figure:
172 Matplotlib Figure object holding all plot elements.
173 cmap:
174 Matplotlib colormap.
175 grid_size: (int, int), optional
176 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
177 subplot: int, optional
178 Subplot index if multiple spatial subplots in figure.
180 Returns
181 -------
182 axes:
183 Matplotlib GeoAxes definition.
184 """
185 # Identify min/max plot bounds.
186 try:
187 lat_axis, lon_axis = get_cube_yxcoordname(cube)
188 xmin = np.nanmin(cube.coord(lon_axis).points)
189 xmax = np.nanmax(cube.coord(lon_axis).points)
190 ymin = np.nanmin(cube.coord(lat_axis).points)
191 ymax = np.nanmax(cube.coord(lat_axis).points)
193 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
194 if np.abs(xmax - xmin) > 180.0:
195 xmin = xmin - 180.0
196 xmax = xmax - 180.0
197 logger.debug("Adjusting plot bounds to fit global extent.")
199 # Consider map projection orientation.
200 # Adapting orientation enables plotting across international dateline.
201 # Users can adapt the default central_longitude if alternative projections views.
202 if xmax > 180.0 or xmin < -180.0:
203 central_longitude = 180.0
204 else:
205 central_longitude = 0.0
207 # Define spatial map projection.
208 coord_system = cube.coord(lat_axis).coord_system
209 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
210 # Define rotated pole map projection for rotated pole inputs.
211 projection = ccrs.RotatedPole(
212 pole_longitude=coord_system.grid_north_pole_longitude,
213 pole_latitude=coord_system.grid_north_pole_latitude,
214 central_rotated_longitude=central_longitude,
215 )
216 crs = projection
217 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 217 ↛ 219line 217 didn't jump to line 219 because the condition on line 217 was never true
218 # Define Transverse Mercator projection for TM inputs.
219 projection = ccrs.TransverseMercator(
220 central_longitude=coord_system.longitude_of_central_meridian,
221 central_latitude=coord_system.latitude_of_projection_origin,
222 false_easting=coord_system.false_easting,
223 false_northing=coord_system.false_northing,
224 scale_factor=coord_system.scale_factor_at_central_meridian,
225 )
226 crs = projection
227 else:
228 # Assume polar projection for regional grids encompassing N. Pole
229 if ymin > 20.0 and ymax > 80.0:
230 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
231 elif ymin < -80.0 and ymax < -20.0:
232 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
233 # Define regular map projection for non-rotated pole inputs.
234 # Alternatives might include e.g. for global model outputs:
235 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
236 # projection = ccrs.NearsidePerspective(
237 # central_longitude=180.0,
238 # central_latitude=0,
239 # satellite_height=35785831,
240 # )
241 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
242 else:
243 projection = ccrs.PlateCarree(central_longitude=central_longitude)
244 crs = ccrs.PlateCarree()
246 # Define axes for plot (or subplot) with required map projection.
247 if subplot is not None:
248 axes = figure.add_subplot(
249 grid_size[0], grid_size[1], subplot, projection=projection
250 )
251 else:
252 axes = figure.add_subplot(projection=projection)
254 # Add coastlines and borderlines if cube contains x and y map coordinates.
255 # Avoid adding lines for 2D masked data or specific fixed ancillary spatial plots.
256 if (cube.ndim > 1 and iris.util.is_masked(cube.data)) or any(
257 name in cube.name() for name in ["land_", "orography", "altitude"]
258 ):
259 pass
260 else:
261 if cmap.name in ["viridis", "Greys"]:
262 coastcol = "magenta"
263 else:
264 coastcol = "black"
265 logger.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
266 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
267 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
269 # Add gridlines.
270 gl = axes.gridlines(
271 alpha=0.3,
272 draw_labels=True,
273 dms=False,
274 x_inline=False,
275 y_inline=False,
276 )
277 gl.top_labels = False
278 gl.right_labels = False
279 if subplot:
280 gl.bottom_labels = False
281 gl.left_labels = False
282 if subplot % grid_size[1] == 1:
283 gl.left_labels = True
284 if subplot > ((grid_size[0] - 1) * grid_size[1]): 284 ↛ 289line 284 didn't jump to line 289 because the condition on line 284 was always true
285 gl.bottom_labels = True
287 # If is lat/lon spatial map, fix extent to keep plot tight.
288 # Specifying crs within set_extent helps ensure only data region is shown.
289 if isinstance(
290 coord_system, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
291 ):
292 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
294 except ValueError:
295 # Skip if not both x and y map coordinates.
296 axes = figure.gca()
298 return axes
301def _get_plot_resolution() -> int:
302 """Get resolution of rasterised plots in pixels per inch."""
303 return get_recipe_metadata().get("plot_resolution", 100)
306def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
307 """Return title and filename based on start and end points or bounds."""
308 if use_bounds and seq_coord.has_bounds():
309 vals = seq_coord.bounds.flatten()
310 else:
311 vals = seq_coord.points
312 start = seq_coord.units.title(vals[0])
313 end = seq_coord.units.title(vals[-1])
315 if start == end:
316 sequence_title = f"\n [{start}]"
317 sequence_fname = f"_{filename_slugify(start)}"
318 else:
319 sequence_title = f"\n [{start} to {end}]"
320 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
322 # Do not include time if coord set to zero.
323 if (
324 seq_coord.units == "hours since 0001-01-01 00:00:00"
325 and vals[0] == 0
326 and vals[-1] == 0
327 ):
328 sequence_title = ""
329 sequence_fname = ""
331 return sequence_title, sequence_fname
334def _set_title_and_filename(
335 seq_coord: iris.coords.Coord,
336 nplot: int,
337 recipe_title: str,
338 filename: str,
339):
340 """Set plot title and filename based on cube coordinate.
342 Parameters
343 ----------
344 sequence_coordinate: iris.coords.Coord
345 Coordinate about which to make a plot sequence.
346 nplot: int
347 Number of output plots to generate - controls title/naming.
348 recipe_title: str
349 Default plot title, potentially to update.
350 filename: str
351 Input plot filename, potentially to update.
353 Returns
354 -------
355 plot_title: str
356 Output formatted plot title string, based on plotted data.
357 plot_filename: str
358 Output formatted plot filename string.
359 """
360 ndim = seq_coord.ndim
361 npoints = np.size(seq_coord.points)
362 sequence_title = ""
363 sequence_fname = ""
365 # Case 1: Multiple dimension sequence input - list number of aggregated cases
366 # (e.g. aggregation histogram plots)
367 if ndim > 1:
368 ncase = np.shape(seq_coord)[0]
369 sequence_title = f"\n [{ncase} cases]"
370 sequence_fname = f"_{ncase}cases"
372 # Case 2: Single dimension input
373 else:
374 # Single sequence point
375 if npoints == 1:
376 if nplot > 1:
377 # Default labels for sequence inputs
378 sequence_value = seq_coord.units.title(seq_coord.points[0])
379 sequence_value = sequence_value.replace(" unknown", "")
380 sequence_title = f"\n [{sequence_value}]"
381 sequence_fname = f"_{filename_slugify(sequence_value)}"
382 else:
383 # Aggregated attribute available where input collapsed over aggregation
384 try:
385 ncase = seq_coord.attributes["number_reference_times"]
386 sequence_title = f"\n [{ncase} cases]"
387 sequence_fname = f"_{ncase}cases"
388 except KeyError:
389 sequence_title, sequence_fname = _get_start_end_strings(
390 seq_coord, use_bounds=seq_coord.has_bounds()
391 )
392 # Multiple sequence (e.g. time) points
393 else:
394 sequence_title, sequence_fname = _get_start_end_strings(
395 seq_coord, use_bounds=False
396 )
398 # Set plot title and filename
399 plot_title = f"{recipe_title}{sequence_title}"
401 # Set plot filename, defaulting to user input if provided.
402 if filename is None:
403 filename = slugify(recipe_title)
404 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
405 else:
406 if nplot > 1:
407 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
408 else:
409 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
411 return plot_title, plot_filename
414def _select_series_coord(cube, series_coordinate):
415 """Determine the grid coordinates to use to calculate grid spacing."""
416 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
417 if series_coordinate in spacing_coordinates: 417 ↛ 423line 417 didn't jump to line 423 because the condition on line 417 was always true
418 # Try the requested coordinate first then the fallbacks in order.
419 fallbacks = [series_coordinate] + [
420 c for c in spacing_coordinates if c != series_coordinate
421 ]
422 else:
423 fallbacks = {series_coordinate}
425 # Try each possible coordinate.
426 for coord in fallbacks:
427 try:
428 return cube.coord(coord)
429 except iris.exceptions.CoordinateNotFoundError:
430 logger.debug("Coordinate %s not found.", coord)
432 # If we get here, none of the fallback options were found.
433 raise iris.exceptions.CoordinateNotFoundError(
434 f"No valid coordinate found for '{series_coordinate}' "
435 f"or fallback options {fallbacks}"
436 )
439def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
440 """Control postage stamp plot output titles based on stamp coordinate."""
441 if stamp_coord.name() == "realization":
442 mtitle = "Member"
443 else:
444 mtitle = stamp_coord.name().capitalize()
446 if stamp_coord.name() == "time":
447 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
448 else:
449 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
451 return mtitle
454def _set_axis_range(cubes):
455 """Get minimum and maximum from levels information."""
456 levels = None
457 for cube in cubes: 457 ↛ 473line 457 didn't jump to line 473 because the loop on line 457 didn't complete
458 # First check if user-specified "auto" range variable.
459 # This maintains the value of levels as None, so proceed.
460 _, levels, _ = colorbar_map_levels(cube, axis="y")
461 if levels is None:
462 break
463 # If levels is changed, recheck to use the vmin,vmax or
464 # levels-based ranges for histogram plots.
465 _, levels, _ = colorbar_map_levels(cube)
466 logger.debug("levels: %s", levels)
467 if levels is not None: 467 ↛ 457line 467 didn't jump to line 457 because the condition on line 467 was always true
468 vmin = min(levels)
469 vmax = max(levels)
470 logger.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
471 break
473 if levels is None:
474 vmin = min(cb.data.min() for cb in cubes)
475 vmax = max(cb.data.max() for cb in cubes)
477 return vmin, vmax
480def _find_matched_slices(cubes, sequence_coordinate):
481 """Identify matched cubes in CubeList by sequence_coordinate values.
483 Ensures common points are compared for multiple cube inputs.
484 """
485 all_points = sorted(
486 set(
487 itertools.chain.from_iterable(
488 cb.coord(sequence_coordinate).points for cb in cubes
489 )
490 )
491 )
492 all_slices = list(
493 itertools.chain.from_iterable(
494 cb.slices_over(sequence_coordinate) for cb in cubes
495 )
496 )
497 # Matched slices (matched by seq coord point; it may happen that
498 # evaluated models do not cover the same seq coord range, hence matching
499 # necessary)
500 cube_iterables = [
501 iris.cube.CubeList(
502 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
503 )
504 for point in all_points
505 ]
507 return cube_iterables
510def _plot_and_save_spatial_plot(
511 cube: iris.cube.Cube,
512 filename: str,
513 title: str,
514 method: Literal["contourf", "pcolormesh", "scatter"],
515 overlay_cube: iris.cube.Cube | None = None,
516 contour_cube: iris.cube.Cube | None = None,
517 point_cube: iris.cube.Cube | None = None,
518 **kwargs,
519):
520 """Plot and save a spatial plot.
522 Parameters
523 ----------
524 cube: Cube
525 2 dimensional (lat and lon) Cube of the data to plot.
526 filename: str
527 Filename of the plot to write.
528 title: str
529 Plot title.
530 method: "contourf" | "pcolormesh" | "scatter"
531 The plotting method to use
532 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
533 overlay_cube: Cube, optional
534 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
535 contour_cube: Cube, optional
536 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
537 point_cube: Cube, optional
538 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
539 """
540 # Setup plot details, size, resolution, etc.
541 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
543 # Specify the color bar
544 cmap, levels, norm = colorbar_map_levels(cube)
546 # If overplotting, set required colorbars
547 if overlay_cube:
548 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
549 if contour_cube:
550 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
552 # Setup plot map projection, extent and coastlines and borderlines.
553 axes = _setup_spatial_map(cube, fig, cmap)
555 # Set colorscale bounds
556 try:
557 vmin = min(levels)
558 vmax = max(levels)
559 except TypeError:
560 vmin, vmax = None, None
561 # Ensure to use norm and not vmin/vmax if levels are defined.
562 if norm is not None:
563 vmin = None
564 vmax = None
565 logger.debug("Plotting using defined levels.")
567 # Plot the field.
568 if method == "contourf":
569 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
570 elif method == "pcolormesh":
571 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
572 elif method == "scatter":
573 # Scatter plot of the field. The marker size is chosen to give
574 # symbols that decrease in size as the number of data points
575 # increases, although the fraction of the figure covered by
576 # symbols increases roughly as N^(1/2), disregarding overlaps,
577 # and has been selected for the default figure size of (10, 10).
578 # Should this be changed, the marker size should be adjusted in
579 # proportion to the area of the figure.
580 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
581 lat_axis, lon_axis = get_cube_yxcoordname(cube)
582 plot = iplt.scatter(
583 cube.coord(lon_axis),
584 cube.coord(lat_axis),
585 c=cube.data[:],
586 s=mrk_size,
587 cmap=cmap,
588 edgecolors="k",
589 norm=norm,
590 vmin=vmin,
591 vmax=vmax,
592 )
593 else:
594 raise ValueError(f"Unknown plotting method: {method}")
596 # Overplot overlay field, if required
597 if overlay_cube:
598 try:
599 over_vmin = min(over_levels)
600 over_vmax = max(over_levels)
601 except TypeError:
602 over_vmin, over_vmax = None, None
603 if over_norm is not None: 603 ↛ 604line 603 didn't jump to line 604 because the condition on line 603 was never true
604 over_vmin = None
605 over_vmax = None
606 overlay = iplt.pcolormesh(
607 overlay_cube,
608 cmap=over_cmap,
609 norm=over_norm,
610 alpha=0.8,
611 vmin=over_vmin,
612 vmax=over_vmax,
613 )
614 # Overplot contour field, if required, with contour labelling.
615 if contour_cube:
616 contour = iplt.contour(
617 contour_cube,
618 colors="darkgray",
619 levels=cntr_levels,
620 norm=cntr_norm,
621 alpha=0.5,
622 linestyles="--",
623 linewidths=1,
624 )
625 plt.clabel(contour)
626 # Overplot valid elements of point-based field, if required.
627 # Check for non-masked points only to avoid plotting missing data.
628 if point_cube:
629 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
630 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
631 lon_coord = point_cube.coord(lon_axis)
632 lat_coord = point_cube.coord(lat_axis)
633 valid = ~point_cube.data.mask
634 valid_lon = iris.coords.AuxCoord(
635 lon_coord.points[valid],
636 standard_name=lon_coord.standard_name,
637 units=lon_coord.units,
638 coord_system=lon_coord.coord_system,
639 )
640 valid_lat = iris.coords.AuxCoord(
641 lat_coord.points[valid],
642 standard_name=lat_coord.standard_name,
643 units=lat_coord.units,
644 coord_system=lat_coord.coord_system,
645 )
646 iplt.scatter(
647 valid_lon,
648 valid_lat,
649 c=point_cube.data[valid],
650 s=mrk_size,
651 cmap=cmap,
652 edgecolors="k",
653 norm=norm,
654 vmin=vmin,
655 vmax=vmax,
656 )
658 # Check to see if transect, and if so, adjust y axis.
659 if is_transect(cube):
660 if "pressure" in [coord.name() for coord in cube.coords()]:
661 axes.invert_yaxis()
662 axes.set_yscale("log")
663 axes.set_ylim(1100, 100)
664 # If both model_level_number and level_height exists, iplt can construct
665 # plot as a function of height above orography (NOT sea level).
666 elif {"model_level_number", "level_height"}.issubset( 666 ↛ 671line 666 didn't jump to line 671 because the condition on line 666 was always true
667 {coord.name() for coord in cube.coords()}
668 ):
669 axes.set_yscale("log")
671 axes.set_title(
672 f"{title}\n"
673 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
674 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
675 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
676 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
677 fontsize=16,
678 )
680 # Inset code
681 axins = inset_axes(
682 axes,
683 width="20%",
684 height="20%",
685 loc="upper right",
686 axes_class=GeoAxes,
687 axes_kwargs={"map_projection": ccrs.PlateCarree()},
688 )
690 # Slightly transparent to reduce plot blocking.
691 axins.patch.set_alpha(0.4)
693 axins.coastlines(resolution="50m")
694 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
696 SLat, SLon, ELat, ELon = (
697 float(coord) for coord in cube.attributes["transect_coords"].split("_")
698 )
700 # Draw line between them
701 axins.plot(
702 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
703 )
705 # Plot points (note: lon, lat order for Cartopy)
706 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
707 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
709 lon_min, lon_max = sorted([SLon, ELon])
710 lat_min, lat_max = sorted([SLat, ELat])
712 # Midpoints
713 lon_mid = (lon_min + lon_max) / 2
714 lat_mid = (lat_min + lat_max) / 2
716 # Maximum half-range
717 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
718 if half_range == 0: # points identical → provide small default 718 ↛ 722line 718 didn't jump to line 722 because the condition on line 718 was always true
719 half_range = 1
721 # Set square extent
722 axins.set_extent(
723 [
724 lon_mid - half_range,
725 lon_mid + half_range,
726 lat_mid - half_range,
727 lat_mid + half_range,
728 ],
729 crs=ccrs.PlateCarree(),
730 )
732 # Ensure square aspect
733 axins.set_aspect("equal")
735 else:
736 # Add title.
737 axes.set_title(title, fontsize=16)
739 # Adjust padding if spatial plot or transect
740 if is_transect(cube):
741 yinfopad = -0.1
742 ycbarpad = 0.1
743 else:
744 yinfopad = 0.01
745 ycbarpad = 0.042
747 # Add watermark with min/max/mean. Currently not user togglable.
748 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
749 axes.annotate(
750 f"Min: {np.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(cube.data):.3g}",
751 xy=(0.025, yinfopad),
752 xycoords="axes fraction",
753 xytext=(-5, 5),
754 textcoords="offset points",
755 ha="left",
756 va="bottom",
757 size=11,
758 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
759 )
761 # Add secondary colour bar for overlay_cube field if required.
762 if overlay_cube:
763 cbarB = fig.colorbar(
764 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
765 )
766 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
767 # add ticks and tick_labels for every levels if less than 20 levels exist
768 if over_levels is not None and len(over_levels) < 20: 768 ↛ 769line 768 didn't jump to line 769 because the condition on line 768 was never true
769 cbarB.set_ticks(over_levels)
770 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
771 if any(
772 var in overlay_cube.name()
773 for var in ("rainfall", "snowfall", "visibility")
774 ):
775 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
776 logger.debug("Set secondary colorbar ticks and labels.")
778 # Add main colour bar.
779 cbar = fig.colorbar(
780 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
781 )
783 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
784 # add ticks and tick_labels for every levels if less than 20 levels exist
785 if levels is not None and len(levels) < 20:
786 cbar.set_ticks(levels)
787 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
788 if any(var in cube.name() for var in ("rainfall", "snowfall", "visibility")): 788 ↛ 791line 788 didn't jump to line 791 because the condition on line 788 was always true
789 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
790 # Tick labels for rainfall rates from Nimrod radar data.
791 if "rainfall rate composite" in cube.name(): 791 ↛ 792line 791 didn't jump to line 792 because the condition on line 791 was never true
792 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
793 # Tick labels for rain accumulations from Nimrod radar data.
794 if "rain accumulation" in cube.name(): 794 ↛ 795line 794 didn't jump to line 795 because the condition on line 794 was never true
795 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
796 if "wts accumulation" in cube.name(): 796 ↛ 797line 796 didn't jump to line 797 because the condition on line 796 was never true
797 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
798 cbar.minorticks_off()
799 cbar.set_ticks(tick_levels)
800 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
801 cbar.set_label(label=f"{cube.name()}", size=14)
802 # Tick labels for model rainfall data.
803 if "surface_microphysical" in cube.name(): 803 ↛ 806line 803 didn't jump to line 806 because the condition on line 803 was always true
804 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
805 # Tick labels for Nimrod weights data.
806 logger.debug("Set colorbar ticks and labels.")
808 # Save plot.
809 _save_close_figure(fig, "spatial", filename)
812def _plot_and_save_postage_stamp_spatial_plot(
813 cube: iris.cube.Cube,
814 filename: str,
815 stamp_coordinate: str,
816 title: str,
817 method: Literal["contourf", "pcolormesh"],
818 overlay_cube: iris.cube.Cube | None = None,
819 contour_cube: iris.cube.Cube | None = None,
820 **kwargs,
821):
822 """Plot postage stamp spatial plots from an ensemble.
824 Parameters
825 ----------
826 cube: Cube
827 Iris cube of data to be plotted. It must have the stamp coordinate.
828 filename: str
829 Filename of the plot to write.
830 stamp_coordinate: str
831 Coordinate that becomes different plots.
832 method: "contourf" | "pcolormesh"
833 The plotting method to use.
834 overlay_cube: Cube, optional
835 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
836 contour_cube: Cube, optional
837 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
839 Raises
840 ------
841 ValueError
842 If the cube doesn't have the right dimensions.
843 """
844 # Use the smallest square grid that will fit the members.
845 nmember = len(cube.coord(stamp_coordinate).points)
846 grid_rows = int(math.sqrt(nmember))
847 grid_size = math.ceil(nmember / grid_rows)
849 fig = plt.figure(
850 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
851 )
853 # Specify the color bar
854 cmap, levels, norm = colorbar_map_levels(cube)
855 # If overplotting, set required colorbars
856 if overlay_cube: 856 ↛ 857line 856 didn't jump to line 857 because the condition on line 856 was never true
857 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
858 if contour_cube: 858 ↛ 859line 858 didn't jump to line 859 because the condition on line 858 was never true
859 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
861 # Make a subplot for each member.
862 for member, subplot in zip(
863 cube.slices_over(stamp_coordinate),
864 range(1, grid_size * grid_rows + 1),
865 strict=False,
866 ):
867 # Setup subplot map projection, extent and coastlines and borderlines.
868 axes = _setup_spatial_map(
869 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
870 )
871 if method == "contourf":
872 # Filled contour plot of the field.
873 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
874 elif method == "pcolormesh":
875 if levels is not None:
876 vmin = min(levels)
877 vmax = max(levels)
878 else:
879 raise TypeError("Unknown vmin and vmax range.")
880 vmin, vmax = None, None
881 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
882 # if levels are defined.
883 if norm is not None: 883 ↛ 884line 883 didn't jump to line 884 because the condition on line 883 was never true
884 vmin = None
885 vmax = None
886 # pcolormesh plot of the field.
887 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
888 else:
889 raise ValueError(f"Unknown plotting method: {method}")
891 # Overplot overlay field, if required
892 if overlay_cube: 892 ↛ 893line 892 didn't jump to line 893 because the condition on line 892 was never true
893 try:
894 over_vmin = min(over_levels)
895 over_vmax = max(over_levels)
896 except TypeError:
897 over_vmin, over_vmax = None, None
898 if over_norm is not None:
899 over_vmin = None
900 over_vmax = None
901 iplt.pcolormesh(
902 overlay_cube[member.coord(stamp_coordinate).points[0]],
903 cmap=over_cmap,
904 norm=over_norm,
905 alpha=0.6,
906 vmin=over_vmin,
907 vmax=over_vmax,
908 )
909 # Overplot contour field, if required
910 if contour_cube: 910 ↛ 911line 910 didn't jump to line 911 because the condition on line 910 was never true
911 iplt.contour(
912 contour_cube[member.coord(stamp_coordinate).points[0]],
913 colors="darkgray",
914 levels=cntr_levels,
915 norm=cntr_norm,
916 alpha=0.6,
917 linestyles="--",
918 linewidths=1,
919 )
920 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
921 axes.set_title(f"{mtitle}")
923 # Put the shared colorbar in its own axes.
924 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
925 colorbar = fig.colorbar(
926 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
927 )
928 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
930 # Overall figure title.
931 fig.suptitle(title, fontsize=16)
933 # Save plot.
934 _save_close_figure(fig, "contour postate stamp", filename)
937def _plot_and_save_line_series(
938 cubes: iris.cube.CubeList,
939 coords: list[iris.coords.Coord],
940 ensemble_coord: str,
941 filename: str,
942 title: str,
943 **kwargs,
944):
945 """Plot and save a 1D line series.
947 Parameters
948 ----------
949 cubes: Cube or CubeList
950 Cube or CubeList containing the cubes to plot on the y-axis.
951 coords: list[Coord]
952 Coordinates to plot on the x-axis, one per cube.
953 ensemble_coord: str
954 Ensemble coordinate in the cube.
955 filename: str
956 Filename of the plot to write.
957 title: str
958 Plot title.
959 """
960 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
962 model_colors_map = get_model_colors_map(cubes)
964 # Store min/max ranges.
965 y_levels = []
967 # Check match-up across sequence coords gives consistent sizes
968 validate_cubes_coords(cubes, coords)
970 for cube, coord in zip(cubes, coords, strict=True):
971 label = None
972 color = "black"
973 if model_colors_map:
974 label = cube.attributes.get("model_name")
975 color = model_colors_map.get(label)
976 if not cube.coords(ensemble_coord): 976 ↛ 978line 976 didn't jump to line 978 because the condition on line 976 was never true
977 # No ensemble coordinate — plot the cube directly as a single line.
978 iplt.plot(coord, cube, color=color, marker="o", ls="-", lw=3, label=label)
979 else:
980 for cube_slice in cube.slices_over(ensemble_coord):
981 # Label with (control) if part of an ensemble or not otherwise.
982 if cube_slice.coord(ensemble_coord).points == [0]:
983 iplt.plot(
984 coord,
985 cube_slice,
986 color=color,
987 marker="o",
988 ls="-",
989 lw=3,
990 label=f"{label} (control)"
991 if len(cube.coord(ensemble_coord).points) > 1
992 else label,
993 )
994 # Label with (perturbed) if part of an ensemble and not the control.
995 else:
996 iplt.plot(
997 coord,
998 cube_slice,
999 color=color,
1000 ls="-",
1001 lw=1.5,
1002 alpha=0.75,
1003 label=f"{label} (member)",
1004 )
1006 # Calculate the global min/max if multiple cubes are given.
1007 _, levels, _ = colorbar_map_levels(cube, axis="y")
1008 if levels is not None: 1008 ↛ 1009line 1008 didn't jump to line 1009 because the condition on line 1008 was never true
1009 y_levels.append(min(levels))
1010 y_levels.append(max(levels))
1012 # Get the current axes.
1013 ax = plt.gca()
1015 # Add some labels and tweak the style.
1016 # check if cubes[0] works for single cube if not CubeList
1017 if coords[0].name() == "time":
1018 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1019 else:
1020 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1021 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1022 ax.set_title(title, fontsize=16)
1024 ax.ticklabel_format(axis="y", useOffset=False)
1025 ax.tick_params(axis="x", labelrotation=15)
1026 ax.tick_params(axis="both", labelsize=12)
1028 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1029 if y_levels: 1029 ↛ 1030line 1029 didn't jump to line 1030 because the condition on line 1029 was never true
1030 ax.set_ylim(min(y_levels), max(y_levels))
1031 logger.debug("Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels))
1032 else:
1033 ax.autoscale()
1035 # Add gridlines
1036 ax.grid(linestyle="--", color="grey", linewidth=1)
1037 # Add zero line
1038 ymin, ymax = ax.get_ylim()
1039 if ymin < 0.0 and ymax > 0.0:
1040 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1041 # Identify unique labels for legend
1042 handles = list(
1043 {
1044 label: handle
1045 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1046 }.values()
1047 )
1048 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1050 # Save plot.
1051 _save_close_figure(fig, "line", filename)
1054def _plot_and_save_line_power_spectrum_series(
1055 cubes: iris.cube.Cube | iris.cube.CubeList,
1056 coords: list[iris.coords.Coord],
1057 ensemble_coord: str,
1058 filename: str,
1059 title: str,
1060 series_coordinate: str,
1061 **kwargs,
1062):
1063 """Plot and save a 1D line series.
1065 Parameters
1066 ----------
1067 cubes: Cube or CubeList
1068 Cube or CubeList containing the cubes to plot on the y-axis.
1069 coords: list[Coord]
1070 Coordinates to plot on the x-axis, one per cube.
1071 ensemble_coord: str
1072 Ensemble coordinate in the cube.
1073 filename: str
1074 Filename of the plot to write.
1075 title: str
1076 Plot title.
1077 series_coordinate: str
1078 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1079 """
1080 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1081 model_colors_map = get_model_colors_map(cubes)
1082 ax = plt.gca()
1084 # Store min/max ranges.
1085 y_levels = []
1087 line_marker = None
1088 line_width = 1
1090 for cube in iter_maybe(cubes):
1091 # next 2 lines replace chunk of code.
1092 xcoord = _select_series_coord(cube, series_coordinate)
1093 xname = xcoord.points
1095 yfield = cube.data # power spectrum
1096 label = None
1097 color = "black"
1098 if model_colors_map: 1098 ↛ 1101line 1098 didn't jump to line 1101 because the condition on line 1098 was always true
1099 label = cube.attributes.get("model_name")
1100 color = model_colors_map.get(label)
1101 for cube_slice in cube.slices_over(ensemble_coord):
1102 # Label with (control) if part of an ensemble or not otherwise.
1103 if cube_slice.coord(ensemble_coord).points == [0]: 1103 ↛ 1117line 1103 didn't jump to line 1117 because the condition on line 1103 was always true
1104 ax.plot(
1105 xname,
1106 yfield,
1107 color=color,
1108 marker=line_marker,
1109 ls="-",
1110 lw=line_width,
1111 label=f"{label} (control)"
1112 if len(cube.coord(ensemble_coord).points) > 1
1113 else label,
1114 )
1115 # Label with (perturbed) if part of an ensemble and not the control.
1116 else:
1117 ax.plot(
1118 xname,
1119 yfield,
1120 color=color,
1121 ls="-",
1122 lw=1.5,
1123 alpha=0.75,
1124 label=f"{label} (member)",
1125 )
1127 # Calculate the global min/max if multiple cubes are given.
1128 _, levels, _ = colorbar_map_levels(cube, axis="y")
1129 if levels is not None: 1129 ↛ 1130line 1129 didn't jump to line 1130 because the condition on line 1129 was never true
1130 y_levels.append(min(levels))
1131 y_levels.append(max(levels))
1133 # Add some labels and tweak the style.
1135 title = f"{title}"
1136 ax.set_title(title, fontsize=16)
1138 # Set appropriate x-axis label based on coordinate
1139 if series_coordinate == "wavelength" or ( 1139 ↛ 1142line 1139 didn't jump to line 1142 because the condition on line 1139 was never true
1140 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1141 ):
1142 ax.set_xlabel("Wavelength (km)", fontsize=14)
1143 elif series_coordinate == "physical_wavenumber" or ( 1143 ↛ 1146line 1143 didn't jump to line 1146 because the condition on line 1143 was never true
1144 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1145 ):
1146 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1147 else: # frequency or check units
1148 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1148 ↛ 1149line 1148 didn't jump to line 1149 because the condition on line 1148 was never true
1149 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1150 else:
1151 ax.set_xlabel("Wavenumber", fontsize=14)
1153 ax.set_ylabel("Power Spectral Density", fontsize=14)
1154 ax.tick_params(axis="both", labelsize=12)
1156 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1158 # Set log-log scale
1159 ax.set_xscale("log")
1160 ax.set_yscale("log")
1162 # Add gridlines
1163 ax.grid(linestyle="--", color="grey", linewidth=1)
1164 # Ientify unique labels for legend
1165 handles = list(
1166 {
1167 label: handle
1168 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1169 }.values()
1170 )
1171 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1173 # Save plot.
1174 _save_close_figure(fig, "line power spectrum", filename)
1177def _plot_and_save_vertical_line_series(
1178 cubes: iris.cube.CubeList,
1179 coords: list[iris.coords.Coord],
1180 ensemble_coord: str,
1181 filename: str,
1182 series_coordinate: str,
1183 title: str,
1184 vmin: float,
1185 vmax: float,
1186 **kwargs,
1187):
1188 """Plot and save a 1D line series in vertical.
1190 Parameters
1191 ----------
1192 cubes: CubeList
1193 1 dimensional Cube or CubeList of the data to plot on x-axis.
1194 coord: list[Coord]
1195 Coordinates to plot on the y-axis, one per cube.
1196 ensemble_coord: str
1197 Ensemble coordinate in the cube.
1198 filename: str
1199 Filename of the plot to write.
1200 series_coordinate: str
1201 Coordinate to use as vertical axis.
1202 title: str
1203 Plot title.
1204 vmin: float
1205 Minimum value for the x-axis.
1206 vmax: float
1207 Maximum value for the x-axis.
1208 """
1209 # plot the vertical pressure axis using log scale
1210 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1212 model_colors_map = get_model_colors_map(cubes)
1214 # Check match-up across sequence coords gives consistent sizes
1215 validate_cubes_coords(cubes, coords)
1217 for cube, coord in zip(cubes, coords, strict=True):
1218 label = None
1219 color = "black"
1220 if model_colors_map: 1220 ↛ 1221line 1220 didn't jump to line 1221 because the condition on line 1220 was never true
1221 label = cube.attributes.get("model_name")
1222 color = model_colors_map.get(label)
1224 for cube_slice in cube.slices_over(ensemble_coord):
1225 # If ensemble data given plot control member with (control)
1226 # unless single forecast.
1227 if cube_slice.coord(ensemble_coord).points == [0]:
1228 iplt.plot(
1229 cube_slice,
1230 coord,
1231 color=color,
1232 marker="o",
1233 ls="-",
1234 lw=3,
1235 label=f"{label} (control)"
1236 if len(cube.coord(ensemble_coord).points) > 1
1237 else label,
1238 )
1239 # If ensemble data given plot perturbed members with (perturbed).
1240 else:
1241 iplt.plot(
1242 cube_slice,
1243 coord,
1244 color=color,
1245 ls="-",
1246 lw=1.5,
1247 alpha=0.75,
1248 label=f"{label} (member)",
1249 )
1251 # Get the current axis
1252 ax = plt.gca()
1254 # Special handling for pressure level data.
1255 if series_coordinate == "pressure": 1255 ↛ 1277line 1255 didn't jump to line 1277 because the condition on line 1255 was always true
1256 # Invert y-axis and set to log scale.
1257 ax.invert_yaxis()
1258 ax.set_yscale("log")
1260 # Define y-ticks and labels for pressure log axis.
1261 y_tick_labels = [
1262 "1000",
1263 "850",
1264 "700",
1265 "500",
1266 "300",
1267 "200",
1268 "100",
1269 ]
1270 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1272 # Set y-axis limits and ticks.
1273 ax.set_ylim(1100, 100)
1275 # Test if series_coordinate is model level data. The UM data uses
1276 # model_level_number and lfric uses full_levels as coordinate.
1277 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1278 # Define y-ticks and labels for vertical axis.
1279 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1280 y_tick_labels = [str(int(i)) for i in y_ticks]
1281 ax.set_ylim(min(y_ticks), max(y_ticks))
1283 ax.set_yticks(y_ticks)
1284 ax.set_yticklabels(y_tick_labels)
1286 # Set x-axis limits.
1287 ax.set_xlim(vmin, vmax)
1288 # Mark y=0 if present in plot.
1289 if vmin < 0.0 and vmax > 0.0: 1289 ↛ 1290line 1289 didn't jump to line 1290 because the condition on line 1289 was never true
1290 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1292 # Add some labels and tweak the style.
1293 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1294 ax.set_xlabel(
1295 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1296 )
1297 ax.set_title(title, fontsize=16)
1298 ax.ticklabel_format(axis="x")
1299 ax.tick_params(axis="y")
1300 ax.tick_params(axis="both", labelsize=12)
1302 # Add gridlines
1303 ax.grid(linestyle="--", color="grey", linewidth=1)
1304 # Ientify unique labels for legend
1305 handles = list(
1306 {
1307 label: handle
1308 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1309 }.values()
1310 )
1311 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1313 # Save plot.
1314 _save_close_figure(fig, "vertical line", filename)
1317def _plot_and_save_scatter_plot(
1318 cube_x: iris.cube.Cube | iris.cube.CubeList,
1319 cube_y: iris.cube.Cube | iris.cube.CubeList,
1320 filename: str,
1321 title: str,
1322 one_to_one: bool,
1323 model_names: list[str] | None = None,
1324 **kwargs,
1325):
1326 """Plot and save a 2D scatter plot.
1328 Parameters
1329 ----------
1330 cube_x: Cube | CubeList
1331 1 dimensional Cube or CubeList of the data to plot on x-axis.
1332 cube_y: Cube | CubeList
1333 1 dimensional Cube or CubeList of the data to plot on y-axis.
1334 filename: str
1335 Filename of the plot to write.
1336 title: str
1337 Plot title.
1338 one_to_one: bool
1339 Whether a 1:1 line is plotted.
1340 """
1341 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1342 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1343 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1344 # over the pairs simultaneously.
1346 # Ensure cube_x and cube_y are iterable
1347 cube_x_iterable = iter_maybe(cube_x)
1348 cube_y_iterable = iter_maybe(cube_y)
1350 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1351 iplt.scatter(cube_x_iter, cube_y_iter)
1352 if one_to_one is True:
1353 plt.plot(
1354 [
1355 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1356 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1357 ],
1358 [
1359 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1360 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1361 ],
1362 "k",
1363 linestyle="--",
1364 )
1365 ax = plt.gca()
1367 # Add some labels and tweak the style.
1368 if model_names is None:
1369 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1370 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1371 else:
1372 # Add the model names, these should be order of base (x) and other (y).
1373 ax.set_xlabel(
1374 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1375 )
1376 ax.set_ylabel(
1377 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1378 )
1379 ax.set_title(title, fontsize=16)
1380 ax.ticklabel_format(axis="y", useOffset=False)
1381 ax.tick_params(axis="x", labelrotation=15)
1382 ax.tick_params(axis="both", labelsize=12)
1383 ax.autoscale()
1385 # Save plot.
1386 _save_close_figure(fig, "scatter", filename)
1389def _plot_and_save_vector_plot(
1390 cube_u: iris.cube.Cube,
1391 cube_v: iris.cube.Cube,
1392 filename: str,
1393 title: str,
1394 method: Literal["contourf", "pcolormesh"],
1395 **kwargs,
1396):
1397 """Plot and save a 2D vector plot.
1399 Parameters
1400 ----------
1401 cube_u: Cube
1402 2 dimensional Cube of u component of the data.
1403 cube_v: Cube
1404 2 dimensional Cube of v component of the data.
1405 filename: str
1406 Filename of the plot to write.
1407 title: str
1408 Plot title.
1409 """
1410 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1411 # Create a cube containing the magnitude of the vector field.
1412 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1413 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1414 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1415 cube_vec_mag.rename(
1416 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1417 )
1419 # Specify the color bar
1420 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1422 # Setup plot map projection, extent and coastlines and borderlines.
1423 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1425 if method == "contourf":
1426 # Filled contour plot of the field.
1427 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1428 elif method == "pcolormesh":
1429 try:
1430 vmin = min(levels)
1431 vmax = max(levels)
1432 except TypeError:
1433 vmin, vmax = None, None
1434 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1435 # if levels are defined.
1436 if norm is not None:
1437 vmin = None
1438 vmax = None
1439 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1440 else:
1441 raise ValueError(f"Unknown plotting method: {method}")
1443 # Check to see if transect, and if so, adjust y axis.
1444 if is_transect(cube_vec_mag):
1445 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1446 axes.invert_yaxis()
1447 axes.set_yscale("log")
1448 axes.set_ylim(1100, 100)
1449 # If both model_level_number and level_height exists, iplt can construct
1450 # plot as a function of height above orography (NOT sea level).
1451 elif {"model_level_number", "level_height"}.issubset(
1452 {coord.name() for coord in cube_vec_mag.coords()}
1453 ):
1454 axes.set_yscale("log")
1456 axes.set_title(
1457 f"{title}\n"
1458 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1459 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1460 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1461 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1462 fontsize=16,
1463 )
1465 else:
1466 # Add title.
1467 axes.set_title(title, fontsize=16)
1469 # Add watermark with min/max/mean. Currently not user togglable.
1470 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1471 axes.annotate(
1472 f"Min: {np.nanmin(cube_vec_mag.data):.3g} Max: {np.nanmax(cube_vec_mag.data):.3g} Mean: {np.nanmean(cube_vec_mag.data):.3g}",
1473 xy=(0.05, -0.05),
1474 xycoords="axes fraction",
1475 xytext=(-5, 5),
1476 textcoords="offset points",
1477 ha="right",
1478 va="bottom",
1479 size=11,
1480 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1481 )
1483 # Add colour bar.
1484 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1485 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1486 # add ticks and tick_labels for every levels if less than 20 levels exist
1487 if levels is not None and len(levels) < 20:
1488 cbar.set_ticks(levels)
1489 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1491 # 30 barbs along the longest axis of the plot, or a barb per point for data
1492 # with less than 30 points.
1493 step = max(max(cube_u.shape) // 30, 1)
1494 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1496 # Save plot.
1497 _save_close_figure(fig, "vector", filename)
1500def _plot_and_save_histogram_series(
1501 cubes: iris.cube.Cube | iris.cube.CubeList,
1502 filename: str,
1503 title: str,
1504 vmin: float,
1505 vmax: float,
1506 **kwargs,
1507):
1508 """Plot and save a histogram series.
1510 Parameters
1511 ----------
1512 cubes: Cube or CubeList
1513 2 dimensional Cube or CubeList of the data to plot as histogram.
1514 filename: str
1515 Filename of the plot to write.
1516 title: str
1517 Plot title.
1518 vmin: float
1519 minimum for colorbar
1520 vmax: float
1521 maximum for colorbar
1522 """
1523 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1524 ax = plt.gca()
1526 model_colors_map = get_model_colors_map(cubes)
1528 # Set default that histograms will produce probability density function
1529 # at each bin (integral over range sums to 1).
1530 density = True
1532 for cube in iter_maybe(cubes):
1533 # Easier to check title (where var name originates)
1534 # than seeing if long names exist etc.
1535 # Exception case, where distribution better fits log scales/bins.
1536 if (
1537 ("surface_microphysical" in title)
1538 or ("rain accumulation" in title)
1539 or ("Rainfall rate Composite" in title)
1540 or ("Nimrod_5min" in title)
1541 ):
1542 if "amount" in title:
1543 # Compute histogram following Klingaman et al. (2017): ASoP
1544 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1545 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1546 density = False
1547 else:
1548 bins = 10.0 ** (
1549 np.arange(-10, 27, 1) / 10.0
1550 ) # Suggestion from RMED toolbox.
1551 bins = np.insert(bins, 0, 0)
1552 ax.set_yscale("log")
1553 vmin = bins[1]
1554 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1555 ax.set_xscale("log")
1556 elif "lightning" in title:
1557 bins = [0, 1, 2, 3, 4, 5]
1558 else:
1559 bins = np.linspace(vmin, vmax, 51)
1560 logger.debug(
1561 "Plotting histogram with %s bins %s - %s.",
1562 np.size(bins),
1563 np.min(bins),
1564 np.max(bins),
1565 )
1567 # Reshape cube data into a single array to allow for a single histogram.
1568 # Otherwise we plot xdim histograms stacked.
1569 cube_data_1d = (cube.data).flatten()
1571 label = None
1572 color = "black"
1573 if model_colors_map:
1574 label = cube.attributes.get("model_name")
1575 color = model_colors_map[label]
1576 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1578 # Compute area under curve.
1579 if (
1580 ("surface_microphysical" in title and "amount" in title)
1581 or ("rain_accumulation" in title)
1582 or ("Rainfall rate Composite" in title)
1583 or ("Nimrod_5min" in title)
1584 ):
1585 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1586 x = x * bin_mean / x.sum()
1587 x = x[1:]
1588 y = y[1:]
1590 ax.plot(
1591 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1592 )
1594 # Add some labels and tweak the style.
1595 ax.set_title(title, fontsize=16)
1596 ax.set_xlabel(
1597 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1598 )
1599 ax.set_ylabel("Normalised probability density", fontsize=14)
1600 if (
1601 ("surface_microphysical" in title and "amount" in title)
1602 or ("rain accumulation" in title)
1603 or ("Nimrod_5min" in title)
1604 ):
1605 ax.set_ylabel(
1606 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1607 )
1608 try:
1609 ax.set_xlim(vmin, vmax)
1610 except ValueError:
1611 pass
1612 ax.tick_params(axis="both", labelsize=12)
1614 # Overlay grid-lines onto histogram plot.
1615 ax.grid(linestyle="--", color="grey", linewidth=1)
1616 if model_colors_map:
1617 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1619 # Save plot.
1620 _save_close_figure(fig, "histogram", filename)
1623def _plot_and_save_postage_stamp_histogram_series(
1624 cube: iris.cube.Cube,
1625 filename: str,
1626 title: str,
1627 stamp_coordinate: str,
1628 vmin: float,
1629 vmax: float,
1630 **kwargs,
1631):
1632 """Plot and save postage (ensemble members) stamps for a histogram series.
1634 Parameters
1635 ----------
1636 cube: Cube
1637 2 dimensional Cube of the data to plot as histogram.
1638 filename: str
1639 Filename of the plot to write.
1640 title: str
1641 Plot title.
1642 stamp_coordinate: str
1643 Coordinate that becomes different plots.
1644 vmin: float
1645 minimum for pdf x-axis
1646 vmax: float
1647 maximum for pdf x-axis
1648 """
1649 # Use the smallest square grid that will fit the members.
1650 nmember = len(cube.coord(stamp_coordinate).points)
1651 grid_rows = int(math.sqrt(nmember))
1652 grid_size = math.ceil(nmember / grid_rows)
1654 fig = plt.figure(
1655 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1656 )
1657 # Make a subplot for each member.
1658 for member, subplot in zip(
1659 cube.slices_over(stamp_coordinate),
1660 range(1, grid_size * grid_rows + 1),
1661 strict=False,
1662 ):
1663 # Implicit interface is much easier here, due to needing to have the
1664 # cartopy GeoAxes generated.
1665 plt.subplot(grid_rows, grid_size, subplot)
1666 # Reshape cube data into a single array to allow for a single histogram.
1667 # Otherwise we plot xdim histograms stacked.
1668 member_data_1d = (member.data).flatten()
1669 plt.hist(member_data_1d, density=True, stacked=True)
1670 axes = plt.gca()
1671 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1672 axes.set_title(f"{mtitle}")
1673 axes.set_xlim(vmin, vmax)
1675 # Overall figure title.
1676 fig.suptitle(title, fontsize=16)
1678 # Save plot.
1679 _save_close_figure(fig, "histogram postage stamp", filename)
1682def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1683 cube: iris.cube.Cube,
1684 filename: str,
1685 title: str,
1686 stamp_coordinate: str,
1687 vmin: float,
1688 vmax: float,
1689 **kwargs,
1690):
1691 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1692 ax.set_title(title, fontsize=16)
1693 ax.set_xlim(vmin, vmax)
1694 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1695 ax.set_ylabel("normalised probability density", fontsize=14)
1696 # Loop over all slices along the stamp_coordinate
1697 for member in cube.slices_over(stamp_coordinate):
1698 # Flatten the member data to 1D
1699 member_data_1d = member.data.flatten()
1700 # Plot the histogram using plt.hist
1701 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1702 plt.hist(
1703 member_data_1d,
1704 density=True,
1705 stacked=True,
1706 label=f"{mtitle}",
1707 )
1709 # Add a legend
1710 ax.legend(fontsize=16)
1712 # Save plot.
1713 _save_close_figure(fig, "histogram postage stamp", filename)
1716def _plot_and_save_scatter_series(
1717 cubes: iris.cube.Cube | iris.cube.CubeList,
1718 filename: str,
1719 title: str,
1720 vmin: float,
1721 vmax: float,
1722 hexbin: bool,
1723 **kwargs,
1724):
1725 """Plot and save a scatter plot series.
1727 Parameters
1728 ----------
1729 cubes: Cube or CubeList
1730 2 dimensional Cube or CubeList of the data to plot as scatter.
1731 filename: str
1732 Filename of the plot to write.
1733 title: str
1734 Plot title.
1735 vmin: float
1736 minimum for colorbar
1737 vmax: float
1738 maximum for colorbar
1739 hexbin: bool
1740 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1741 Else scatter of all points, with potential to overplot many comparisons on same plot.
1742 """
1743 if hexbin:
1744 # Check cubes using same functionality as the difference operator.
1745 if len(cubes) != 2:
1746 raise ValueError(
1747 "Cubes should contain exactly 2 cubes for hexbin plotting."
1748 )
1749 title = title.replace("scatter", "hexbin")
1750 filename = filename.replace("scatter", "hexbin")
1752 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1753 ax = plt.gca()
1755 model_colors_map = get_model_colors_map(cubes)
1757 percentiles = np.arange(0, 100, 5)
1758 percentiles[0] = 1
1759 percentiles[-1] = 99
1760 quantiles = iris.cube.CubeList()
1762 # Loop through all output cubes for both data points and overplotting quantiles.
1763 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1764 for plottype in ["points", "quantiles"]:
1765 nplot = 0
1766 for cube in iter_maybe(cubes):
1767 label = None
1768 color = "black"
1769 if model_colors_map: 1769 ↛ 1774line 1769 didn't jump to line 1774 because the condition on line 1769 was always true
1770 label = cube.attributes.get("model_name")
1771 color = model_colors_map[label]
1773 # Plot all data points
1774 if plottype == "points":
1775 if nplot > 0:
1776 if hexbin:
1777 hb = plt.hexbin(
1778 cubes[0].data.flatten(),
1779 cube.data.flatten(),
1780 alpha=0.3,
1781 gridsize=100,
1782 mincnt=1,
1783 )
1784 else:
1785 plt.scatter(
1786 cubes[0].data.flatten(),
1787 cube.data.flatten(),
1788 color=color,
1789 marker="+",
1790 label=None,
1791 alpha=0.3,
1792 )
1794 elif plottype == "quantiles": 1794 ↛ 1813line 1794 didn't jump to line 1813 because the condition on line 1794 was always true
1795 # Construct Q-Q plot
1796 quantiles.append(
1797 cube.collapsed(
1798 cube.coords(dim_coords=True),
1799 iris.analysis.PERCENTILE,
1800 percent=percentiles,
1801 )
1802 )
1803 if nplot > 0:
1804 iplt.scatter(
1805 quantiles[0],
1806 quantiles[-1],
1807 color=color,
1808 marker="o",
1809 label=label,
1810 edgecolors="black",
1811 )
1813 nplot = nplot + 1
1815 # Add some labels and tweak the style.
1816 ax.set_title(title, fontsize=16)
1817 ax.set_xlabel(
1818 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1819 )
1820 ax.set_ylabel(
1821 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1822 )
1823 ax.tick_params(axis="both", labelsize=12)
1824 ax.autoscale()
1826 # Set 1:1 line and equal axes if scatter plot of common cube names
1827 nameA = iter_maybe(cubes)[0].name()
1828 nameB = iter_maybe(cubes)[1].name()
1829 if any(part in nameB.split("_") for part in nameA.split("_")): 1829 ↛ 1840line 1829 didn't jump to line 1840 because the condition on line 1829 was always true
1830 lims = [
1831 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1832 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1833 ]
1834 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1835 ax.set_aspect("equal")
1836 ax.set_xlim(lims)
1837 ax.set_ylim(lims)
1839 # Overlay grid-lines onto scatter plot.
1840 ax.grid(linestyle="--", color="grey", linewidth=1)
1841 if model_colors_map: 1841 ↛ 1845line 1841 didn't jump to line 1845 because the condition on line 1841 was always true
1842 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1844 # Add colorbar if hexbin output
1845 if hexbin:
1846 cb = plt.colorbar(
1847 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1848 )
1849 cb.set_label("Number of data points", size=12)
1851 # Save plot.
1852 _save_close_figure(fig, "scatter", filename)
1855def _spatial_plot(
1856 method: Literal["contourf", "pcolormesh", "scatter"],
1857 cube: iris.cube.Cube,
1858 filename: str | None,
1859 sequence_coordinate: str,
1860 stamp_coordinate: str,
1861 overlay_cube: iris.cube.Cube | None = None,
1862 contour_cube: iris.cube.Cube | None = None,
1863 point_cube: iris.cube.Cube | None = None,
1864 **kwargs,
1865):
1866 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1868 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1869 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1870 is present then postage stamp plots will be produced.
1872 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1873 be overplotted on the same figure.
1875 Parameters
1876 ----------
1877 method: "contourf" | "pcolormesh" | "scatter"
1878 The plotting method to use.
1879 Select choice of "contourf" or "pcolormesh" for gridded data.
1880 Use "scatter" for point-based data.
1881 cube: Cube
1882 Iris cube of the data to plot. It should have two spatial dimensions,
1883 such as lat and lon, and may also have a another two dimension to be
1884 plotted sequentially and/or as postage stamp plots.
1885 filename: str | None
1886 Name of the plot to write, used as a prefix for plot sequences. If None
1887 uses the recipe name.
1888 sequence_coordinate: str
1889 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1890 This coordinate must exist in the cube.
1891 stamp_coordinate: str
1892 Coordinate about which to plot postage stamp plots. Defaults to
1893 ``"realization"``.
1894 overlay_cube: Cube | None, optional
1895 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1896 contour_cube: Cube | None, optional
1897 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1898 point_cube: Cube | None, optional
1899 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
1901 Raises
1902 ------
1903 ValueError
1904 If the cube doesn't have the right dimensions.
1905 TypeError
1906 If the cube isn't a single cube.
1907 """
1908 # Ensure we've got a single cube.
1909 cube = check_single_cube(cube)
1911 # Set title based on recipe metadata or use cube name
1912 recipe_title = get_recipe_metadata().get("title", cube.name())
1914 # Check if there is a valid stamp coordinate in cube dimensions.
1915 if stamp_coordinate == "realization": 1915 ↛ 1920line 1915 didn't jump to line 1920 because the condition on line 1915 was always true
1916 stamp_coordinate = check_stamp_coordinate(cube)
1918 # Make postage stamp plots if stamp_coordinate exists and has more than a
1919 # single point.
1920 plotting_func = _plot_and_save_spatial_plot
1921 try:
1922 if cube.coord(stamp_coordinate).shape[0] > 1:
1923 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1924 except iris.exceptions.CoordinateNotFoundError:
1925 pass
1927 # Produce a geographical scatter plot if the data have a
1928 # dimension called observation or model_obs_error
1929 if any(
1930 crd.var_name == "station"
1931 or crd.var_name == "Station_Name"
1932 or crd.var_name == "model_obs_error"
1933 for crd in cube.coords()
1934 ):
1935 plotting_func = _plot_and_save_spatial_plot
1936 method = "scatter"
1938 # Must have a sequence coordinate.
1939 try:
1940 cube.coord(sequence_coordinate)
1941 except iris.exceptions.CoordinateNotFoundError as err:
1942 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1944 # Create a plot for each value of the sequence coordinate.
1945 plot_index = []
1946 nplot = np.size(cube.coord(sequence_coordinate).points)
1948 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1949 # Set plot titles and filename
1950 seq_coord = cube_slice.coord(sequence_coordinate)
1951 plot_title, plot_filename = _set_title_and_filename(
1952 seq_coord, nplot, recipe_title, filename
1953 )
1955 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1956 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1957 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1958 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1960 # Do the actual plotting.
1961 plotting_func(
1962 cube_slice,
1963 filename=plot_filename,
1964 stamp_coordinate=stamp_coordinate,
1965 title=plot_title,
1966 method=method,
1967 overlay_cube=overlay_slice,
1968 contour_cube=contour_slice,
1969 point_cube=point_slice,
1970 **kwargs,
1971 )
1972 plot_index.append(plot_filename)
1974 # Add list of plots to plot metadata.
1975 complete_plot_index = _append_to_plot_index(plot_index)
1977 # Make a page to display the plots.
1978 _make_plot_html_page(complete_plot_index)
1981####################
1982# Public functions #
1983####################
1986def spatial_contour_plot(
1987 cube: iris.cube.Cube,
1988 filename: str | None = None,
1989 sequence_coordinate: str = "time",
1990 stamp_coordinate: str = "realization",
1991 **kwargs,
1992) -> iris.cube.Cube:
1993 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1995 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1996 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1997 is present then postage stamp plots will be produced.
1999 Parameters
2000 ----------
2001 cube: Cube
2002 Iris cube of the data to plot. It should have two spatial dimensions,
2003 such as lat and lon, and may also have a another two dimension to be
2004 plotted sequentially and/or as postage stamp plots.
2005 filename: str, optional
2006 Name of the plot to write, used as a prefix for plot sequences. Defaults
2007 to the recipe name.
2008 sequence_coordinate: str, optional
2009 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2010 This coordinate must exist in the cube.
2011 stamp_coordinate: str, optional
2012 Coordinate about which to plot postage stamp plots. Defaults to
2013 ``"realization"``.
2015 Returns
2016 -------
2017 Cube
2018 The original cube (so further operations can be applied).
2020 Raises
2021 ------
2022 ValueError
2023 If the cube doesn't have the right dimensions.
2024 TypeError
2025 If the cube isn't a single cube.
2026 """
2027 _spatial_plot(
2028 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2029 )
2030 return cube
2033def spatial_pcolormesh_plot(
2034 cube: iris.cube.Cube,
2035 filename: str | None = None,
2036 sequence_coordinate: str = "time",
2037 stamp_coordinate: str = "realization",
2038 **kwargs,
2039) -> iris.cube.Cube:
2040 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2042 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2043 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2044 is present then postage stamp plots will be produced.
2046 This function is significantly faster than ``spatial_contour_plot``,
2047 especially at high resolutions, and should be preferred unless contiguous
2048 contour areas are important.
2050 Parameters
2051 ----------
2052 cube: Cube
2053 Iris cube of the data to plot. It should have two spatial dimensions,
2054 such as lat and lon, and may also have a another two dimension to be
2055 plotted sequentially and/or as postage stamp plots.
2056 filename: str, optional
2057 Name of the plot to write, used as a prefix for plot sequences. Defaults
2058 to the recipe name.
2059 sequence_coordinate: str, optional
2060 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2061 This coordinate must exist in the cube.
2062 stamp_coordinate: str, optional
2063 Coordinate about which to plot postage stamp plots. Defaults to
2064 ``"realization"``.
2066 Returns
2067 -------
2068 Cube
2069 The original cube (so further operations can be applied).
2071 Raises
2072 ------
2073 ValueError
2074 If the cube doesn't have the right dimensions.
2075 TypeError
2076 If the cube isn't a single cube.
2077 """
2078 _spatial_plot(
2079 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2080 )
2081 return cube
2084def spatial_multi_pcolormesh_plot(
2085 cube: iris.cube.Cube,
2086 overlay_cube: iris.cube.Cube | None = None,
2087 contour_cube: iris.cube.Cube | None = None,
2088 point_cube: iris.cube.Cube | None = None,
2089 filename: str | None = None,
2090 sequence_coordinate: str = "time",
2091 stamp_coordinate: str = "realization",
2092 **kwargs,
2093) -> iris.cube.Cube:
2094 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2096 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2097 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2098 is present then postage stamp plots will be produced.
2100 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2102 If specified, contours of a contour_cube can be overplotted on top of those.
2104 If specified, a spatial scatter map of point_cube can be overplotted.
2106 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2108 This function is significantly faster than ``spatial_contour_plot``,
2109 especially at high resolutions, and should be preferred unless contiguous
2110 contour areas are important.
2112 Parameters
2113 ----------
2114 cube: Cube
2115 Iris cube of the data to plot. It should have two spatial dimensions,
2116 such as lat and lon, and may also have two additional dimensions to be
2117 plotted sequentially and/or as postage stamp plots.
2118 overlay_cube: Cube, optional
2119 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2120 such as lat and lon, and may also have two additional dimensions to be
2121 plotted sequentially and/or as postage stamp plots. This is likely to be a masked cube in order not to hide the underlying basis cube.
2122 If not provided, output plot generated without overlay cube.
2123 contour_cube: Cube, optional
2124 Iris cube of the data to plot as a contour overlay on top of basis cube (and overlay_cube). It should have two spatial dimensions,
2125 such as lat and lon, and may also have two additional dimensions to be
2126 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2127 point_cube: Cube, optional
2128 Iris cube of the data to plot as a scatter map overlay on top of basis cube (overlay_cube and/or contour_cube). It should have two
2129 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2130 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2131 postage stamp plots. If not provided, output plot generated without point-based layer.
2132 filename: str, optional
2133 Name of the plot to write, used as a prefix for plot sequences. Defaults
2134 to the recipe name.
2135 sequence_coordinate: str, optional
2136 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2137 This coordinate must exist in the cube.
2138 stamp_coordinate: str, optional
2139 Coordinate about which to plot postage stamp plots. Defaults to
2140 ``"realization"``.
2142 Returns
2143 -------
2144 Cube
2145 The original cube (so further operations can be applied).
2147 Raises
2148 ------
2149 ValueError
2150 If the cube doesn't have the right dimensions.
2151 TypeError
2152 If the cube isn't a single cube.
2153 """
2154 _spatial_plot(
2155 "pcolormesh",
2156 cube,
2157 filename,
2158 sequence_coordinate,
2159 stamp_coordinate,
2160 overlay_cube=overlay_cube,
2161 contour_cube=contour_cube,
2162 point_cube=point_cube,
2163 )
2164 return cube, overlay_cube, contour_cube, point_cube
2167# TODO: Expand function to handle ensemble data.
2168# line_coordinate: str, optional
2169# Coordinate about which to plot multiple lines. Defaults to
2170# ``"realization"``.
2171def plot_line_series(
2172 cube: iris.cube.Cube | iris.cube.CubeList,
2173 filename: str | None = None,
2174 series_coordinate: str = "time",
2175 sequence_coordinate: str = "time",
2176 # add the following for ensembles
2177 stamp_coordinate: str = "realization",
2178 single_plot: bool = False,
2179 **kwargs,
2180) -> iris.cube.Cube | iris.cube.CubeList:
2181 """Plot a line plot for the specified coordinate.
2183 The Cube or CubeList must be 1D.
2185 Parameters
2186 ----------
2187 iris.cube | iris.cube.CubeList
2188 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2189 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2190 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2191 filename: str, optional
2192 Name of the plot to write, used as a prefix for plot sequences. Defaults
2193 to the recipe name.
2194 series_coordinate: str, optional
2195 Coordinate about which to make a series. Defaults to ``"time"``. This
2196 coordinate must exist in the cube.
2198 Returns
2199 -------
2200 iris.cube.Cube | iris.cube.CubeList
2201 The original Cube or CubeList (so further operations can be applied).
2203 Raises
2204 ------
2205 ValueError
2206 If the cubes don't have the right dimensions.
2207 TypeError
2208 If the cube isn't a Cube or CubeList.
2209 """
2210 # Ensure we have a name for the plot file.
2211 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2213 num_models = get_num_models(cube)
2215 validate_cube_shape(cube, num_models)
2217 # Iterate over all cubes and extract coordinate to plot.
2218 cubes = iris.cube.CubeList(iter_maybe(cube))
2219 coords = []
2220 for model_cube in cubes:
2221 try:
2222 coords.append(model_cube.coord(series_coordinate))
2223 except iris.exceptions.CoordinateNotFoundError as err:
2224 raise ValueError(
2225 f"Cube must have a {series_coordinate} coordinate."
2226 ) from err
2227 if model_cube.coords("realization") and model_cube.ndim > 2:
2228 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2230 plot_index = []
2232 # Check if this is a spectral plot by looking for spectral coordinates
2233 is_spectral_plot = series_coordinate in [
2234 "frequency",
2235 "physical_wavenumber",
2236 "wavelength",
2237 ]
2239 if is_spectral_plot:
2240 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2241 # coordinate frequency/wavenumber.
2242 # If several power spectra are plotted with time as sequence_coordinate for the
2243 # time slider option.
2245 # Internal plotting function.
2246 plotting_func = _plot_and_save_line_power_spectrum_series
2248 for model_cube in cubes:
2249 try:
2250 model_cube.coord(sequence_coordinate)
2251 except iris.exceptions.CoordinateNotFoundError as err:
2252 raise ValueError(
2253 f"Cube must have a {sequence_coordinate} coordinate."
2254 ) from err
2256 if num_models == 1: 2256 ↛ 2271line 2256 didn't jump to line 2271 because the condition on line 2256 was always true
2257 # check for ensembles
2258 if ( 2258 ↛ 2262line 2258 didn't jump to line 2262 because the condition on line 2258 was never true
2259 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2260 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2261 ):
2262 if single_plot:
2263 # Plot spectra, mean and ensemble spread on 1 plot
2264 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2265 else:
2266 # Plot postage stamps
2267 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2268 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2269 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2270 else:
2271 all_points = sorted(
2272 set(
2273 itertools.chain.from_iterable(
2274 cb.coord(sequence_coordinate).points for cb in cubes
2275 )
2276 )
2277 )
2278 all_slices = list(
2279 itertools.chain.from_iterable(
2280 cb.slices_over(sequence_coordinate) for cb in cubes
2281 )
2282 )
2283 # Matched slices (matched by seq coord point; it may happen that
2284 # evaluated models do not cover the same seq coord range, hence matching
2285 # necessary)
2286 cube_iterables = [
2287 iris.cube.CubeList(
2288 s
2289 for s in all_slices
2290 if s.coord(sequence_coordinate).points[0] == point
2291 )
2292 for point in all_points
2293 ]
2294 nplot = len(all_points)
2296 # Create a plot for each value of the sequence coordinate. Allowing for
2297 # multiple cubes in a CubeList to be plotted in the same plot for similar
2298 # sequence values. Passing a CubeList into the internal plotting function
2299 # for similar values of the sequence coordinate. cube_slice can be an
2300 # iris.cube.Cube or an iris.cube.CubeList.
2302 for cube_slice in cube_iterables:
2303 # Normalize cube_slice to a list of cubes
2304 if isinstance(cube_slice, iris.cube.CubeList): 2304 ↛ 2305line 2304 didn't jump to line 2305 because the condition on line 2304 was never true
2305 cubes = list(cube_slice)
2306 elif isinstance(cube_slice, iris.cube.Cube): 2306 ↛ 2309line 2306 didn't jump to line 2309 because the condition on line 2306 was always true
2307 cubes = [cube_slice]
2308 else:
2309 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2311 # Use sequence value so multiple sequences can merge.
2312 seq_coord = cube_slice[0].coord(sequence_coordinate)
2313 plot_title, plot_filename = _set_title_and_filename(
2314 seq_coord, nplot, recipe_title, filename
2315 )
2317 # Format the coordinate value in a unit appropriate way.
2318 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2320 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2321 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2321 ↛ 2322line 2321 didn't jump to line 2322 because the condition on line 2321 was never true
2322 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2324 # Do the actual plotting.
2325 plotting_func(
2326 cube_slice,
2327 coords,
2328 stamp_coordinate,
2329 plot_filename,
2330 title,
2331 series_coordinate,
2332 )
2334 plot_index.append(plot_filename)
2335 else:
2336 # Format the title and filename using plotted series coordinate
2337 nplot = 1
2338 seq_coord = coords[0]
2339 plot_title, plot_filename = _set_title_and_filename(
2340 seq_coord, nplot, recipe_title, filename
2341 )
2343 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2344 if (
2345 "station" in [c.name() for c in cubes[0].coords()]
2346 and len(cubes[0].coord("station").points) > 1
2347 ):
2348 for station in cubes[0].coord("station").points:
2349 station_cubes = cubes.extract(iris.Constraint(station=station))
2350 station_name = station_cubes[0].coord("Station_Name").points[0]
2351 station_plotname = plot_filename.replace(
2352 ".png", "_" + station_name + ".png"
2353 )
2354 _plot_and_save_line_series(
2355 station_cubes,
2356 coords,
2357 "realization",
2358 station_plotname,
2359 f"{plot_title} {station_name}",
2360 )
2361 plot_index.append(station_plotname)
2363 else:
2364 # Do the actual plotting for all other series coordinate options.
2365 _plot_and_save_line_series(
2366 cubes, coords, stamp_coordinate, plot_filename, plot_title
2367 )
2369 plot_index.append(plot_filename)
2371 # append plot to list of plots
2372 complete_plot_index = _append_to_plot_index(plot_index)
2374 # Make a page to display the plots.
2375 _make_plot_html_page(complete_plot_index)
2377 return cube
2380def plot_vertical_line_series(
2381 cubes: iris.cube.Cube | iris.cube.CubeList,
2382 filename: str | None = None,
2383 series_coordinate: str = "model_level_number",
2384 sequence_coordinate: str = "time",
2385 # line_coordinate: str = "realization",
2386 **kwargs,
2387) -> iris.cube.Cube | iris.cube.CubeList:
2388 """Plot a line plot against a type of vertical coordinate.
2390 The Cube or CubeList must be 1D.
2392 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2393 then a sequence of plots will be produced.
2395 Parameters
2396 ----------
2397 iris.cube | iris.cube.CubeList
2398 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2399 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2400 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2401 filename: str, optional
2402 Name of the plot to write, used as a prefix for plot sequences. Defaults
2403 to the recipe name.
2404 series_coordinate: str, optional
2405 Coordinate to plot on the y-axis. Can be ``pressure`` or
2406 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2407 for LFRic. Defaults to ``model_level_number``.
2408 This coordinate must exist in the cube.
2409 sequence_coordinate: str, optional
2410 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2411 This coordinate must exist in the cube.
2413 Returns
2414 -------
2415 iris.cube.Cube | iris.cube.CubeList
2416 The original Cube or CubeList (so further operations can be applied).
2417 Plotted data.
2419 Raises
2420 ------
2421 ValueError
2422 If the cubes doesn't have the right dimensions.
2423 TypeError
2424 If the cube isn't a Cube or CubeList.
2425 """
2426 # Ensure we have a name for the plot file.
2427 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2429 cubes = iter_maybe(cubes)
2430 # Initialise empty list to hold all data from all cubes in a CubeList
2431 all_data = []
2433 # Store min/max ranges for x range.
2434 x_levels = []
2436 num_models = get_num_models(cubes)
2438 validate_cube_shape(cubes, num_models)
2440 # Iterate over all cubes in cube or CubeList and plot.
2441 coords = []
2442 for cube in cubes:
2443 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2444 try:
2445 coords.append(cube.coord(series_coordinate))
2446 except iris.exceptions.CoordinateNotFoundError as err:
2447 raise ValueError(
2448 f"Cube must have a {series_coordinate} coordinate."
2449 ) from err
2451 try:
2452 if cube.ndim > 1 or not cube.coords("realization"): 2452 ↛ 2460line 2452 didn't jump to line 2460 because the condition on line 2452 was always true
2453 cube.coord(sequence_coordinate)
2454 except iris.exceptions.CoordinateNotFoundError as err:
2455 raise ValueError(
2456 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2457 ) from err
2459 # Get minimum and maximum from levels information.
2460 _, levels, _ = colorbar_map_levels(cube, axis="x")
2461 if levels is not None: 2461 ↛ 2465line 2461 didn't jump to line 2465 because the condition on line 2461 was always true
2462 x_levels.append(min(levels))
2463 x_levels.append(max(levels))
2464 else:
2465 all_data.append(cube.data)
2467 if len(x_levels) == 0: 2467 ↛ 2469line 2467 didn't jump to line 2469 because the condition on line 2467 was never true
2468 # Combine all data into a single NumPy array
2469 combined_data = np.concatenate(all_data)
2471 # Set the lower and upper limit for the x-axis to ensure all plots have
2472 # same range. This needs to read the whole cube over the range of the
2473 # sequence and if applicable postage stamp coordinate.
2474 vmin = np.floor(combined_data.min())
2475 vmax = np.ceil(combined_data.max())
2476 else:
2477 vmin = min(x_levels)
2478 vmax = max(x_levels)
2480 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2481 # a single profile directly without iterating over a sequence.
2482 sequence_coords = [
2483 cube.coord(sequence_coordinate)
2484 for cube in cubes
2485 if cube.coords(sequence_coordinate)
2486 ]
2487 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2488 np.size(coord.points) > 1 for coord in sequence_coords
2489 )
2490 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2491 np.size(coord.points) == 1 for coord in sequence_coords
2492 )
2494 plot_index = []
2495 if has_sequence_coord: 2495 ↛ 2520line 2495 didn't jump to line 2520 because the condition on line 2495 was always true
2496 # Matching the slices (matching by seq coord point; it may happen that
2497 # evaluated models do not cover the same seq coord range, hence matching
2498 # necessary)
2499 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2500 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2501 for cubes_slice in cube_iterables:
2502 # Format the coordinate value in a unit appropriate way.
2503 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2504 plot_title, plot_filename = _set_title_and_filename(
2505 seq_coord, nplot, recipe_title, filename
2506 )
2508 # Do the actual plotting.
2509 _plot_and_save_vertical_line_series(
2510 cubes_slice,
2511 coords,
2512 "realization",
2513 plot_filename,
2514 series_coordinate,
2515 title=plot_title,
2516 vmin=vmin,
2517 vmax=vmax,
2518 )
2519 plot_index.append(plot_filename)
2520 elif has_scalar_sequence_coord:
2521 # Scalar sequence coordinate (typically aggregated time bounds):
2522 # make one plot and include sequence period in title/filename.
2523 plot_title, plot_filename = _set_title_and_filename(
2524 sequence_coords[0], 1, recipe_title, filename
2525 )
2527 _plot_and_save_vertical_line_series(
2528 cubes,
2529 coords,
2530 "realization",
2531 plot_filename,
2532 series_coordinate,
2533 title=plot_title,
2534 vmin=vmin,
2535 vmax=vmax,
2536 )
2537 plot_index.append(plot_filename)
2538 else:
2539 # 1D case: no sequence coordinate, plot a single profile.
2540 plot_title = recipe_title
2541 if filename:
2542 plot_filename = filename
2543 else:
2544 plot_filename = f"{slugify(plot_title)}.png"
2546 _plot_and_save_vertical_line_series(
2547 cubes,
2548 coords,
2549 "realization",
2550 plot_filename,
2551 series_coordinate,
2552 title=plot_title,
2553 vmin=vmin,
2554 vmax=vmax,
2555 )
2556 plot_index.append(plot_filename)
2558 # Add list of plots to plot metadata.
2559 complete_plot_index = _append_to_plot_index(plot_index)
2561 # Make a page to display the plots.
2562 _make_plot_html_page(complete_plot_index)
2564 return cubes
2567def qq_plot(
2568 cubes: iris.cube.CubeList,
2569 coordinates: list[str],
2570 percentiles: list[float],
2571 model_names: list[str],
2572 filename: str | None = None,
2573 one_to_one: bool = True,
2574 **kwargs,
2575) -> iris.cube.CubeList:
2576 """Plot a Quantile-Quantile plot between two models for common time points.
2578 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2579 collapsed within the operator over all specified coordinates such as
2580 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2581 ensemble members to ensure a 1D cube (array).
2583 Parameters
2584 ----------
2585 cubes: iris.cube.CubeList
2586 Two cubes of the same variable with different models.
2587 coordinate: list[str]
2588 The list of coordinates to collapse over. This list should be
2589 every coordinate within the cube to result in a 1D cube around
2590 the percentile coordinate.
2591 percent: list[float]
2592 A list of percentiles to appear in the plot.
2593 model_names: list[str]
2594 A list of model names to appear on the axis of the plot.
2595 filename: str, optional
2596 Filename of the plot to write.
2597 one_to_one: bool, optional
2598 If True a 1:1 line is plotted; if False it is not. Default is True.
2600 Raises
2601 ------
2602 ValueError
2603 When the cubes are not compatible.
2605 Notes
2606 -----
2607 The quantile-quantile plot is a variant on the scatter plot representing
2608 two datasets by their quantiles (percentiles) for common time points.
2609 This plot does not use a theoretical distribution to compare against, but
2610 compares percentiles of two datasets. This plot does
2611 not use all raw data points, but plots the selected percentiles (quantiles) of
2612 each variable instead for the two datasets, thereby normalising the data for a
2613 direct comparison between the selected percentiles of the two dataset distributions.
2615 Quantile-quantile plots are valuable for comparing against
2616 observations and other models. Identical percentiles between the variables
2617 will lie on the one-to-one line implying the values correspond well to each
2618 other. Where there is a deviation from the one-to-one line a range of
2619 possibilities exist depending on how and where the data is shifted (e.g.,
2620 Wilks 2011 [Wilks2011]_).
2622 For distributions above the one-to-one line the distribution is left-skewed;
2623 below is right-skewed. A distinct break implies a bimodal distribution, and
2624 closer values/values further apart at the tails imply poor representation of
2625 the extremes.
2627 """
2628 # Check cubes using same functionality as the difference operator.
2629 if len(cubes) != 2:
2630 raise ValueError("cubes should contain exactly 2 cubes.")
2631 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2632 other: Cube = cubes.extract_cube(
2633 iris.Constraint(
2634 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2635 )
2636 )
2638 # Get spatial coord names.
2639 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2640 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2642 # Ensure cubes to compare are on common differencing grid.
2643 # This is triggered if either
2644 # i) latitude and longitude shapes are not the same. Note grid points
2645 # are not compared directly as these can differ through rounding
2646 # errors.
2647 # ii) or variables are known to often sit on different grid staggering
2648 # in different models (e.g. cell center vs cell edge), as is the case
2649 # for UM and LFRic comparisons.
2650 # In future greater choice of regridding method might be applied depending
2651 # on variable type. Linear regridding can in general be appropriate for smooth
2652 # variables. Care should be taken with interpretation of differences
2653 # given this dependency on regridding.
2654 if (
2655 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2656 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2657 ) or (
2658 base.long_name
2659 in [
2660 "eastward_wind_at_10m",
2661 "northward_wind_at_10m",
2662 "northward_wind_at_cell_centres",
2663 "eastward_wind_at_cell_centres",
2664 "zonal_wind_at_pressure_levels",
2665 "meridional_wind_at_pressure_levels",
2666 "potential_vorticity_at_pressure_levels",
2667 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2668 ]
2669 ):
2670 logger.debug("Linear regridding base cube to other grid to compute differences")
2671 base = regrid_onto_cube(base, other, method="Linear")
2673 # Extract just common time points.
2674 base, other = _extract_common_time_points(base, other)
2676 # Equalise attributes so we can merge.
2677 fully_equalise_attributes([base, other])
2678 logger.debug("Base: %s\nOther: %s", base, other)
2680 # Collapse cubes.
2681 base = collapse(
2682 base,
2683 coordinate=coordinates,
2684 method="PERCENTILE",
2685 additional_percent=percentiles,
2686 )
2687 other = collapse(
2688 other,
2689 coordinate=coordinates,
2690 method="PERCENTILE",
2691 additional_percent=percentiles,
2692 )
2694 # Ensure we have a name for the plot file.
2695 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2696 title = f"{recipe_title}"
2698 if filename is None:
2699 filename = slugify(recipe_title)
2701 # Add file extension.
2702 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2704 # Do the actual plotting on a scatter plot
2705 _plot_and_save_scatter_plot(
2706 base, other, plot_filename, title, one_to_one, model_names
2707 )
2709 # Add list of plots to plot metadata.
2710 plot_index = _append_to_plot_index([plot_filename])
2712 # Make a page to display the plots.
2713 _make_plot_html_page(plot_index)
2715 return iris.cube.CubeList([base, other])
2718def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2719 """
2720 Plot a Hinton style triangle/scorecard plot.
2722 This plot type can be useful for summarising high level information, such as comparing
2723 how 'skillful' two models are when verified against observations for a variety of metrics,
2724 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2725 than customisable by the user as input arguments; many have been designed to automatically
2726 scale the plot depending on the number of x and y components.
2728 Parameters
2729 ----------
2730 change: np.ndarray
2731 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2732 size/direction.
2733 signif: np.ndarray
2734 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2735 xaxis_labels: list
2736 List of labels for the xaxis (must match the second dimension length of signif and change,
2737 along with magnitude if not None).
2738 yaxis_labels: list
2739 List of labels for the yaxis (must match the first dimension length of signif and change,
2740 along with magnitude if not None).
2741 magnitude: np.ndarray | None
2742 Optional 2D array, matching the shape of change, signif, which contains numerical values
2743 the user wishes to display under each respective triangle.
2745 Returns
2746 -------
2747 matplotlib axes object to either display or do further modifications to.
2748 """
2749 # Setup colors of triangles
2750 color_pos = "#7CAE00"
2751 color_neg = "#7B68EE"
2753 # Setup cell/text size ratios
2754 figsize = None
2755 cell_size_in = 0.35
2756 text_row_ratio = 0.25
2758 # Ensure arrays, and change to bool for sig.
2759 change = np.asarray(change)
2760 signif = np.asarray(signif).astype(bool)
2761 if magnitude is not None: 2761 ↛ 2762line 2761 didn't jump to line 2762 because the condition on line 2761 was never true
2762 magnitude = np.asarray(magnitude)
2764 # Get the number of x and y elements
2765 ny, nx = change.shape
2767 # Build non-uniform y coordinates
2768 tri_height = 1.0
2769 txt_height = text_row_ratio
2771 tri_y = []
2772 txt_y = []
2773 y_edges = [0.0]
2775 y = 0.0
2776 for _j in range(ny):
2777 tri_y.append(y + tri_height / 2)
2778 y += tri_height
2779 y_edges.append(y)
2781 if magnitude is not None: 2781 ↛ 2782line 2781 didn't jump to line 2782 because the condition on line 2781 was never true
2782 txt_y.append(y + txt_height / 2)
2783 y += txt_height
2784 y_edges.append(y)
2786 total_height = y
2788 # Dynamic figure size
2789 if figsize is None: 2789 ↛ 2794line 2789 didn't jump to line 2794 because the condition on line 2789 was always true
2790 width = nx * cell_size_in
2791 height = total_height * cell_size_in + 2
2792 figsize = (width, height)
2794 fig, ax = plt.subplots(figsize=figsize)
2796 # Setup axes and grid.
2797 ax.set_aspect("equal", adjustable="box")
2798 ax.set_xlim(-0.5, nx - 0.5)
2799 ax.set_ylim(0, total_height)
2801 ax.set_xticks(np.arange(nx))
2802 ax.set_xticklabels(xaxis_labels, rotation=90)
2804 ax.set_yticks(tri_y)
2805 ax.set_yticklabels(yaxis_labels)
2807 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2808 ax.set_yticks(y_edges, minor=True)
2810 ax.set_axisbelow(True)
2811 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2812 ax.grid(False, which="major")
2813 ax.tick_params(which="minor", length=0)
2815 ax.invert_yaxis()
2817 # Compute marker scaling (fixed overlap)
2818 fig.canvas.draw()
2820 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2821 width_in, height_in = bbox.width, bbox.height
2823 cell_w = (width_in * fig.dpi) / nx
2824 cell_h = (height_in * fig.dpi) / total_height
2825 cell_pixels = min(cell_w, cell_h)
2827 max_marker_size = (0.6 * cell_pixels) ** 2
2829 text_fontsize = cell_pixels * 0.15
2831 # Plot triangles + text
2832 for j in range(ny):
2833 for i in range(nx):
2834 val = change[j, i]
2835 if np.isnan(val): 2835 ↛ 2836line 2835 didn't jump to line 2836 because the condition on line 2835 was never true
2836 continue
2838 if abs(val) < 0.01: 2838 ↛ 2839line 2838 didn't jump to line 2839 because the condition on line 2838 was never true
2839 continue
2841 sig = signif[j, i]
2842 size = max_marker_size * abs(val)
2844 # Triangle style
2845 if val >= 0:
2846 marker = "^"
2847 color = color_pos
2848 else:
2849 marker = "v"
2850 color = color_neg
2852 if sig:
2853 edgecolor = "black"
2854 linewidth = 0.6
2855 else:
2856 edgecolor = "none"
2857 linewidth = 0.0
2859 # Triangle
2860 ax.scatter(
2861 i,
2862 tri_y[j],
2863 s=size,
2864 marker=marker,
2865 c=color,
2866 edgecolors=edgecolor,
2867 linewidths=linewidth,
2868 zorder=3,
2869 clip_on=True, # ensures no rendering bleed
2870 )
2872 # Text row
2873 if magnitude is not None: 2873 ↛ 2874line 2873 didn't jump to line 2874 because the condition on line 2873 was never true
2874 mag_val = magnitude[j, i]
2876 if not np.isnan(mag_val):
2877 ax.text(
2878 i,
2879 txt_y[j],
2880 f"{mag_val:.1f}",
2881 ha="center",
2882 va="center",
2883 fontsize=text_fontsize,
2884 color="black",
2885 zorder=4,
2886 )
2888 plt.tight_layout()
2889 return fig, ax
2892def scatter_plot(
2893 cube_x: iris.cube.Cube | iris.cube.CubeList,
2894 cube_y: iris.cube.Cube | iris.cube.CubeList,
2895 filename: str | None = None,
2896 one_to_one: bool = True,
2897 **kwargs,
2898) -> iris.cube.CubeList:
2899 """Plot a scatter plot between two variables.
2901 Both cubes must be 1D.
2903 Parameters
2904 ----------
2905 cube_x: Cube | CubeList
2906 1 dimensional Cube of the data to plot on y-axis.
2907 cube_y: Cube | CubeList
2908 1 dimensional Cube of the data to plot on x-axis.
2909 filename: str, optional
2910 Filename of the plot to write.
2911 one_to_one: bool, optional
2912 If True a 1:1 line is plotted; if False it is not. Default is True.
2914 Returns
2915 -------
2916 cubes: CubeList
2917 CubeList of the original x and y cubes for further processing.
2919 Raises
2920 ------
2921 ValueError
2922 If the cube doesn't have the right dimensions and cubes not the same
2923 size.
2924 TypeError
2925 If the cube isn't a single cube.
2927 Notes
2928 -----
2929 Scatter plots are used for determining if there is a relationship between
2930 two variables. Positive relations have a slope going from bottom left to top
2931 right; Negative relations have a slope going from top left to bottom right.
2932 """
2933 # Iterate over all cubes in cube or CubeList and plot.
2934 for cube_iter in iter_maybe(cube_x):
2935 # Check cubes are correct shape.
2936 cube_iter = check_single_cube(cube_iter)
2937 if cube_iter.ndim > 1:
2938 raise ValueError("cube_x must be 1D.")
2940 # Iterate over all cubes in cube or CubeList and plot.
2941 for cube_iter in iter_maybe(cube_y):
2942 # Check cubes are correct shape.
2943 cube_iter = check_single_cube(cube_iter)
2944 if cube_iter.ndim > 1:
2945 raise ValueError("cube_y must be 1D.")
2947 # Ensure we have a name for the plot file.
2948 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2949 title = f"{recipe_title}"
2951 if filename is None:
2952 filename = slugify(recipe_title)
2954 # Add file extension.
2955 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2957 # Do the actual plotting.
2958 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2960 # Add list of plots to plot metadata.
2961 plot_index = _append_to_plot_index([plot_filename])
2963 # Make a page to display the plots.
2964 _make_plot_html_page(plot_index)
2966 return iris.cube.CubeList([cube_x, cube_y])
2969def vector_plot(
2970 cube_u: iris.cube.Cube,
2971 cube_v: iris.cube.Cube,
2972 filename: str | None = None,
2973 sequence_coordinate: str = "time",
2974 **kwargs,
2975) -> iris.cube.CubeList:
2976 """Plot a vector plot based on the input u and v components."""
2977 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2979 # Cubes must have a matching sequence coordinate.
2980 try:
2981 # Check that the u and v cubes have the same sequence coordinate.
2982 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2982 ↛ anywhereline 2982 didn't jump anywhere: it always raised an exception.
2983 raise ValueError("Coordinates do not match.")
2984 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2985 raise ValueError(
2986 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2987 ) from err
2989 # Create a plot for each value of the sequence coordinate.
2990 plot_index = []
2991 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2992 for cube_u_slice, cube_v_slice in zip(
2993 cube_u.slices_over(sequence_coordinate),
2994 cube_v.slices_over(sequence_coordinate),
2995 strict=True,
2996 ):
2997 # Format the coordinate value in a unit appropriate way.
2998 seq_coord = cube_u_slice.coord(sequence_coordinate)
2999 plot_title, plot_filename = _set_title_and_filename(
3000 seq_coord, nplot, recipe_title, filename
3001 )
3003 # Do the actual plotting.
3004 _plot_and_save_vector_plot(
3005 cube_u_slice,
3006 cube_v_slice,
3007 filename=plot_filename,
3008 title=plot_title,
3009 method="pcolormesh",
3010 )
3011 plot_index.append(plot_filename)
3013 # Add list of plots to plot metadata.
3014 complete_plot_index = _append_to_plot_index(plot_index)
3016 # Make a page to display the plots.
3017 _make_plot_html_page(complete_plot_index)
3019 return iris.cube.CubeList([cube_u, cube_v])
3022def plot_histogram_series(
3023 cubes: iris.cube.Cube | iris.cube.CubeList,
3024 filename: str | None = None,
3025 sequence_coordinate: str = "time",
3026 stamp_coordinate: str = "realization",
3027 single_plot: bool = False,
3028 **kwargs,
3029) -> iris.cube.Cube | iris.cube.CubeList:
3030 """Plot a histogram plot for each vertical level provided.
3032 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3033 is present then a sequence of plots will be produced using the time slider
3034 functionality to scroll through histograms against time. If a
3035 stamp_coordinate is present then postage stamp plots will be produced. If
3036 stamp_coordinate and single_plot is True, all postage stamp plots will be
3037 plotted in a single plot instead of separate postage stamp plots.
3039 Parameters
3040 ----------
3041 cubes: Cube | iris.cube.CubeList
3042 Iris cube or CubeList of the data to plot. It should have a single dimension other
3043 than the stamp coordinate.
3044 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3045 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3046 filename: str, optional
3047 Name of the plot to write, used as a prefix for plot sequences. Defaults
3048 to the recipe name.
3049 sequence_coordinate: str, optional
3050 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3051 This coordinate must exist in the cube and will be used for the time
3052 slider.
3053 stamp_coordinate: str, optional
3054 Coordinate about which to plot postage stamp plots. Defaults to
3055 ``"realization"``.
3056 single_plot: bool, optional
3057 If True, all postage stamp plots will be plotted in a single plot. If
3058 False, each postage stamp plot will be plotted separately. Is only valid
3059 if stamp_coordinate exists and has more than a single point.
3061 Returns
3062 -------
3063 iris.cube.Cube | iris.cube.CubeList
3064 The original Cube or CubeList (so further operations can be applied).
3065 Plotted data.
3067 Raises
3068 ------
3069 ValueError
3070 If the cube doesn't have the right dimensions.
3071 TypeError
3072 If the cube isn't a Cube or CubeList.
3073 """
3074 recipe_title = get_recipe_metadata().get("title", "Histogram")
3076 cubes = iter_maybe(cubes)
3078 # Internal plotting function.
3079 plotting_func = _plot_and_save_histogram_series
3081 num_models = get_num_models(cubes)
3083 validate_cube_shape(cubes, num_models)
3085 # If several histograms are plotted, check sequence_coordinate
3086 check_sequence_coordinate(cubes, sequence_coordinate)
3088 # Get axis minimum and maximum from levels information.
3089 # If no levels set, derive minima and maxima from data in CubeList.
3090 vmin, vmax = _set_axis_range(cubes)
3092 # Make postage stamp plots if stamp_coordinate exists and has more than a
3093 # single point. If single_plot is True:
3094 # -- all postage stamp plots will be plotted in a single plot instead of
3095 # separate postage stamp plots.
3096 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3097 # produced per single model only
3098 if num_models == 1:
3099 if ( 3099 ↛ 3103line 3099 didn't jump to line 3103 because the condition on line 3099 was never true
3100 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3101 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3102 ):
3103 if single_plot:
3104 plotting_func = (
3105 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3106 )
3107 else:
3108 plotting_func = _plot_and_save_postage_stamp_histogram_series
3109 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3110 else:
3111 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3113 plot_index = []
3114 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3115 # Create a plot for each value of the sequence coordinate. Allowing for
3116 # multiple cubes in a CubeList to be plotted in the same plot for similar
3117 # sequence values. Passing a CubeList into the internal plotting function
3118 # for similar values of the sequence coordinate. cube_slice can be an
3119 # iris.cube.Cube or an iris.cube.CubeList.
3120 for cube_slice in cube_iterables:
3121 single_cube = cube_slice
3122 if isinstance(cube_slice, iris.cube.CubeList):
3123 single_cube = cube_slice[0]
3125 # Ensure valid stamp coordinate in cube dimensions
3126 if stamp_coordinate == "realization": 3126 ↛ 3129line 3126 didn't jump to line 3129 because the condition on line 3126 was always true
3127 stamp_coordinate = check_stamp_coordinate(single_cube)
3128 # Set plot titles and filename, based on sequence coordinate
3129 seq_coord = single_cube.coord(sequence_coordinate)
3130 # Use time coordinate in title and filename if single histogram output.
3131 if sequence_coordinate == "realization" and nplot == 1: 3131 ↛ 3132line 3131 didn't jump to line 3132 because the condition on line 3131 was never true
3132 seq_coord = single_cube.coord("time")
3133 # Use station name in title and filename if model vs obs comparison
3134 if sequence_coordinate == "station": 3134 ↛ 3135line 3134 didn't jump to line 3135 because the condition on line 3134 was never true
3135 seq_coord = single_cube.coord("Station_Name")
3137 plot_title, plot_filename = _set_title_and_filename(
3138 seq_coord, nplot, recipe_title, filename
3139 )
3141 # Do the actual plotting.
3142 plotting_func(
3143 cube_slice,
3144 filename=plot_filename,
3145 stamp_coordinate=stamp_coordinate,
3146 title=plot_title,
3147 vmin=vmin,
3148 vmax=vmax,
3149 )
3150 plot_index.append(plot_filename)
3152 # Add list of plots to plot metadata.
3153 complete_plot_index = _append_to_plot_index(plot_index)
3155 # Make a page to display the plots.
3156 _make_plot_html_page(complete_plot_index)
3158 return cubes
3161def plot_scatter_series(
3162 cubes: iris.cube.Cube | iris.cube.CubeList,
3163 filename: str | None = None,
3164 sequence_coordinate: str = "time",
3165 stamp_coordinate: str = "realization",
3166 hexbin: bool = False,
3167 **kwargs,
3168) -> iris.cube.Cube | iris.cube.CubeList:
3169 """Plot a scatter plot for each sequence coordinate provided.
3171 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3172 is present then a sequence of plots will be produced using the time slider
3173 functionality to scroll through scatter against time. If a
3174 stamp_coordinate is present then postage stamp plots will be produced. If
3175 stamp_coordinate and single_plot is True, all postage stamp plots will be
3176 plotted in a single plot instead of separate postage stamp plots.
3178 Parameters
3179 ----------
3180 cubes: Cube | iris.cube.CubeList
3181 Iris cube or CubeList of the data to plot. It should have a single dimension other
3182 than the stamp coordinate.
3183 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3184 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3185 filename: str, optional
3186 Name of the plot to write, used as a prefix for plot sequences. Defaults
3187 to the recipe name.
3188 sequence_coordinate: str, optional
3189 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3190 This coordinate must exist in the cube and will be used for the time
3191 slider.
3192 stamp_coordinate: str, optional
3193 Coordinate about which to plot postage stamp plots. Defaults to
3194 ``"realization"``.
3195 hexbin: bool, optional
3196 If True, generate hexbin comparison plot.
3197 If False, generate point-by-point scatter plot.
3199 Returns
3200 -------
3201 iris.cube.Cube | iris.cube.CubeList
3202 The original Cube or CubeList (so further operations can be applied).
3203 Plotted data.
3205 Raises
3206 ------
3207 ValueError
3208 If the cube doesn't have the right dimensions.
3209 TypeError
3210 If the cube isn't a Cube or CubeList.
3211 """
3212 recipe_title = get_recipe_metadata().get("title", "Scatter")
3214 cubes = iter_maybe(cubes)
3216 # Internal plotting function.
3217 plotting_func = _plot_and_save_scatter_series
3219 num_models = get_num_models(cubes)
3221 validate_cube_shape(cubes, num_models)
3223 check_sequence_coordinate(cubes, sequence_coordinate)
3225 vmin, vmax = _set_axis_range(cubes)
3227 # Require >1 models to compare on scatter plot
3228 if num_models > 1:
3229 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3230 else:
3231 raise ValueError(
3232 "Scatter plot series requires multiple number of models in input data."
3233 )
3235 plot_index = []
3236 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3237 # Create a plot for each value of the sequence coordinate. Allowing for
3238 # multiple cubes in a CubeList to be plotted in the same plot for similar
3239 # sequence values. Passing a CubeList into the internal plotting function
3240 # for similar values of the sequence coordinate. cube_slice can be an
3241 # iris.cube.Cube or an iris.cube.CubeList.
3242 for cube_slice in cube_iterables:
3243 single_cube = cube_slice
3244 if isinstance(cube_slice, iris.cube.CubeList): 3244 ↛ 3248line 3244 didn't jump to line 3248 because the condition on line 3244 was always true
3245 single_cube = cube_slice[0]
3247 # Ensure valid stamp coordinate in cube dimensions
3248 if stamp_coordinate == "realization": 3248 ↛ 3251line 3248 didn't jump to line 3251 because the condition on line 3248 was always true
3249 stamp_coordinate = check_stamp_coordinate(single_cube)
3250 # Set plot titles and filename, based on sequence coordinate
3251 seq_coord = single_cube.coord(sequence_coordinate)
3252 # Use time coordinate in title and filename if single histogram output.
3253 if sequence_coordinate == "realization" and nplot == 1:
3254 seq_coord = single_cube.coord("time")
3255 # Use station name in title and filename if model vs obs comparison
3256 if sequence_coordinate == "station":
3257 seq_coord = single_cube.coord("Station_Name")
3259 plot_title, plot_filename = _set_title_and_filename(
3260 seq_coord, nplot, recipe_title, filename
3261 )
3263 # Do the actual plotting.
3264 plotting_func(
3265 cube_slice,
3266 filename=plot_filename,
3267 stamp_coordinate=stamp_coordinate,
3268 title=plot_title,
3269 vmin=vmin,
3270 vmax=vmax,
3271 hexbin=hexbin,
3272 )
3273 plot_index.append(plot_filename)
3275 # Add list of plots to plot metadata.
3276 complete_plot_index = _append_to_plot_index(plot_index)
3278 # Make a page to display the plots.
3279 _make_plot_html_page(complete_plot_index)
3281 return cubes
3284def _plot_and_save_postage_stamp_power_spectrum_series(
3285 cubes: iris.cube.Cube,
3286 coords: list[iris.coords.Coord],
3287 stamp_coordinate: str,
3288 filename: str,
3289 title: str,
3290 series_coordinate: str | None = None,
3291 **kwargs,
3292):
3293 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3295 Parameters
3296 ----------
3297 cubes: Cube or CubeList
3298 Cube or Cubelist of the power spectrum data.
3299 coords: list[Coord]
3300 Coordinates to plot on the x-axis, one per cube.
3301 stamp_coordinate: str
3302 Coordinate that becomes different plots.
3303 filename: str
3304 Filename of the plot to write.
3305 title: str
3306 Plot title.
3307 series_coordinate: str, optional
3308 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3310 """
3311 # Use the smallest square grid that will fit the members.
3312 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3314 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3315 model_colors_map = get_model_colors_map(cubes)
3316 # ax = plt.gca()
3317 # Make a subplot for each member.
3318 for member, subplot in zip(
3319 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3320 ):
3321 ax = plt.subplot(grid_size, grid_size, subplot)
3323 # Store min/max ranges.
3324 y_levels = []
3326 line_marker = None
3327 line_width = 1
3329 for cube in iter_maybe(member):
3330 xcoord = _select_series_coord(cube, series_coordinate)
3331 xname = xcoord.points
3333 yfield = cube.data # power spectrum
3334 label = None
3335 color = "black"
3336 if model_colors_map: 3336 ↛ 3337line 3336 didn't jump to line 3337 because the condition on line 3336 was never true
3337 label = cube.attributes.get("model_name")
3338 color = model_colors_map.get(label)
3340 if member.coord(stamp_coordinate).points == [0]:
3341 ax.plot(
3342 xname,
3343 yfield,
3344 color=color,
3345 marker=line_marker,
3346 ls="-",
3347 lw=line_width,
3348 label=f"{label} (control)"
3349 if len(cube.coord(stamp_coordinate).points) > 1
3350 else label,
3351 )
3352 # Label with member if part of an ensemble and not the control.
3353 else:
3354 ax.plot(
3355 xname,
3356 yfield,
3357 color=color,
3358 ls="-",
3359 lw=1.5,
3360 alpha=0.75,
3361 label=f"{label} (member)",
3362 )
3364 # Calculate the global min/max if multiple cubes are given.
3365 _, levels, _ = colorbar_map_levels(cube, axis="y")
3366 if levels is not None: 3366 ↛ 3367line 3366 didn't jump to line 3367 because the condition on line 3366 was never true
3367 y_levels.append(min(levels))
3368 y_levels.append(max(levels))
3370 # Add some labels and tweak the style.
3371 title = f"{title}"
3372 ax.set_title(title, fontsize=16)
3374 # Set appropriate x-axis label based on coordinate
3375 if series_coordinate == "wavelength" or ( 3375 ↛ 3378line 3375 didn't jump to line 3378 because the condition on line 3375 was never true
3376 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3377 ):
3378 ax.set_xlabel("Wavelength (km)", fontsize=14)
3379 elif series_coordinate == "physical_wavenumber" or ( 3379 ↛ 3384line 3379 didn't jump to line 3384 because the condition on line 3379 was always true
3380 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3381 ):
3382 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3383 else: # frequency or check units
3384 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3385 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3386 else:
3387 ax.set_xlabel("Wavenumber", fontsize=14)
3389 ax.set_ylabel("Power Spectral Density", fontsize=14)
3390 ax.tick_params(axis="both", labelsize=12)
3392 # Set log-log scale
3393 ax.set_xscale("log")
3394 ax.set_yscale("log")
3396 # Add gridlines
3397 ax.grid(linestyle="--", color="grey", linewidth=1)
3398 # Ientify unique labels for legend
3399 handles = list(
3400 {
3401 label: handle
3402 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3403 }.values()
3404 )
3405 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3407 ax = plt.gca()
3408 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3410 # Save plot.
3411 _save_close_figure(fig, "histogram postage stamp", filename)
3414def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3415 cubes: iris.cube.Cube,
3416 coords: list[iris.coords.Coord],
3417 stamp_coordinate: str,
3418 filename: str,
3419 title: str,
3420 series_coordinate: str | None = None,
3421 **kwargs,
3422):
3423 """Plot and save power spectra for ensemble members in single plot.
3425 Parameters
3426 ----------
3427 cubes: Cube or CubeList
3428 Cube or Cubelist of the power spectrum data.
3429 coords: list[Coord]
3430 Coordinates to plot on the x-axis, one per cube.
3431 stamp_coordinate: str
3432 Coordinate that becomes different plots.
3433 filename: str
3434 Filename of the plot to write.
3435 title: str
3436 Plot title.
3437 series_coordinate: str, optional
3438 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3440 """
3441 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3442 model_colors_map = get_model_colors_map(cubes)
3444 line_marker = None
3445 line_width = 1
3447 # Compute ensemble statistics to show spread
3448 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3449 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3450 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3452 xcoord_global = mean_cube.coord(series_coordinate)
3453 x_global = xcoord_global.points
3455 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3456 xcoord = _select_series_coord(member, series_coordinate)
3457 xname = xcoord.points
3459 yfield = member.data # power spectrum
3460 color = "black"
3461 if model_colors_map: 3461 ↛ 3465line 3461 didn't jump to line 3465 because the condition on line 3461 was always true
3462 label = member.attributes.get("model_name") if i == 0 else None
3463 color = model_colors_map.get(label)
3465 if member.coord(stamp_coordinate).points == [0]:
3466 ax.plot(
3467 xname,
3468 yfield,
3469 color=color,
3470 marker=line_marker,
3471 ls="-",
3472 lw=line_width,
3473 label=f"{label} (control)"
3474 if len(member.coord(stamp_coordinate).points) > 1
3475 else label,
3476 )
3477 # Label with member number if part of an ensemble and not the control.
3478 else:
3479 ax.plot(
3480 xname,
3481 yfield,
3482 color=color,
3483 ls="-",
3484 lw=1.5,
3485 alpha=0.75,
3486 label=label,
3487 )
3489 # Set appropriate x-axis label based on coordinate
3490 if series_coordinate == "wavelength" or ( 3490 ↛ 3493line 3490 didn't jump to line 3493 because the condition on line 3490 was never true
3491 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3492 ):
3493 ax.set_xlabel("Wavelength (km)", fontsize=14)
3494 elif series_coordinate == "physical_wavenumber" or ( 3494 ↛ 3499line 3494 didn't jump to line 3499 because the condition on line 3494 was always true
3495 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3496 ):
3497 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3498 else: # frequency or check units
3499 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3500 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3501 else:
3502 ax.set_xlabel("Wavenumber", fontsize=14)
3504 # Add ensemble spread shading
3505 ax.fill_between(
3506 x_global,
3507 min_cube.data,
3508 max_cube.data,
3509 color="grey",
3510 alpha=0.3,
3511 label="Ensemble spread",
3512 )
3514 # Add ensemble mean line
3515 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3517 ax.set_ylabel("Power Spectral Density", fontsize=14)
3518 ax.tick_params(axis="both", labelsize=12)
3520 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3521 # Set log-log scale
3522 ax.set_xscale("log")
3523 ax.set_yscale("log")
3525 # Add gridlines
3526 ax.grid(linestyle="--", color="grey", linewidth=1)
3527 # Identify unique labels for legend
3528 handles = list(
3529 {
3530 label: handle
3531 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3532 }.values()
3533 )
3534 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3536 # Figure title.
3537 ax.set_title(title, fontsize=16)
3539 # Save plot.
3540 _save_close_figure(fig, "power spectra postage stamp", filename)