Coverage for src/CSET/operators/plot.py: 83%
1100 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-06 10:09 +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.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(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.min(cube_vec_mag.data):.3g} Max: {np.max(cube_vec_mag.data):.3g} Mean: {np.mean(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 ax.set_xlim(vmin, vmax)
1609 ax.tick_params(axis="both", labelsize=12)
1611 # Overlay grid-lines onto histogram plot.
1612 ax.grid(linestyle="--", color="grey", linewidth=1)
1613 if model_colors_map:
1614 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1616 # Save plot.
1617 _save_close_figure(fig, "histogram", filename)
1620def _plot_and_save_postage_stamp_histogram_series(
1621 cube: iris.cube.Cube,
1622 filename: str,
1623 title: str,
1624 stamp_coordinate: str,
1625 vmin: float,
1626 vmax: float,
1627 **kwargs,
1628):
1629 """Plot and save postage (ensemble members) stamps for a histogram series.
1631 Parameters
1632 ----------
1633 cube: Cube
1634 2 dimensional Cube of the data to plot as histogram.
1635 filename: str
1636 Filename of the plot to write.
1637 title: str
1638 Plot title.
1639 stamp_coordinate: str
1640 Coordinate that becomes different plots.
1641 vmin: float
1642 minimum for pdf x-axis
1643 vmax: float
1644 maximum for pdf x-axis
1645 """
1646 # Use the smallest square grid that will fit the members.
1647 nmember = len(cube.coord(stamp_coordinate).points)
1648 grid_rows = int(math.sqrt(nmember))
1649 grid_size = math.ceil(nmember / grid_rows)
1651 fig = plt.figure(
1652 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1653 )
1654 # Make a subplot for each member.
1655 for member, subplot in zip(
1656 cube.slices_over(stamp_coordinate),
1657 range(1, grid_size * grid_rows + 1),
1658 strict=False,
1659 ):
1660 # Implicit interface is much easier here, due to needing to have the
1661 # cartopy GeoAxes generated.
1662 plt.subplot(grid_rows, grid_size, subplot)
1663 # Reshape cube data into a single array to allow for a single histogram.
1664 # Otherwise we plot xdim histograms stacked.
1665 member_data_1d = (member.data).flatten()
1666 plt.hist(member_data_1d, density=True, stacked=True)
1667 axes = plt.gca()
1668 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1669 axes.set_title(f"{mtitle}")
1670 axes.set_xlim(vmin, vmax)
1672 # Overall figure title.
1673 fig.suptitle(title, fontsize=16)
1675 # Save plot.
1676 _save_close_figure(fig, "histogram postage stamp", filename)
1679def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1680 cube: iris.cube.Cube,
1681 filename: str,
1682 title: str,
1683 stamp_coordinate: str,
1684 vmin: float,
1685 vmax: float,
1686 **kwargs,
1687):
1688 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1689 ax.set_title(title, fontsize=16)
1690 ax.set_xlim(vmin, vmax)
1691 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1692 ax.set_ylabel("normalised probability density", fontsize=14)
1693 # Loop over all slices along the stamp_coordinate
1694 for member in cube.slices_over(stamp_coordinate):
1695 # Flatten the member data to 1D
1696 member_data_1d = member.data.flatten()
1697 # Plot the histogram using plt.hist
1698 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1699 plt.hist(
1700 member_data_1d,
1701 density=True,
1702 stacked=True,
1703 label=f"{mtitle}",
1704 )
1706 # Add a legend
1707 ax.legend(fontsize=16)
1709 # Save plot.
1710 _save_close_figure(fig, "histogram postage stamp", filename)
1713def _plot_and_save_scatter_series(
1714 cubes: iris.cube.Cube | iris.cube.CubeList,
1715 filename: str,
1716 title: str,
1717 vmin: float,
1718 vmax: float,
1719 hexbin: bool,
1720 **kwargs,
1721):
1722 """Plot and save a scatter plot series.
1724 Parameters
1725 ----------
1726 cubes: Cube or CubeList
1727 2 dimensional Cube or CubeList of the data to plot as scatter.
1728 filename: str
1729 Filename of the plot to write.
1730 title: str
1731 Plot title.
1732 vmin: float
1733 minimum for colorbar
1734 vmax: float
1735 maximum for colorbar
1736 hexbin: bool
1737 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1738 Else scatter of all points, with potential to overplot many comparisons on same plot.
1739 """
1740 if hexbin:
1741 # Check cubes using same functionality as the difference operator.
1742 if len(cubes) != 2:
1743 raise ValueError(
1744 "Cubes should contain exactly 2 cubes for hexbin plotting."
1745 )
1746 title = title.replace("scatter", "hexbin")
1747 filename = filename.replace("scatter", "hexbin")
1749 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1750 ax = plt.gca()
1752 model_colors_map = get_model_colors_map(cubes)
1754 percentiles = np.arange(0, 100, 5)
1755 percentiles[0] = 1
1756 percentiles[-1] = 99
1757 quantiles = iris.cube.CubeList()
1759 # Loop through all output cubes for both data points and overplotting quantiles.
1760 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1761 for plottype in ["points", "quantiles"]:
1762 nplot = 0
1763 for cube in iter_maybe(cubes):
1764 label = None
1765 color = "black"
1766 if model_colors_map: 1766 ↛ 1771line 1766 didn't jump to line 1771 because the condition on line 1766 was always true
1767 label = cube.attributes.get("model_name")
1768 color = model_colors_map[label]
1770 # Plot all data points
1771 if plottype == "points":
1772 if nplot > 0:
1773 if hexbin:
1774 hb = plt.hexbin(
1775 cubes[0].data.flatten(),
1776 cube.data.flatten(),
1777 alpha=0.3,
1778 gridsize=100,
1779 mincnt=1,
1780 )
1781 else:
1782 plt.scatter(
1783 cubes[0].data.flatten(),
1784 cube.data.flatten(),
1785 color=color,
1786 marker="+",
1787 label=None,
1788 alpha=0.3,
1789 )
1791 elif plottype == "quantiles": 1791 ↛ 1810line 1791 didn't jump to line 1810 because the condition on line 1791 was always true
1792 # Construct Q-Q plot
1793 quantiles.append(
1794 cube.collapsed(
1795 cube.coords(dim_coords=True),
1796 iris.analysis.PERCENTILE,
1797 percent=percentiles,
1798 )
1799 )
1800 if nplot > 0:
1801 iplt.scatter(
1802 quantiles[0],
1803 quantiles[-1],
1804 color=color,
1805 marker="o",
1806 label=label,
1807 edgecolors="black",
1808 )
1810 nplot = nplot + 1
1812 # Add some labels and tweak the style.
1813 ax.set_title(title, fontsize=16)
1814 ax.set_xlabel(
1815 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1816 )
1817 ax.set_ylabel(
1818 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1819 )
1820 ax.tick_params(axis="both", labelsize=12)
1821 ax.autoscale()
1823 # Set 1:1 line and equal axes if scatter plot of common cube names
1824 nameA = iter_maybe(cubes)[0].name()
1825 nameB = iter_maybe(cubes)[1].name()
1826 if any(part in nameB.split("_") for part in nameA.split("_")): 1826 ↛ 1837line 1826 didn't jump to line 1837 because the condition on line 1826 was always true
1827 lims = [
1828 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1829 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1830 ]
1831 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1832 ax.set_aspect("equal")
1833 ax.set_xlim(lims)
1834 ax.set_ylim(lims)
1836 # Overlay grid-lines onto scatter plot.
1837 ax.grid(linestyle="--", color="grey", linewidth=1)
1838 if model_colors_map: 1838 ↛ 1842line 1838 didn't jump to line 1842 because the condition on line 1838 was always true
1839 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1841 # Add colorbar if hexbin output
1842 if hexbin:
1843 cb = plt.colorbar(
1844 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1845 )
1846 cb.set_label("Number of data points", size=12)
1848 # Save plot.
1849 _save_close_figure(fig, "scatter", filename)
1852def _spatial_plot(
1853 method: Literal["contourf", "pcolormesh", "scatter"],
1854 cube: iris.cube.Cube,
1855 filename: str | None,
1856 sequence_coordinate: str,
1857 stamp_coordinate: str,
1858 overlay_cube: iris.cube.Cube | None = None,
1859 contour_cube: iris.cube.Cube | None = None,
1860 point_cube: iris.cube.Cube | None = None,
1861 **kwargs,
1862):
1863 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1865 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1866 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1867 is present then postage stamp plots will be produced.
1869 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1870 be overplotted on the same figure.
1872 Parameters
1873 ----------
1874 method: "contourf" | "pcolormesh" | "scatter"
1875 The plotting method to use.
1876 Select choice of "contourf" or "pcolormesh" for gridded data.
1877 Use "scatter" for point-based data.
1878 cube: Cube
1879 Iris cube of the data to plot. It should have two spatial dimensions,
1880 such as lat and lon, and may also have a another two dimension to be
1881 plotted sequentially and/or as postage stamp plots.
1882 filename: str | None
1883 Name of the plot to write, used as a prefix for plot sequences. If None
1884 uses the recipe name.
1885 sequence_coordinate: str
1886 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1887 This coordinate must exist in the cube.
1888 stamp_coordinate: str
1889 Coordinate about which to plot postage stamp plots. Defaults to
1890 ``"realization"``.
1891 overlay_cube: Cube | None, optional
1892 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1893 contour_cube: Cube | None, optional
1894 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1895 point_cube: Cube | None, optional
1896 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
1898 Raises
1899 ------
1900 ValueError
1901 If the cube doesn't have the right dimensions.
1902 TypeError
1903 If the cube isn't a single cube.
1904 """
1905 # Ensure we've got a single cube.
1906 cube = check_single_cube(cube)
1908 # Set title based on recipe metadata or use cube name
1909 recipe_title = get_recipe_metadata().get("title", cube.name())
1911 # Check if there is a valid stamp coordinate in cube dimensions.
1912 if stamp_coordinate == "realization": 1912 ↛ 1917line 1912 didn't jump to line 1917 because the condition on line 1912 was always true
1913 stamp_coordinate = check_stamp_coordinate(cube)
1915 # Make postage stamp plots if stamp_coordinate exists and has more than a
1916 # single point.
1917 plotting_func = _plot_and_save_spatial_plot
1918 try:
1919 if cube.coord(stamp_coordinate).shape[0] > 1:
1920 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1921 except iris.exceptions.CoordinateNotFoundError:
1922 pass
1924 # Produce a geographical scatter plot if the data have a
1925 # dimension called observation or model_obs_error
1926 if any(
1927 crd.var_name == "station"
1928 or crd.var_name == "Station_Name"
1929 or crd.var_name == "model_obs_error"
1930 for crd in cube.coords()
1931 ):
1932 plotting_func = _plot_and_save_spatial_plot
1933 method = "scatter"
1935 # Must have a sequence coordinate.
1936 try:
1937 cube.coord(sequence_coordinate)
1938 except iris.exceptions.CoordinateNotFoundError as err:
1939 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1941 # Create a plot for each value of the sequence coordinate.
1942 plot_index = []
1943 nplot = np.size(cube.coord(sequence_coordinate).points)
1945 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1946 # Set plot titles and filename
1947 seq_coord = cube_slice.coord(sequence_coordinate)
1948 plot_title, plot_filename = _set_title_and_filename(
1949 seq_coord, nplot, recipe_title, filename
1950 )
1952 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1953 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1954 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1955 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1957 # Do the actual plotting.
1958 plotting_func(
1959 cube_slice,
1960 filename=plot_filename,
1961 stamp_coordinate=stamp_coordinate,
1962 title=plot_title,
1963 method=method,
1964 overlay_cube=overlay_slice,
1965 contour_cube=contour_slice,
1966 point_cube=point_slice,
1967 **kwargs,
1968 )
1969 plot_index.append(plot_filename)
1971 # Add list of plots to plot metadata.
1972 complete_plot_index = _append_to_plot_index(plot_index)
1974 # Make a page to display the plots.
1975 _make_plot_html_page(complete_plot_index)
1978####################
1979# Public functions #
1980####################
1983def spatial_contour_plot(
1984 cube: iris.cube.Cube,
1985 filename: str | None = None,
1986 sequence_coordinate: str = "time",
1987 stamp_coordinate: str = "realization",
1988 **kwargs,
1989) -> iris.cube.Cube:
1990 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1992 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1993 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1994 is present then postage stamp plots will be produced.
1996 Parameters
1997 ----------
1998 cube: Cube
1999 Iris cube of the data to plot. It should have two spatial dimensions,
2000 such as lat and lon, and may also have a another two dimension to be
2001 plotted sequentially and/or as postage stamp plots.
2002 filename: str, optional
2003 Name of the plot to write, used as a prefix for plot sequences. Defaults
2004 to the recipe name.
2005 sequence_coordinate: str, optional
2006 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2007 This coordinate must exist in the cube.
2008 stamp_coordinate: str, optional
2009 Coordinate about which to plot postage stamp plots. Defaults to
2010 ``"realization"``.
2012 Returns
2013 -------
2014 Cube
2015 The original cube (so further operations can be applied).
2017 Raises
2018 ------
2019 ValueError
2020 If the cube doesn't have the right dimensions.
2021 TypeError
2022 If the cube isn't a single cube.
2023 """
2024 _spatial_plot(
2025 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2026 )
2027 return cube
2030def spatial_pcolormesh_plot(
2031 cube: iris.cube.Cube,
2032 filename: str | None = None,
2033 sequence_coordinate: str = "time",
2034 stamp_coordinate: str = "realization",
2035 **kwargs,
2036) -> iris.cube.Cube:
2037 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2039 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2040 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2041 is present then postage stamp plots will be produced.
2043 This function is significantly faster than ``spatial_contour_plot``,
2044 especially at high resolutions, and should be preferred unless contiguous
2045 contour areas are important.
2047 Parameters
2048 ----------
2049 cube: Cube
2050 Iris cube of the data to plot. It should have two spatial dimensions,
2051 such as lat and lon, and may also have a another two dimension to be
2052 plotted sequentially and/or as postage stamp plots.
2053 filename: str, optional
2054 Name of the plot to write, used as a prefix for plot sequences. Defaults
2055 to the recipe name.
2056 sequence_coordinate: str, optional
2057 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2058 This coordinate must exist in the cube.
2059 stamp_coordinate: str, optional
2060 Coordinate about which to plot postage stamp plots. Defaults to
2061 ``"realization"``.
2063 Returns
2064 -------
2065 Cube
2066 The original cube (so further operations can be applied).
2068 Raises
2069 ------
2070 ValueError
2071 If the cube doesn't have the right dimensions.
2072 TypeError
2073 If the cube isn't a single cube.
2074 """
2075 _spatial_plot(
2076 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2077 )
2078 return cube
2081def spatial_multi_pcolormesh_plot(
2082 cube: iris.cube.Cube,
2083 overlay_cube: iris.cube.Cube | None = None,
2084 contour_cube: iris.cube.Cube | None = None,
2085 point_cube: iris.cube.Cube | None = None,
2086 filename: str | None = None,
2087 sequence_coordinate: str = "time",
2088 stamp_coordinate: str = "realization",
2089 **kwargs,
2090) -> iris.cube.Cube:
2091 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2093 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2094 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2095 is present then postage stamp plots will be produced.
2097 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2099 If specified, contours of a contour_cube can be overplotted on top of those.
2101 If specified, a spatial scatter map of point_cube can be overplotted.
2103 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2105 This function is significantly faster than ``spatial_contour_plot``,
2106 especially at high resolutions, and should be preferred unless contiguous
2107 contour areas are important.
2109 Parameters
2110 ----------
2111 cube: Cube
2112 Iris cube of the data to plot. It should have two spatial dimensions,
2113 such as lat and lon, and may also have two additional dimensions to be
2114 plotted sequentially and/or as postage stamp plots.
2115 overlay_cube: Cube, optional
2116 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2117 such as lat and lon, and may also have two additional dimensions to be
2118 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.
2119 If not provided, output plot generated without overlay cube.
2120 contour_cube: Cube, optional
2121 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,
2122 such as lat and lon, and may also have two additional dimensions to be
2123 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2124 point_cube: Cube, optional
2125 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
2126 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2127 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2128 postage stamp plots. If not provided, output plot generated without point-based layer.
2129 filename: str, optional
2130 Name of the plot to write, used as a prefix for plot sequences. Defaults
2131 to the recipe name.
2132 sequence_coordinate: str, optional
2133 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2134 This coordinate must exist in the cube.
2135 stamp_coordinate: str, optional
2136 Coordinate about which to plot postage stamp plots. Defaults to
2137 ``"realization"``.
2139 Returns
2140 -------
2141 Cube
2142 The original cube (so further operations can be applied).
2144 Raises
2145 ------
2146 ValueError
2147 If the cube doesn't have the right dimensions.
2148 TypeError
2149 If the cube isn't a single cube.
2150 """
2151 _spatial_plot(
2152 "pcolormesh",
2153 cube,
2154 filename,
2155 sequence_coordinate,
2156 stamp_coordinate,
2157 overlay_cube=overlay_cube,
2158 contour_cube=contour_cube,
2159 point_cube=point_cube,
2160 )
2161 return cube, overlay_cube, contour_cube, point_cube
2164# TODO: Expand function to handle ensemble data.
2165# line_coordinate: str, optional
2166# Coordinate about which to plot multiple lines. Defaults to
2167# ``"realization"``.
2168def plot_line_series(
2169 cube: iris.cube.Cube | iris.cube.CubeList,
2170 filename: str | None = None,
2171 series_coordinate: str = "time",
2172 sequence_coordinate: str = "time",
2173 # add the following for ensembles
2174 stamp_coordinate: str = "realization",
2175 single_plot: bool = False,
2176 **kwargs,
2177) -> iris.cube.Cube | iris.cube.CubeList:
2178 """Plot a line plot for the specified coordinate.
2180 The Cube or CubeList must be 1D.
2182 Parameters
2183 ----------
2184 iris.cube | iris.cube.CubeList
2185 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2186 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2187 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2188 filename: str, optional
2189 Name of the plot to write, used as a prefix for plot sequences. Defaults
2190 to the recipe name.
2191 series_coordinate: str, optional
2192 Coordinate about which to make a series. Defaults to ``"time"``. This
2193 coordinate must exist in the cube.
2195 Returns
2196 -------
2197 iris.cube.Cube | iris.cube.CubeList
2198 The original Cube or CubeList (so further operations can be applied).
2200 Raises
2201 ------
2202 ValueError
2203 If the cubes don't have the right dimensions.
2204 TypeError
2205 If the cube isn't a Cube or CubeList.
2206 """
2207 # Ensure we have a name for the plot file.
2208 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2210 num_models = get_num_models(cube)
2212 validate_cube_shape(cube, num_models)
2214 # Iterate over all cubes and extract coordinate to plot.
2215 cubes = iris.cube.CubeList(iter_maybe(cube))
2216 coords = []
2217 for model_cube in cubes:
2218 try:
2219 coords.append(model_cube.coord(series_coordinate))
2220 except iris.exceptions.CoordinateNotFoundError as err:
2221 raise ValueError(
2222 f"Cube must have a {series_coordinate} coordinate."
2223 ) from err
2224 if model_cube.coords("realization") and model_cube.ndim > 2:
2225 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2227 plot_index = []
2229 # Check if this is a spectral plot by looking for spectral coordinates
2230 is_spectral_plot = series_coordinate in [
2231 "frequency",
2232 "physical_wavenumber",
2233 "wavelength",
2234 ]
2236 if is_spectral_plot:
2237 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2238 # coordinate frequency/wavenumber.
2239 # If several power spectra are plotted with time as sequence_coordinate for the
2240 # time slider option.
2242 # Internal plotting function.
2243 plotting_func = _plot_and_save_line_power_spectrum_series
2245 for model_cube in cubes:
2246 try:
2247 model_cube.coord(sequence_coordinate)
2248 except iris.exceptions.CoordinateNotFoundError as err:
2249 raise ValueError(
2250 f"Cube must have a {sequence_coordinate} coordinate."
2251 ) from err
2253 if num_models == 1: 2253 ↛ 2268line 2253 didn't jump to line 2268 because the condition on line 2253 was always true
2254 # check for ensembles
2255 if ( 2255 ↛ 2259line 2255 didn't jump to line 2259 because the condition on line 2255 was never true
2256 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2257 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2258 ):
2259 if single_plot:
2260 # Plot spectra, mean and ensemble spread on 1 plot
2261 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2262 else:
2263 # Plot postage stamps
2264 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2265 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2266 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2267 else:
2268 all_points = sorted(
2269 set(
2270 itertools.chain.from_iterable(
2271 cb.coord(sequence_coordinate).points for cb in cubes
2272 )
2273 )
2274 )
2275 all_slices = list(
2276 itertools.chain.from_iterable(
2277 cb.slices_over(sequence_coordinate) for cb in cubes
2278 )
2279 )
2280 # Matched slices (matched by seq coord point; it may happen that
2281 # evaluated models do not cover the same seq coord range, hence matching
2282 # necessary)
2283 cube_iterables = [
2284 iris.cube.CubeList(
2285 s
2286 for s in all_slices
2287 if s.coord(sequence_coordinate).points[0] == point
2288 )
2289 for point in all_points
2290 ]
2291 nplot = len(all_points)
2293 # Create a plot for each value of the sequence coordinate. Allowing for
2294 # multiple cubes in a CubeList to be plotted in the same plot for similar
2295 # sequence values. Passing a CubeList into the internal plotting function
2296 # for similar values of the sequence coordinate. cube_slice can be an
2297 # iris.cube.Cube or an iris.cube.CubeList.
2299 for cube_slice in cube_iterables:
2300 # Normalize cube_slice to a list of cubes
2301 if isinstance(cube_slice, iris.cube.CubeList): 2301 ↛ 2302line 2301 didn't jump to line 2302 because the condition on line 2301 was never true
2302 cubes = list(cube_slice)
2303 elif isinstance(cube_slice, iris.cube.Cube): 2303 ↛ 2306line 2303 didn't jump to line 2306 because the condition on line 2303 was always true
2304 cubes = [cube_slice]
2305 else:
2306 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2308 # Use sequence value so multiple sequences can merge.
2309 seq_coord = cube_slice[0].coord(sequence_coordinate)
2310 plot_title, plot_filename = _set_title_and_filename(
2311 seq_coord, nplot, recipe_title, filename
2312 )
2314 # Format the coordinate value in a unit appropriate way.
2315 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2317 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2318 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2318 ↛ 2319line 2318 didn't jump to line 2319 because the condition on line 2318 was never true
2319 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2321 # Do the actual plotting.
2322 plotting_func(
2323 cube_slice,
2324 coords,
2325 stamp_coordinate,
2326 plot_filename,
2327 title,
2328 series_coordinate,
2329 )
2331 plot_index.append(plot_filename)
2332 else:
2333 # Format the title and filename using plotted series coordinate
2334 nplot = 1
2335 seq_coord = coords[0]
2336 plot_title, plot_filename = _set_title_and_filename(
2337 seq_coord, nplot, recipe_title, filename
2338 )
2340 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2341 if (
2342 "station" in [c.name() for c in cubes[0].coords()]
2343 and len(cubes[0].coord("station").points) > 1
2344 ):
2345 for station in cubes[0].coord("station").points:
2346 station_cubes = cubes.extract(iris.Constraint(station=station))
2347 station_name = station_cubes[0].coord("Station_Name").points[0]
2348 station_plotname = plot_filename.replace(
2349 ".png", "_" + station_name + ".png"
2350 )
2351 _plot_and_save_line_series(
2352 station_cubes,
2353 coords,
2354 "realization",
2355 station_plotname,
2356 f"{plot_title} {station_name}",
2357 )
2358 plot_index.append(station_plotname)
2360 else:
2361 # Do the actual plotting for all other series coordinate options.
2362 _plot_and_save_line_series(
2363 cubes, coords, stamp_coordinate, plot_filename, plot_title
2364 )
2366 plot_index.append(plot_filename)
2368 # append plot to list of plots
2369 complete_plot_index = _append_to_plot_index(plot_index)
2371 # Make a page to display the plots.
2372 _make_plot_html_page(complete_plot_index)
2374 return cube
2377def plot_vertical_line_series(
2378 cubes: iris.cube.Cube | iris.cube.CubeList,
2379 filename: str | None = None,
2380 series_coordinate: str = "model_level_number",
2381 sequence_coordinate: str = "time",
2382 # line_coordinate: str = "realization",
2383 **kwargs,
2384) -> iris.cube.Cube | iris.cube.CubeList:
2385 """Plot a line plot against a type of vertical coordinate.
2387 The Cube or CubeList must be 1D.
2389 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2390 then a sequence of plots will be produced.
2392 Parameters
2393 ----------
2394 iris.cube | iris.cube.CubeList
2395 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2396 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2397 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2398 filename: str, optional
2399 Name of the plot to write, used as a prefix for plot sequences. Defaults
2400 to the recipe name.
2401 series_coordinate: str, optional
2402 Coordinate to plot on the y-axis. Can be ``pressure`` or
2403 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2404 for LFRic. Defaults to ``model_level_number``.
2405 This coordinate must exist in the cube.
2406 sequence_coordinate: str, optional
2407 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2408 This coordinate must exist in the cube.
2410 Returns
2411 -------
2412 iris.cube.Cube | iris.cube.CubeList
2413 The original Cube or CubeList (so further operations can be applied).
2414 Plotted data.
2416 Raises
2417 ------
2418 ValueError
2419 If the cubes doesn't have the right dimensions.
2420 TypeError
2421 If the cube isn't a Cube or CubeList.
2422 """
2423 # Ensure we have a name for the plot file.
2424 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2426 cubes = iter_maybe(cubes)
2427 # Initialise empty list to hold all data from all cubes in a CubeList
2428 all_data = []
2430 # Store min/max ranges for x range.
2431 x_levels = []
2433 num_models = get_num_models(cubes)
2435 validate_cube_shape(cubes, num_models)
2437 # Iterate over all cubes in cube or CubeList and plot.
2438 coords = []
2439 for cube in cubes:
2440 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2441 try:
2442 coords.append(cube.coord(series_coordinate))
2443 except iris.exceptions.CoordinateNotFoundError as err:
2444 raise ValueError(
2445 f"Cube must have a {series_coordinate} coordinate."
2446 ) from err
2448 try:
2449 if cube.ndim > 1 or not cube.coords("realization"): 2449 ↛ 2457line 2449 didn't jump to line 2457 because the condition on line 2449 was always true
2450 cube.coord(sequence_coordinate)
2451 except iris.exceptions.CoordinateNotFoundError as err:
2452 raise ValueError(
2453 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2454 ) from err
2456 # Get minimum and maximum from levels information.
2457 _, levels, _ = colorbar_map_levels(cube, axis="x")
2458 if levels is not None: 2458 ↛ 2462line 2458 didn't jump to line 2462 because the condition on line 2458 was always true
2459 x_levels.append(min(levels))
2460 x_levels.append(max(levels))
2461 else:
2462 all_data.append(cube.data)
2464 if len(x_levels) == 0: 2464 ↛ 2466line 2464 didn't jump to line 2466 because the condition on line 2464 was never true
2465 # Combine all data into a single NumPy array
2466 combined_data = np.concatenate(all_data)
2468 # Set the lower and upper limit for the x-axis to ensure all plots have
2469 # same range. This needs to read the whole cube over the range of the
2470 # sequence and if applicable postage stamp coordinate.
2471 vmin = np.floor(combined_data.min())
2472 vmax = np.ceil(combined_data.max())
2473 else:
2474 vmin = min(x_levels)
2475 vmax = max(x_levels)
2477 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2478 # a single profile directly without iterating over a sequence.
2479 sequence_coords = [
2480 cube.coord(sequence_coordinate)
2481 for cube in cubes
2482 if cube.coords(sequence_coordinate)
2483 ]
2484 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2485 np.size(coord.points) > 1 for coord in sequence_coords
2486 )
2487 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2488 np.size(coord.points) == 1 for coord in sequence_coords
2489 )
2491 plot_index = []
2492 if has_sequence_coord: 2492 ↛ 2517line 2492 didn't jump to line 2517 because the condition on line 2492 was always true
2493 # Matching the slices (matching by seq coord point; it may happen that
2494 # evaluated models do not cover the same seq coord range, hence matching
2495 # necessary)
2496 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2497 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2498 for cubes_slice in cube_iterables:
2499 # Format the coordinate value in a unit appropriate way.
2500 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2501 plot_title, plot_filename = _set_title_and_filename(
2502 seq_coord, nplot, recipe_title, filename
2503 )
2505 # Do the actual plotting.
2506 _plot_and_save_vertical_line_series(
2507 cubes_slice,
2508 coords,
2509 "realization",
2510 plot_filename,
2511 series_coordinate,
2512 title=plot_title,
2513 vmin=vmin,
2514 vmax=vmax,
2515 )
2516 plot_index.append(plot_filename)
2517 elif has_scalar_sequence_coord:
2518 # Scalar sequence coordinate (typically aggregated time bounds):
2519 # make one plot and include sequence period in title/filename.
2520 plot_title, plot_filename = _set_title_and_filename(
2521 sequence_coords[0], 1, recipe_title, filename
2522 )
2524 _plot_and_save_vertical_line_series(
2525 cubes,
2526 coords,
2527 "realization",
2528 plot_filename,
2529 series_coordinate,
2530 title=plot_title,
2531 vmin=vmin,
2532 vmax=vmax,
2533 )
2534 plot_index.append(plot_filename)
2535 else:
2536 # 1D case: no sequence coordinate, plot a single profile.
2537 plot_title = recipe_title
2538 if filename:
2539 plot_filename = filename
2540 else:
2541 plot_filename = f"{slugify(plot_title)}.png"
2543 _plot_and_save_vertical_line_series(
2544 cubes,
2545 coords,
2546 "realization",
2547 plot_filename,
2548 series_coordinate,
2549 title=plot_title,
2550 vmin=vmin,
2551 vmax=vmax,
2552 )
2553 plot_index.append(plot_filename)
2555 # Add list of plots to plot metadata.
2556 complete_plot_index = _append_to_plot_index(plot_index)
2558 # Make a page to display the plots.
2559 _make_plot_html_page(complete_plot_index)
2561 return cubes
2564def qq_plot(
2565 cubes: iris.cube.CubeList,
2566 coordinates: list[str],
2567 percentiles: list[float],
2568 model_names: list[str],
2569 filename: str | None = None,
2570 one_to_one: bool = True,
2571 **kwargs,
2572) -> iris.cube.CubeList:
2573 """Plot a Quantile-Quantile plot between two models for common time points.
2575 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2576 collapsed within the operator over all specified coordinates such as
2577 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2578 ensemble members to ensure a 1D cube (array).
2580 Parameters
2581 ----------
2582 cubes: iris.cube.CubeList
2583 Two cubes of the same variable with different models.
2584 coordinate: list[str]
2585 The list of coordinates to collapse over. This list should be
2586 every coordinate within the cube to result in a 1D cube around
2587 the percentile coordinate.
2588 percent: list[float]
2589 A list of percentiles to appear in the plot.
2590 model_names: list[str]
2591 A list of model names to appear on the axis of the plot.
2592 filename: str, optional
2593 Filename of the plot to write.
2594 one_to_one: bool, optional
2595 If True a 1:1 line is plotted; if False it is not. Default is True.
2597 Raises
2598 ------
2599 ValueError
2600 When the cubes are not compatible.
2602 Notes
2603 -----
2604 The quantile-quantile plot is a variant on the scatter plot representing
2605 two datasets by their quantiles (percentiles) for common time points.
2606 This plot does not use a theoretical distribution to compare against, but
2607 compares percentiles of two datasets. This plot does
2608 not use all raw data points, but plots the selected percentiles (quantiles) of
2609 each variable instead for the two datasets, thereby normalising the data for a
2610 direct comparison between the selected percentiles of the two dataset distributions.
2612 Quantile-quantile plots are valuable for comparing against
2613 observations and other models. Identical percentiles between the variables
2614 will lie on the one-to-one line implying the values correspond well to each
2615 other. Where there is a deviation from the one-to-one line a range of
2616 possibilities exist depending on how and where the data is shifted (e.g.,
2617 Wilks 2011 [Wilks2011]_).
2619 For distributions above the one-to-one line the distribution is left-skewed;
2620 below is right-skewed. A distinct break implies a bimodal distribution, and
2621 closer values/values further apart at the tails imply poor representation of
2622 the extremes.
2624 References
2625 ----------
2626 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2627 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2628 """
2629 # Check cubes using same functionality as the difference operator.
2630 if len(cubes) != 2:
2631 raise ValueError("cubes should contain exactly 2 cubes.")
2632 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2633 other: Cube = cubes.extract_cube(
2634 iris.Constraint(
2635 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2636 )
2637 )
2639 # Get spatial coord names.
2640 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2641 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2643 # Ensure cubes to compare are on common differencing grid.
2644 # This is triggered if either
2645 # i) latitude and longitude shapes are not the same. Note grid points
2646 # are not compared directly as these can differ through rounding
2647 # errors.
2648 # ii) or variables are known to often sit on different grid staggering
2649 # in different models (e.g. cell center vs cell edge), as is the case
2650 # for UM and LFRic comparisons.
2651 # In future greater choice of regridding method might be applied depending
2652 # on variable type. Linear regridding can in general be appropriate for smooth
2653 # variables. Care should be taken with interpretation of differences
2654 # given this dependency on regridding.
2655 if (
2656 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2657 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2658 ) or (
2659 base.long_name
2660 in [
2661 "eastward_wind_at_10m",
2662 "northward_wind_at_10m",
2663 "northward_wind_at_cell_centres",
2664 "eastward_wind_at_cell_centres",
2665 "zonal_wind_at_pressure_levels",
2666 "meridional_wind_at_pressure_levels",
2667 "potential_vorticity_at_pressure_levels",
2668 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2669 ]
2670 ):
2671 logger.debug("Linear regridding base cube to other grid to compute differences")
2672 base = regrid_onto_cube(base, other, method="Linear")
2674 # Extract just common time points.
2675 base, other = _extract_common_time_points(base, other)
2677 # Equalise attributes so we can merge.
2678 fully_equalise_attributes([base, other])
2679 logger.debug("Base: %s\nOther: %s", base, other)
2681 # Collapse cubes.
2682 base = collapse(
2683 base,
2684 coordinate=coordinates,
2685 method="PERCENTILE",
2686 additional_percent=percentiles,
2687 )
2688 other = collapse(
2689 other,
2690 coordinate=coordinates,
2691 method="PERCENTILE",
2692 additional_percent=percentiles,
2693 )
2695 # Ensure we have a name for the plot file.
2696 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2697 title = f"{recipe_title}"
2699 if filename is None:
2700 filename = slugify(recipe_title)
2702 # Add file extension.
2703 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2705 # Do the actual plotting on a scatter plot
2706 _plot_and_save_scatter_plot(
2707 base, other, plot_filename, title, one_to_one, model_names
2708 )
2710 # Add list of plots to plot metadata.
2711 plot_index = _append_to_plot_index([plot_filename])
2713 # Make a page to display the plots.
2714 _make_plot_html_page(plot_index)
2716 return iris.cube.CubeList([base, other])
2719def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2720 """
2721 Plot a Hinton style triangle/scorecard plot.
2723 This plot type can be useful for summarising high level information, such as comparing
2724 how 'skillful' two models are when verified against observations for a variety of metrics,
2725 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2726 than customisable by the user as input arguments; many have been designed to automatically
2727 scale the plot depending on the number of x and y components.
2729 Parameters
2730 ----------
2731 change: np.ndarray
2732 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2733 size/direction.
2734 signif: np.ndarray
2735 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2736 xaxis_labels: list
2737 List of labels for the xaxis (must match the second dimension length of signif and change,
2738 along with magnitude if not None).
2739 yaxis_labels: list
2740 List of labels for the yaxis (must match the first dimension length of signif and change,
2741 along with magnitude if not None).
2742 magnitude: np.ndarray | None
2743 Optional 2D array, matching the shape of change, signif, which contains numerical values
2744 the user wishes to display under each respective triangle.
2746 Returns
2747 -------
2748 matplotlib axes object to either display or do further modifications to.
2749 """
2750 # Setup colors of triangles
2751 color_pos = "#7CAE00"
2752 color_neg = "#7B68EE"
2754 # Setup cell/text size ratios
2755 figsize = None
2756 cell_size_in = 0.35
2757 text_row_ratio = 0.25
2759 # Ensure arrays, and change to bool for sig.
2760 change = np.asarray(change)
2761 signif = np.asarray(signif).astype(bool)
2762 if magnitude is not None: 2762 ↛ 2763line 2762 didn't jump to line 2763 because the condition on line 2762 was never true
2763 magnitude = np.asarray(magnitude)
2765 # Get the number of x and y elements
2766 ny, nx = change.shape
2768 # Build non-uniform y coordinates
2769 tri_height = 1.0
2770 txt_height = text_row_ratio
2772 tri_y = []
2773 txt_y = []
2774 y_edges = [0.0]
2776 y = 0.0
2777 for _j in range(ny):
2778 tri_y.append(y + tri_height / 2)
2779 y += tri_height
2780 y_edges.append(y)
2782 if magnitude is not None: 2782 ↛ 2783line 2782 didn't jump to line 2783 because the condition on line 2782 was never true
2783 txt_y.append(y + txt_height / 2)
2784 y += txt_height
2785 y_edges.append(y)
2787 total_height = y
2789 # Dynamic figure size
2790 if figsize is None: 2790 ↛ 2795line 2790 didn't jump to line 2795 because the condition on line 2790 was always true
2791 width = nx * cell_size_in
2792 height = total_height * cell_size_in + 2
2793 figsize = (width, height)
2795 fig, ax = plt.subplots(figsize=figsize)
2797 # Setup axes and grid.
2798 ax.set_aspect("equal", adjustable="box")
2799 ax.set_xlim(-0.5, nx - 0.5)
2800 ax.set_ylim(0, total_height)
2802 ax.set_xticks(np.arange(nx))
2803 ax.set_xticklabels(xaxis_labels, rotation=90)
2805 ax.set_yticks(tri_y)
2806 ax.set_yticklabels(yaxis_labels)
2808 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2809 ax.set_yticks(y_edges, minor=True)
2811 ax.set_axisbelow(True)
2812 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2813 ax.grid(False, which="major")
2814 ax.tick_params(which="minor", length=0)
2816 ax.invert_yaxis()
2818 # Compute marker scaling (fixed overlap)
2819 fig.canvas.draw()
2821 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2822 width_in, height_in = bbox.width, bbox.height
2824 cell_w = (width_in * fig.dpi) / nx
2825 cell_h = (height_in * fig.dpi) / total_height
2826 cell_pixels = min(cell_w, cell_h)
2828 max_marker_size = (0.6 * cell_pixels) ** 2
2830 text_fontsize = cell_pixels * 0.15
2832 # Plot triangles + text
2833 for j in range(ny):
2834 for i in range(nx):
2835 val = change[j, i]
2836 if np.isnan(val): 2836 ↛ 2837line 2836 didn't jump to line 2837 because the condition on line 2836 was never true
2837 continue
2839 if abs(val) < 0.01: 2839 ↛ 2840line 2839 didn't jump to line 2840 because the condition on line 2839 was never true
2840 continue
2842 sig = signif[j, i]
2843 size = max_marker_size * abs(val)
2845 # Triangle style
2846 if val >= 0:
2847 marker = "^"
2848 color = color_pos
2849 else:
2850 marker = "v"
2851 color = color_neg
2853 if sig:
2854 edgecolor = "black"
2855 linewidth = 0.6
2856 else:
2857 edgecolor = "none"
2858 linewidth = 0.0
2860 # Triangle
2861 ax.scatter(
2862 i,
2863 tri_y[j],
2864 s=size,
2865 marker=marker,
2866 c=color,
2867 edgecolors=edgecolor,
2868 linewidths=linewidth,
2869 zorder=3,
2870 clip_on=True, # ensures no rendering bleed
2871 )
2873 # Text row
2874 if magnitude is not None: 2874 ↛ 2875line 2874 didn't jump to line 2875 because the condition on line 2874 was never true
2875 mag_val = magnitude[j, i]
2877 if not np.isnan(mag_val):
2878 ax.text(
2879 i,
2880 txt_y[j],
2881 f"{mag_val:.1f}",
2882 ha="center",
2883 va="center",
2884 fontsize=text_fontsize,
2885 color="black",
2886 zorder=4,
2887 )
2889 plt.tight_layout()
2890 return fig, ax
2893def scatter_plot(
2894 cube_x: iris.cube.Cube | iris.cube.CubeList,
2895 cube_y: iris.cube.Cube | iris.cube.CubeList,
2896 filename: str | None = None,
2897 one_to_one: bool = True,
2898 **kwargs,
2899) -> iris.cube.CubeList:
2900 """Plot a scatter plot between two variables.
2902 Both cubes must be 1D.
2904 Parameters
2905 ----------
2906 cube_x: Cube | CubeList
2907 1 dimensional Cube of the data to plot on y-axis.
2908 cube_y: Cube | CubeList
2909 1 dimensional Cube of the data to plot on x-axis.
2910 filename: str, optional
2911 Filename of the plot to write.
2912 one_to_one: bool, optional
2913 If True a 1:1 line is plotted; if False it is not. Default is True.
2915 Returns
2916 -------
2917 cubes: CubeList
2918 CubeList of the original x and y cubes for further processing.
2920 Raises
2921 ------
2922 ValueError
2923 If the cube doesn't have the right dimensions and cubes not the same
2924 size.
2925 TypeError
2926 If the cube isn't a single cube.
2928 Notes
2929 -----
2930 Scatter plots are used for determining if there is a relationship between
2931 two variables. Positive relations have a slope going from bottom left to top
2932 right; Negative relations have a slope going from top left to bottom right.
2933 """
2934 # Iterate over all cubes in cube or CubeList and plot.
2935 for cube_iter in iter_maybe(cube_x):
2936 # Check cubes are correct shape.
2937 cube_iter = check_single_cube(cube_iter)
2938 if cube_iter.ndim > 1:
2939 raise ValueError("cube_x must be 1D.")
2941 # Iterate over all cubes in cube or CubeList and plot.
2942 for cube_iter in iter_maybe(cube_y):
2943 # Check cubes are correct shape.
2944 cube_iter = check_single_cube(cube_iter)
2945 if cube_iter.ndim > 1:
2946 raise ValueError("cube_y must be 1D.")
2948 # Ensure we have a name for the plot file.
2949 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2950 title = f"{recipe_title}"
2952 if filename is None:
2953 filename = slugify(recipe_title)
2955 # Add file extension.
2956 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2958 # Do the actual plotting.
2959 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2961 # Add list of plots to plot metadata.
2962 plot_index = _append_to_plot_index([plot_filename])
2964 # Make a page to display the plots.
2965 _make_plot_html_page(plot_index)
2967 return iris.cube.CubeList([cube_x, cube_y])
2970def vector_plot(
2971 cube_u: iris.cube.Cube,
2972 cube_v: iris.cube.Cube,
2973 filename: str | None = None,
2974 sequence_coordinate: str = "time",
2975 **kwargs,
2976) -> iris.cube.CubeList:
2977 """Plot a vector plot based on the input u and v components."""
2978 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2980 # Cubes must have a matching sequence coordinate.
2981 try:
2982 # Check that the u and v cubes have the same sequence coordinate.
2983 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2983 ↛ anywhereline 2983 didn't jump anywhere: it always raised an exception.
2984 raise ValueError("Coordinates do not match.")
2985 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2986 raise ValueError(
2987 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2988 ) from err
2990 # Create a plot for each value of the sequence coordinate.
2991 plot_index = []
2992 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2993 for cube_u_slice, cube_v_slice in zip(
2994 cube_u.slices_over(sequence_coordinate),
2995 cube_v.slices_over(sequence_coordinate),
2996 strict=True,
2997 ):
2998 # Format the coordinate value in a unit appropriate way.
2999 seq_coord = cube_u_slice.coord(sequence_coordinate)
3000 plot_title, plot_filename = _set_title_and_filename(
3001 seq_coord, nplot, recipe_title, filename
3002 )
3004 # Do the actual plotting.
3005 _plot_and_save_vector_plot(
3006 cube_u_slice,
3007 cube_v_slice,
3008 filename=plot_filename,
3009 title=plot_title,
3010 method="pcolormesh",
3011 )
3012 plot_index.append(plot_filename)
3014 # Add list of plots to plot metadata.
3015 complete_plot_index = _append_to_plot_index(plot_index)
3017 # Make a page to display the plots.
3018 _make_plot_html_page(complete_plot_index)
3020 return iris.cube.CubeList([cube_u, cube_v])
3023def plot_histogram_series(
3024 cubes: iris.cube.Cube | iris.cube.CubeList,
3025 filename: str | None = None,
3026 sequence_coordinate: str = "time",
3027 stamp_coordinate: str = "realization",
3028 single_plot: bool = False,
3029 **kwargs,
3030) -> iris.cube.Cube | iris.cube.CubeList:
3031 """Plot a histogram plot for each vertical level provided.
3033 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3034 is present then a sequence of plots will be produced using the time slider
3035 functionality to scroll through histograms against time. If a
3036 stamp_coordinate is present then postage stamp plots will be produced. If
3037 stamp_coordinate and single_plot is True, all postage stamp plots will be
3038 plotted in a single plot instead of separate postage stamp plots.
3040 Parameters
3041 ----------
3042 cubes: Cube | iris.cube.CubeList
3043 Iris cube or CubeList of the data to plot. It should have a single dimension other
3044 than the stamp coordinate.
3045 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3046 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3047 filename: str, optional
3048 Name of the plot to write, used as a prefix for plot sequences. Defaults
3049 to the recipe name.
3050 sequence_coordinate: str, optional
3051 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3052 This coordinate must exist in the cube and will be used for the time
3053 slider.
3054 stamp_coordinate: str, optional
3055 Coordinate about which to plot postage stamp plots. Defaults to
3056 ``"realization"``.
3057 single_plot: bool, optional
3058 If True, all postage stamp plots will be plotted in a single plot. If
3059 False, each postage stamp plot will be plotted separately. Is only valid
3060 if stamp_coordinate exists and has more than a single point.
3062 Returns
3063 -------
3064 iris.cube.Cube | iris.cube.CubeList
3065 The original Cube or CubeList (so further operations can be applied).
3066 Plotted data.
3068 Raises
3069 ------
3070 ValueError
3071 If the cube doesn't have the right dimensions.
3072 TypeError
3073 If the cube isn't a Cube or CubeList.
3074 """
3075 recipe_title = get_recipe_metadata().get("title", "Histogram")
3077 cubes = iter_maybe(cubes)
3079 # Internal plotting function.
3080 plotting_func = _plot_and_save_histogram_series
3082 num_models = get_num_models(cubes)
3084 validate_cube_shape(cubes, num_models)
3086 # If several histograms are plotted, check sequence_coordinate
3087 check_sequence_coordinate(cubes, sequence_coordinate)
3089 # Get axis minimum and maximum from levels information.
3090 # If no levels set, derive minima and maxima from data in CubeList.
3091 vmin, vmax = _set_axis_range(cubes)
3093 # Make postage stamp plots if stamp_coordinate exists and has more than a
3094 # single point. If single_plot is True:
3095 # -- all postage stamp plots will be plotted in a single plot instead of
3096 # separate postage stamp plots.
3097 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3098 # produced per single model only
3099 if num_models == 1:
3100 if ( 3100 ↛ 3104line 3100 didn't jump to line 3104 because the condition on line 3100 was never true
3101 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3102 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3103 ):
3104 if single_plot:
3105 plotting_func = (
3106 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3107 )
3108 else:
3109 plotting_func = _plot_and_save_postage_stamp_histogram_series
3110 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3111 else:
3112 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3114 plot_index = []
3115 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3116 # Create a plot for each value of the sequence coordinate. Allowing for
3117 # multiple cubes in a CubeList to be plotted in the same plot for similar
3118 # sequence values. Passing a CubeList into the internal plotting function
3119 # for similar values of the sequence coordinate. cube_slice can be an
3120 # iris.cube.Cube or an iris.cube.CubeList.
3121 for cube_slice in cube_iterables:
3122 single_cube = cube_slice
3123 if isinstance(cube_slice, iris.cube.CubeList):
3124 single_cube = cube_slice[0]
3126 # Ensure valid stamp coordinate in cube dimensions
3127 if stamp_coordinate == "realization": 3127 ↛ 3130line 3127 didn't jump to line 3130 because the condition on line 3127 was always true
3128 stamp_coordinate = check_stamp_coordinate(single_cube)
3129 # Set plot titles and filename, based on sequence coordinate
3130 seq_coord = single_cube.coord(sequence_coordinate)
3131 # Use time coordinate in title and filename if single histogram output.
3132 if sequence_coordinate == "realization" and nplot == 1: 3132 ↛ 3133line 3132 didn't jump to line 3133 because the condition on line 3132 was never true
3133 seq_coord = single_cube.coord("time")
3134 # Use station name in title and filename if model vs obs comparison
3135 if sequence_coordinate == "station": 3135 ↛ 3136line 3135 didn't jump to line 3136 because the condition on line 3135 was never true
3136 seq_coord = single_cube.coord("Station_Name")
3138 plot_title, plot_filename = _set_title_and_filename(
3139 seq_coord, nplot, recipe_title, filename
3140 )
3142 # Do the actual plotting.
3143 plotting_func(
3144 cube_slice,
3145 filename=plot_filename,
3146 stamp_coordinate=stamp_coordinate,
3147 title=plot_title,
3148 vmin=vmin,
3149 vmax=vmax,
3150 )
3151 plot_index.append(plot_filename)
3153 # Add list of plots to plot metadata.
3154 complete_plot_index = _append_to_plot_index(plot_index)
3156 # Make a page to display the plots.
3157 _make_plot_html_page(complete_plot_index)
3159 return cubes
3162def plot_scatter_series(
3163 cubes: iris.cube.Cube | iris.cube.CubeList,
3164 filename: str | None = None,
3165 sequence_coordinate: str = "time",
3166 stamp_coordinate: str = "realization",
3167 hexbin: bool = False,
3168 **kwargs,
3169) -> iris.cube.Cube | iris.cube.CubeList:
3170 """Plot a scatter plot for each sequence coordinate provided.
3172 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3173 is present then a sequence of plots will be produced using the time slider
3174 functionality to scroll through scatter against time. If a
3175 stamp_coordinate is present then postage stamp plots will be produced. If
3176 stamp_coordinate and single_plot is True, all postage stamp plots will be
3177 plotted in a single plot instead of separate postage stamp plots.
3179 Parameters
3180 ----------
3181 cubes: Cube | iris.cube.CubeList
3182 Iris cube or CubeList of the data to plot. It should have a single dimension other
3183 than the stamp coordinate.
3184 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3185 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3186 filename: str, optional
3187 Name of the plot to write, used as a prefix for plot sequences. Defaults
3188 to the recipe name.
3189 sequence_coordinate: str, optional
3190 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3191 This coordinate must exist in the cube and will be used for the time
3192 slider.
3193 stamp_coordinate: str, optional
3194 Coordinate about which to plot postage stamp plots. Defaults to
3195 ``"realization"``.
3196 hexbin: bool, optional
3197 If True, generate hexbin comparison plot.
3198 If False, generate point-by-point scatter plot.
3200 Returns
3201 -------
3202 iris.cube.Cube | iris.cube.CubeList
3203 The original Cube or CubeList (so further operations can be applied).
3204 Plotted data.
3206 Raises
3207 ------
3208 ValueError
3209 If the cube doesn't have the right dimensions.
3210 TypeError
3211 If the cube isn't a Cube or CubeList.
3212 """
3213 recipe_title = get_recipe_metadata().get("title", "Scatter")
3215 cubes = iter_maybe(cubes)
3217 # Internal plotting function.
3218 plotting_func = _plot_and_save_scatter_series
3220 num_models = get_num_models(cubes)
3222 validate_cube_shape(cubes, num_models)
3224 check_sequence_coordinate(cubes, sequence_coordinate)
3226 vmin, vmax = _set_axis_range(cubes)
3228 # Require >1 models to compare on scatter plot
3229 if num_models > 1:
3230 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3231 else:
3232 raise ValueError(
3233 "Scatter plot series requires multiple number of models in input data."
3234 )
3236 plot_index = []
3237 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3238 # Create a plot for each value of the sequence coordinate. Allowing for
3239 # multiple cubes in a CubeList to be plotted in the same plot for similar
3240 # sequence values. Passing a CubeList into the internal plotting function
3241 # for similar values of the sequence coordinate. cube_slice can be an
3242 # iris.cube.Cube or an iris.cube.CubeList.
3243 for cube_slice in cube_iterables:
3244 single_cube = cube_slice
3245 if isinstance(cube_slice, iris.cube.CubeList): 3245 ↛ 3249line 3245 didn't jump to line 3249 because the condition on line 3245 was always true
3246 single_cube = cube_slice[0]
3248 # Ensure valid stamp coordinate in cube dimensions
3249 if stamp_coordinate == "realization": 3249 ↛ 3252line 3249 didn't jump to line 3252 because the condition on line 3249 was always true
3250 stamp_coordinate = check_stamp_coordinate(single_cube)
3251 # Set plot titles and filename, based on sequence coordinate
3252 seq_coord = single_cube.coord(sequence_coordinate)
3253 # Use time coordinate in title and filename if single histogram output.
3254 if sequence_coordinate == "realization" and nplot == 1:
3255 seq_coord = single_cube.coord("time")
3256 # Use station name in title and filename if model vs obs comparison
3257 if sequence_coordinate == "station":
3258 seq_coord = single_cube.coord("Station_Name")
3260 plot_title, plot_filename = _set_title_and_filename(
3261 seq_coord, nplot, recipe_title, filename
3262 )
3264 # Do the actual plotting.
3265 plotting_func(
3266 cube_slice,
3267 filename=plot_filename,
3268 stamp_coordinate=stamp_coordinate,
3269 title=plot_title,
3270 vmin=vmin,
3271 vmax=vmax,
3272 hexbin=hexbin,
3273 )
3274 plot_index.append(plot_filename)
3276 # Add list of plots to plot metadata.
3277 complete_plot_index = _append_to_plot_index(plot_index)
3279 # Make a page to display the plots.
3280 _make_plot_html_page(complete_plot_index)
3282 return cubes
3285def _plot_and_save_postage_stamp_power_spectrum_series(
3286 cubes: iris.cube.Cube,
3287 coords: list[iris.coords.Coord],
3288 stamp_coordinate: str,
3289 filename: str,
3290 title: str,
3291 series_coordinate: str | None = None,
3292 **kwargs,
3293):
3294 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3296 Parameters
3297 ----------
3298 cubes: Cube or CubeList
3299 Cube or Cubelist of the power spectrum data.
3300 coords: list[Coord]
3301 Coordinates to plot on the x-axis, one per cube.
3302 stamp_coordinate: str
3303 Coordinate that becomes different plots.
3304 filename: str
3305 Filename of the plot to write.
3306 title: str
3307 Plot title.
3308 series_coordinate: str, optional
3309 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3311 """
3312 # Use the smallest square grid that will fit the members.
3313 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3315 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3316 model_colors_map = get_model_colors_map(cubes)
3317 # ax = plt.gca()
3318 # Make a subplot for each member.
3319 for member, subplot in zip(
3320 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3321 ):
3322 ax = plt.subplot(grid_size, grid_size, subplot)
3324 # Store min/max ranges.
3325 y_levels = []
3327 line_marker = None
3328 line_width = 1
3330 for cube in iter_maybe(member):
3331 xcoord = _select_series_coord(cube, series_coordinate)
3332 xname = xcoord.points
3334 yfield = cube.data # power spectrum
3335 label = None
3336 color = "black"
3337 if model_colors_map: 3337 ↛ 3338line 3337 didn't jump to line 3338 because the condition on line 3337 was never true
3338 label = cube.attributes.get("model_name")
3339 color = model_colors_map.get(label)
3341 if member.coord(stamp_coordinate).points == [0]:
3342 ax.plot(
3343 xname,
3344 yfield,
3345 color=color,
3346 marker=line_marker,
3347 ls="-",
3348 lw=line_width,
3349 label=f"{label} (control)"
3350 if len(cube.coord(stamp_coordinate).points) > 1
3351 else label,
3352 )
3353 # Label with member if part of an ensemble and not the control.
3354 else:
3355 ax.plot(
3356 xname,
3357 yfield,
3358 color=color,
3359 ls="-",
3360 lw=1.5,
3361 alpha=0.75,
3362 label=f"{label} (member)",
3363 )
3365 # Calculate the global min/max if multiple cubes are given.
3366 _, levels, _ = colorbar_map_levels(cube, axis="y")
3367 if levels is not None: 3367 ↛ 3368line 3367 didn't jump to line 3368 because the condition on line 3367 was never true
3368 y_levels.append(min(levels))
3369 y_levels.append(max(levels))
3371 # Add some labels and tweak the style.
3372 title = f"{title}"
3373 ax.set_title(title, fontsize=16)
3375 # Set appropriate x-axis label based on coordinate
3376 if series_coordinate == "wavelength" or ( 3376 ↛ 3379line 3376 didn't jump to line 3379 because the condition on line 3376 was never true
3377 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3378 ):
3379 ax.set_xlabel("Wavelength (km)", fontsize=14)
3380 elif series_coordinate == "physical_wavenumber" or ( 3380 ↛ 3385line 3380 didn't jump to line 3385 because the condition on line 3380 was always true
3381 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3382 ):
3383 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3384 else: # frequency or check units
3385 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3386 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3387 else:
3388 ax.set_xlabel("Wavenumber", fontsize=14)
3390 ax.set_ylabel("Power Spectral Density", fontsize=14)
3391 ax.tick_params(axis="both", labelsize=12)
3393 # Set log-log scale
3394 ax.set_xscale("log")
3395 ax.set_yscale("log")
3397 # Add gridlines
3398 ax.grid(linestyle="--", color="grey", linewidth=1)
3399 # Ientify unique labels for legend
3400 handles = list(
3401 {
3402 label: handle
3403 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3404 }.values()
3405 )
3406 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3408 ax = plt.gca()
3409 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3411 # Save plot.
3412 _save_close_figure(fig, "histogram postage stamp", filename)
3415def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3416 cubes: iris.cube.Cube,
3417 coords: list[iris.coords.Coord],
3418 stamp_coordinate: str,
3419 filename: str,
3420 title: str,
3421 series_coordinate: str | None = None,
3422 **kwargs,
3423):
3424 """Plot and save power spectra for ensemble members in single plot.
3426 Parameters
3427 ----------
3428 cubes: Cube or CubeList
3429 Cube or Cubelist of the power spectrum data.
3430 coords: list[Coord]
3431 Coordinates to plot on the x-axis, one per cube.
3432 stamp_coordinate: str
3433 Coordinate that becomes different plots.
3434 filename: str
3435 Filename of the plot to write.
3436 title: str
3437 Plot title.
3438 series_coordinate: str, optional
3439 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3441 """
3442 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3443 model_colors_map = get_model_colors_map(cubes)
3445 line_marker = None
3446 line_width = 1
3448 # Compute ensemble statistics to show spread
3449 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3450 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3451 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3453 xcoord_global = mean_cube.coord(series_coordinate)
3454 x_global = xcoord_global.points
3456 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3457 xcoord = _select_series_coord(member, series_coordinate)
3458 xname = xcoord.points
3460 yfield = member.data # power spectrum
3461 color = "black"
3462 if model_colors_map: 3462 ↛ 3466line 3462 didn't jump to line 3466 because the condition on line 3462 was always true
3463 label = member.attributes.get("model_name") if i == 0 else None
3464 color = model_colors_map.get(label)
3466 if member.coord(stamp_coordinate).points == [0]:
3467 ax.plot(
3468 xname,
3469 yfield,
3470 color=color,
3471 marker=line_marker,
3472 ls="-",
3473 lw=line_width,
3474 label=f"{label} (control)"
3475 if len(member.coord(stamp_coordinate).points) > 1
3476 else label,
3477 )
3478 # Label with member number if part of an ensemble and not the control.
3479 else:
3480 ax.plot(
3481 xname,
3482 yfield,
3483 color=color,
3484 ls="-",
3485 lw=1.5,
3486 alpha=0.75,
3487 label=label,
3488 )
3490 # Set appropriate x-axis label based on coordinate
3491 if series_coordinate == "wavelength" or ( 3491 ↛ 3494line 3491 didn't jump to line 3494 because the condition on line 3491 was never true
3492 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3493 ):
3494 ax.set_xlabel("Wavelength (km)", fontsize=14)
3495 elif series_coordinate == "physical_wavenumber" or ( 3495 ↛ 3500line 3495 didn't jump to line 3500 because the condition on line 3495 was always true
3496 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3497 ):
3498 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3499 else: # frequency or check units
3500 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3501 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3502 else:
3503 ax.set_xlabel("Wavenumber", fontsize=14)
3505 # Add ensemble spread shading
3506 ax.fill_between(
3507 x_global,
3508 min_cube.data,
3509 max_cube.data,
3510 color="grey",
3511 alpha=0.3,
3512 label="Ensemble spread",
3513 )
3515 # Add ensemble mean line
3516 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3518 ax.set_ylabel("Power Spectral Density", fontsize=14)
3519 ax.tick_params(axis="both", labelsize=12)
3521 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3522 # Set log-log scale
3523 ax.set_xscale("log")
3524 ax.set_yscale("log")
3526 # Add gridlines
3527 ax.grid(linestyle="--", color="grey", linewidth=1)
3528 # Identify unique labels for legend
3529 handles = list(
3530 {
3531 label: handle
3532 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3533 }.values()
3534 )
3535 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3537 # Figure title.
3538 ax.set_title(title, fontsize=16)
3540 # Save plot.
3541 _save_close_figure(fig, "power spectra postage stamp", filename)