Coverage for src/CSET/operators/plot.py: 83%
1102 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-10 09:05 +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
1097 # If data from power spectra is all np.nans (like T+0h rainfall field which
1098 # might be full of zeros), then set yfield to zeros so it doesn't crash the
1099 # plotting.
1100 if np.all(np.isnan(yfield)):
1101 yfield = np.zeros_like(yfield)
1103 label = None
1104 color = "black"
1105 if model_colors_map: 1105 ↛ 1108line 1105 didn't jump to line 1108 because the condition on line 1105 was always true
1106 label = cube.attributes.get("model_name")
1107 color = model_colors_map.get(label)
1108 for cube_slice in cube.slices_over(ensemble_coord):
1109 # Label with (control) if part of an ensemble or not otherwise.
1110 if cube_slice.coord(ensemble_coord).points == [0]: 1110 ↛ 1124line 1110 didn't jump to line 1124 because the condition on line 1110 was always true
1111 ax.plot(
1112 xname,
1113 yfield,
1114 color=color,
1115 marker=line_marker,
1116 ls="-",
1117 lw=line_width,
1118 label=f"{label} (control)"
1119 if len(cube.coord(ensemble_coord).points) > 1
1120 else label,
1121 )
1122 # Label with (perturbed) if part of an ensemble and not the control.
1123 else:
1124 ax.plot(
1125 xname,
1126 yfield,
1127 color=color,
1128 ls="-",
1129 lw=1.5,
1130 alpha=0.75,
1131 label=f"{label} (member)",
1132 )
1134 # Calculate the global min/max if multiple cubes are given.
1135 _, levels, _ = colorbar_map_levels(cube, axis="y")
1136 if levels is not None: 1136 ↛ 1137line 1136 didn't jump to line 1137 because the condition on line 1136 was never true
1137 y_levels.append(min(levels))
1138 y_levels.append(max(levels))
1140 # Add some labels and tweak the style.
1142 title = f"{title}"
1143 ax.set_title(title, fontsize=16)
1145 # Set appropriate x-axis label based on coordinate
1146 if series_coordinate == "wavelength" or ( 1146 ↛ 1149line 1146 didn't jump to line 1149 because the condition on line 1146 was never true
1147 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1148 ):
1149 ax.set_xlabel("Wavelength (km)", fontsize=14)
1150 elif series_coordinate == "physical_wavenumber" or ( 1150 ↛ 1153line 1150 didn't jump to line 1153 because the condition on line 1150 was never true
1151 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1152 ):
1153 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1154 else: # frequency or check units
1155 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1155 ↛ 1156line 1155 didn't jump to line 1156 because the condition on line 1155 was never true
1156 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1157 else:
1158 ax.set_xlabel("Wavenumber", fontsize=14)
1160 ax.set_ylabel("Power Spectral Density", fontsize=14)
1161 ax.tick_params(axis="both", labelsize=12)
1163 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1165 # Set log-log scale
1166 ax.set_xscale("log")
1167 ax.set_yscale("log")
1169 # Add gridlines
1170 ax.grid(linestyle="--", color="grey", linewidth=1)
1171 # Ientify unique labels for legend
1172 handles = list(
1173 {
1174 label: handle
1175 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1176 }.values()
1177 )
1178 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1180 # Save plot.
1181 _save_close_figure(fig, "line power spectrum", filename)
1184def _plot_and_save_vertical_line_series(
1185 cubes: iris.cube.CubeList,
1186 coords: list[iris.coords.Coord],
1187 ensemble_coord: str,
1188 filename: str,
1189 series_coordinate: str,
1190 title: str,
1191 vmin: float,
1192 vmax: float,
1193 **kwargs,
1194):
1195 """Plot and save a 1D line series in vertical.
1197 Parameters
1198 ----------
1199 cubes: CubeList
1200 1 dimensional Cube or CubeList of the data to plot on x-axis.
1201 coord: list[Coord]
1202 Coordinates to plot on the y-axis, one per cube.
1203 ensemble_coord: str
1204 Ensemble coordinate in the cube.
1205 filename: str
1206 Filename of the plot to write.
1207 series_coordinate: str
1208 Coordinate to use as vertical axis.
1209 title: str
1210 Plot title.
1211 vmin: float
1212 Minimum value for the x-axis.
1213 vmax: float
1214 Maximum value for the x-axis.
1215 """
1216 # plot the vertical pressure axis using log scale
1217 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1219 model_colors_map = get_model_colors_map(cubes)
1221 # Check match-up across sequence coords gives consistent sizes
1222 validate_cubes_coords(cubes, coords)
1224 for cube, coord in zip(cubes, coords, strict=True):
1225 label = None
1226 color = "black"
1227 if model_colors_map: 1227 ↛ 1228line 1227 didn't jump to line 1228 because the condition on line 1227 was never true
1228 label = cube.attributes.get("model_name")
1229 color = model_colors_map.get(label)
1231 for cube_slice in cube.slices_over(ensemble_coord):
1232 # If ensemble data given plot control member with (control)
1233 # unless single forecast.
1234 if cube_slice.coord(ensemble_coord).points == [0]:
1235 iplt.plot(
1236 cube_slice,
1237 coord,
1238 color=color,
1239 marker="o",
1240 ls="-",
1241 lw=3,
1242 label=f"{label} (control)"
1243 if len(cube.coord(ensemble_coord).points) > 1
1244 else label,
1245 )
1246 # If ensemble data given plot perturbed members with (perturbed).
1247 else:
1248 iplt.plot(
1249 cube_slice,
1250 coord,
1251 color=color,
1252 ls="-",
1253 lw=1.5,
1254 alpha=0.75,
1255 label=f"{label} (member)",
1256 )
1258 # Get the current axis
1259 ax = plt.gca()
1261 # Special handling for pressure level data.
1262 if series_coordinate == "pressure": 1262 ↛ 1284line 1262 didn't jump to line 1284 because the condition on line 1262 was always true
1263 # Invert y-axis and set to log scale.
1264 ax.invert_yaxis()
1265 ax.set_yscale("log")
1267 # Define y-ticks and labels for pressure log axis.
1268 y_tick_labels = [
1269 "1000",
1270 "850",
1271 "700",
1272 "500",
1273 "300",
1274 "200",
1275 "100",
1276 ]
1277 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1279 # Set y-axis limits and ticks.
1280 ax.set_ylim(1100, 100)
1282 # Test if series_coordinate is model level data. The UM data uses
1283 # model_level_number and lfric uses full_levels as coordinate.
1284 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1285 # Define y-ticks and labels for vertical axis.
1286 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1287 y_tick_labels = [str(int(i)) for i in y_ticks]
1288 ax.set_ylim(min(y_ticks), max(y_ticks))
1290 ax.set_yticks(y_ticks)
1291 ax.set_yticklabels(y_tick_labels)
1293 # Set x-axis limits.
1294 ax.set_xlim(vmin, vmax)
1295 # Mark y=0 if present in plot.
1296 if vmin < 0.0 and vmax > 0.0: 1296 ↛ 1297line 1296 didn't jump to line 1297 because the condition on line 1296 was never true
1297 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1299 # Add some labels and tweak the style.
1300 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1301 ax.set_xlabel(
1302 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1303 )
1304 ax.set_title(title, fontsize=16)
1305 ax.ticklabel_format(axis="x")
1306 ax.tick_params(axis="y")
1307 ax.tick_params(axis="both", labelsize=12)
1309 # Add gridlines
1310 ax.grid(linestyle="--", color="grey", linewidth=1)
1311 # Ientify unique labels for legend
1312 handles = list(
1313 {
1314 label: handle
1315 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1316 }.values()
1317 )
1318 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1320 # Save plot.
1321 _save_close_figure(fig, "vertical line", filename)
1324def _plot_and_save_scatter_plot(
1325 cube_x: iris.cube.Cube | iris.cube.CubeList,
1326 cube_y: iris.cube.Cube | iris.cube.CubeList,
1327 filename: str,
1328 title: str,
1329 one_to_one: bool,
1330 model_names: list[str] | None = None,
1331 **kwargs,
1332):
1333 """Plot and save a 2D scatter plot.
1335 Parameters
1336 ----------
1337 cube_x: Cube | CubeList
1338 1 dimensional Cube or CubeList of the data to plot on x-axis.
1339 cube_y: Cube | CubeList
1340 1 dimensional Cube or CubeList of the data to plot on y-axis.
1341 filename: str
1342 Filename of the plot to write.
1343 title: str
1344 Plot title.
1345 one_to_one: bool
1346 Whether a 1:1 line is plotted.
1347 """
1348 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1349 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1350 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1351 # over the pairs simultaneously.
1353 # Ensure cube_x and cube_y are iterable
1354 cube_x_iterable = iter_maybe(cube_x)
1355 cube_y_iterable = iter_maybe(cube_y)
1357 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1358 iplt.scatter(cube_x_iter, cube_y_iter)
1359 if one_to_one is True:
1360 plt.plot(
1361 [
1362 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1363 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1364 ],
1365 [
1366 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1367 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1368 ],
1369 "k",
1370 linestyle="--",
1371 )
1372 ax = plt.gca()
1374 # Add some labels and tweak the style.
1375 if model_names is None:
1376 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1377 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1378 else:
1379 # Add the model names, these should be order of base (x) and other (y).
1380 ax.set_xlabel(
1381 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1382 )
1383 ax.set_ylabel(
1384 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1385 )
1386 ax.set_title(title, fontsize=16)
1387 ax.ticklabel_format(axis="y", useOffset=False)
1388 ax.tick_params(axis="x", labelrotation=15)
1389 ax.tick_params(axis="both", labelsize=12)
1390 ax.autoscale()
1392 # Save plot.
1393 _save_close_figure(fig, "scatter", filename)
1396def _plot_and_save_vector_plot(
1397 cube_u: iris.cube.Cube,
1398 cube_v: iris.cube.Cube,
1399 filename: str,
1400 title: str,
1401 method: Literal["contourf", "pcolormesh"],
1402 **kwargs,
1403):
1404 """Plot and save a 2D vector plot.
1406 Parameters
1407 ----------
1408 cube_u: Cube
1409 2 dimensional Cube of u component of the data.
1410 cube_v: Cube
1411 2 dimensional Cube of v component of the data.
1412 filename: str
1413 Filename of the plot to write.
1414 title: str
1415 Plot title.
1416 """
1417 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1418 # Create a cube containing the magnitude of the vector field.
1419 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1420 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1421 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1422 cube_vec_mag.rename(
1423 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1424 )
1426 # Specify the color bar
1427 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1429 # Setup plot map projection, extent and coastlines and borderlines.
1430 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1432 if method == "contourf":
1433 # Filled contour plot of the field.
1434 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1435 elif method == "pcolormesh":
1436 try:
1437 vmin = min(levels)
1438 vmax = max(levels)
1439 except TypeError:
1440 vmin, vmax = None, None
1441 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1442 # if levels are defined.
1443 if norm is not None:
1444 vmin = None
1445 vmax = None
1446 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1447 else:
1448 raise ValueError(f"Unknown plotting method: {method}")
1450 # Check to see if transect, and if so, adjust y axis.
1451 if is_transect(cube_vec_mag):
1452 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1453 axes.invert_yaxis()
1454 axes.set_yscale("log")
1455 axes.set_ylim(1100, 100)
1456 # If both model_level_number and level_height exists, iplt can construct
1457 # plot as a function of height above orography (NOT sea level).
1458 elif {"model_level_number", "level_height"}.issubset(
1459 {coord.name() for coord in cube_vec_mag.coords()}
1460 ):
1461 axes.set_yscale("log")
1463 axes.set_title(
1464 f"{title}\n"
1465 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1466 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1467 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1468 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1469 fontsize=16,
1470 )
1472 else:
1473 # Add title.
1474 axes.set_title(title, fontsize=16)
1476 # Add watermark with min/max/mean. Currently not user togglable.
1477 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1478 axes.annotate(
1479 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}",
1480 xy=(0.05, -0.05),
1481 xycoords="axes fraction",
1482 xytext=(-5, 5),
1483 textcoords="offset points",
1484 ha="right",
1485 va="bottom",
1486 size=11,
1487 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1488 )
1490 # Add colour bar.
1491 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1492 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1493 # add ticks and tick_labels for every levels if less than 20 levels exist
1494 if levels is not None and len(levels) < 20:
1495 cbar.set_ticks(levels)
1496 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1498 # 30 barbs along the longest axis of the plot, or a barb per point for data
1499 # with less than 30 points.
1500 step = max(max(cube_u.shape) // 30, 1)
1501 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1503 # Save plot.
1504 _save_close_figure(fig, "vector", filename)
1507def _plot_and_save_histogram_series(
1508 cubes: iris.cube.Cube | iris.cube.CubeList,
1509 filename: str,
1510 title: str,
1511 vmin: float,
1512 vmax: float,
1513 **kwargs,
1514):
1515 """Plot and save a histogram series.
1517 Parameters
1518 ----------
1519 cubes: Cube or CubeList
1520 2 dimensional Cube or CubeList of the data to plot as histogram.
1521 filename: str
1522 Filename of the plot to write.
1523 title: str
1524 Plot title.
1525 vmin: float
1526 minimum for colorbar
1527 vmax: float
1528 maximum for colorbar
1529 """
1530 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1531 ax = plt.gca()
1533 model_colors_map = get_model_colors_map(cubes)
1535 # Set default that histograms will produce probability density function
1536 # at each bin (integral over range sums to 1).
1537 density = True
1539 for cube in iter_maybe(cubes):
1540 # Easier to check title (where var name originates)
1541 # than seeing if long names exist etc.
1542 # Exception case, where distribution better fits log scales/bins.
1543 if (
1544 ("surface_microphysical" in title)
1545 or ("rain accumulation" in title)
1546 or ("Rainfall rate Composite" in title)
1547 or ("Nimrod_5min" in title)
1548 ):
1549 if "amount" in title:
1550 # Compute histogram following Klingaman et al. (2017): ASoP
1551 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1552 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1553 density = False
1554 else:
1555 bins = 10.0 ** (
1556 np.arange(-10, 27, 1) / 10.0
1557 ) # Suggestion from RMED toolbox.
1558 bins = np.insert(bins, 0, 0)
1559 ax.set_yscale("log")
1560 vmin = bins[1]
1561 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1562 ax.set_xscale("log")
1563 elif "lightning" in title:
1564 bins = [0, 1, 2, 3, 4, 5]
1565 else:
1566 bins = np.linspace(vmin, vmax, 51)
1567 logger.debug(
1568 "Plotting histogram with %s bins %s - %s.",
1569 np.size(bins),
1570 np.min(bins),
1571 np.max(bins),
1572 )
1574 # Reshape cube data into a single array to allow for a single histogram.
1575 # Otherwise we plot xdim histograms stacked.
1576 cube_data_1d = (cube.data).flatten()
1578 label = None
1579 color = "black"
1580 if model_colors_map:
1581 label = cube.attributes.get("model_name")
1582 color = model_colors_map[label]
1583 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1585 # Compute area under curve.
1586 if (
1587 ("surface_microphysical" in title and "amount" in title)
1588 or ("rain_accumulation" in title)
1589 or ("Rainfall rate Composite" in title)
1590 or ("Nimrod_5min" in title)
1591 ):
1592 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1593 x = x * bin_mean / x.sum()
1594 x = x[1:]
1595 y = y[1:]
1597 ax.plot(
1598 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1599 )
1601 # Add some labels and tweak the style.
1602 ax.set_title(title, fontsize=16)
1603 ax.set_xlabel(
1604 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1605 )
1606 ax.set_ylabel("Normalised probability density", fontsize=14)
1607 if (
1608 ("surface_microphysical" in title and "amount" in title)
1609 or ("rain accumulation" in title)
1610 or ("Nimrod_5min" in title)
1611 ):
1612 ax.set_ylabel(
1613 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1614 )
1615 ax.set_xlim(vmin, vmax)
1616 ax.tick_params(axis="both", labelsize=12)
1618 # Overlay grid-lines onto histogram plot.
1619 ax.grid(linestyle="--", color="grey", linewidth=1)
1620 if model_colors_map:
1621 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1623 # Save plot.
1624 _save_close_figure(fig, "histogram", filename)
1627def _plot_and_save_postage_stamp_histogram_series(
1628 cube: iris.cube.Cube,
1629 filename: str,
1630 title: str,
1631 stamp_coordinate: str,
1632 vmin: float,
1633 vmax: float,
1634 **kwargs,
1635):
1636 """Plot and save postage (ensemble members) stamps for a histogram series.
1638 Parameters
1639 ----------
1640 cube: Cube
1641 2 dimensional Cube of the data to plot as histogram.
1642 filename: str
1643 Filename of the plot to write.
1644 title: str
1645 Plot title.
1646 stamp_coordinate: str
1647 Coordinate that becomes different plots.
1648 vmin: float
1649 minimum for pdf x-axis
1650 vmax: float
1651 maximum for pdf x-axis
1652 """
1653 # Use the smallest square grid that will fit the members.
1654 nmember = len(cube.coord(stamp_coordinate).points)
1655 grid_rows = int(math.sqrt(nmember))
1656 grid_size = math.ceil(nmember / grid_rows)
1658 fig = plt.figure(
1659 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1660 )
1661 # Make a subplot for each member.
1662 for member, subplot in zip(
1663 cube.slices_over(stamp_coordinate),
1664 range(1, grid_size * grid_rows + 1),
1665 strict=False,
1666 ):
1667 # Implicit interface is much easier here, due to needing to have the
1668 # cartopy GeoAxes generated.
1669 plt.subplot(grid_rows, grid_size, subplot)
1670 # Reshape cube data into a single array to allow for a single histogram.
1671 # Otherwise we plot xdim histograms stacked.
1672 member_data_1d = (member.data).flatten()
1673 plt.hist(member_data_1d, density=True, stacked=True)
1674 axes = plt.gca()
1675 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1676 axes.set_title(f"{mtitle}")
1677 axes.set_xlim(vmin, vmax)
1679 # Overall figure title.
1680 fig.suptitle(title, fontsize=16)
1682 # Save plot.
1683 _save_close_figure(fig, "histogram postage stamp", filename)
1686def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1687 cube: iris.cube.Cube,
1688 filename: str,
1689 title: str,
1690 stamp_coordinate: str,
1691 vmin: float,
1692 vmax: float,
1693 **kwargs,
1694):
1695 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1696 ax.set_title(title, fontsize=16)
1697 ax.set_xlim(vmin, vmax)
1698 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1699 ax.set_ylabel("normalised probability density", fontsize=14)
1700 # Loop over all slices along the stamp_coordinate
1701 for member in cube.slices_over(stamp_coordinate):
1702 # Flatten the member data to 1D
1703 member_data_1d = member.data.flatten()
1704 # Plot the histogram using plt.hist
1705 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1706 plt.hist(
1707 member_data_1d,
1708 density=True,
1709 stacked=True,
1710 label=f"{mtitle}",
1711 )
1713 # Add a legend
1714 ax.legend(fontsize=16)
1716 # Save plot.
1717 _save_close_figure(fig, "histogram postage stamp", filename)
1720def _plot_and_save_scatter_series(
1721 cubes: iris.cube.Cube | iris.cube.CubeList,
1722 filename: str,
1723 title: str,
1724 vmin: float,
1725 vmax: float,
1726 hexbin: bool,
1727 **kwargs,
1728):
1729 """Plot and save a scatter plot series.
1731 Parameters
1732 ----------
1733 cubes: Cube or CubeList
1734 2 dimensional Cube or CubeList of the data to plot as scatter.
1735 filename: str
1736 Filename of the plot to write.
1737 title: str
1738 Plot title.
1739 vmin: float
1740 minimum for colorbar
1741 vmax: float
1742 maximum for colorbar
1743 hexbin: bool
1744 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1745 Else scatter of all points, with potential to overplot many comparisons on same plot.
1746 """
1747 if hexbin:
1748 # Check cubes using same functionality as the difference operator.
1749 if len(cubes) != 2:
1750 raise ValueError(
1751 "Cubes should contain exactly 2 cubes for hexbin plotting."
1752 )
1753 title = title.replace("scatter", "hexbin")
1754 filename = filename.replace("scatter", "hexbin")
1756 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1757 ax = plt.gca()
1759 model_colors_map = get_model_colors_map(cubes)
1761 percentiles = np.arange(0, 100, 5)
1762 percentiles[0] = 1
1763 percentiles[-1] = 99
1764 quantiles = iris.cube.CubeList()
1766 # Loop through all output cubes for both data points and overplotting quantiles.
1767 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1768 for plottype in ["points", "quantiles"]:
1769 nplot = 0
1770 for cube in iter_maybe(cubes):
1771 label = None
1772 color = "black"
1773 if model_colors_map: 1773 ↛ 1778line 1773 didn't jump to line 1778 because the condition on line 1773 was always true
1774 label = cube.attributes.get("model_name")
1775 color = model_colors_map[label]
1777 # Plot all data points
1778 if plottype == "points":
1779 if nplot > 0:
1780 if hexbin:
1781 hb = plt.hexbin(
1782 cubes[0].data.flatten(),
1783 cube.data.flatten(),
1784 alpha=0.3,
1785 gridsize=100,
1786 mincnt=1,
1787 )
1788 else:
1789 plt.scatter(
1790 cubes[0].data.flatten(),
1791 cube.data.flatten(),
1792 color=color,
1793 marker="+",
1794 label=None,
1795 alpha=0.3,
1796 )
1798 elif plottype == "quantiles": 1798 ↛ 1817line 1798 didn't jump to line 1817 because the condition on line 1798 was always true
1799 # Construct Q-Q plot
1800 quantiles.append(
1801 cube.collapsed(
1802 cube.coords(dim_coords=True),
1803 iris.analysis.PERCENTILE,
1804 percent=percentiles,
1805 )
1806 )
1807 if nplot > 0:
1808 iplt.scatter(
1809 quantiles[0],
1810 quantiles[-1],
1811 color=color,
1812 marker="o",
1813 label=label,
1814 edgecolors="black",
1815 )
1817 nplot = nplot + 1
1819 # Add some labels and tweak the style.
1820 ax.set_title(title, fontsize=16)
1821 ax.set_xlabel(
1822 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1823 )
1824 ax.set_ylabel(
1825 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1826 )
1827 ax.tick_params(axis="both", labelsize=12)
1828 ax.autoscale()
1830 # Set 1:1 line and equal axes if scatter plot of common cube names
1831 nameA = iter_maybe(cubes)[0].name()
1832 nameB = iter_maybe(cubes)[1].name()
1833 if any(part in nameB.split("_") for part in nameA.split("_")): 1833 ↛ 1844line 1833 didn't jump to line 1844 because the condition on line 1833 was always true
1834 lims = [
1835 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1836 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1837 ]
1838 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1839 ax.set_aspect("equal")
1840 ax.set_xlim(lims)
1841 ax.set_ylim(lims)
1843 # Overlay grid-lines onto scatter plot.
1844 ax.grid(linestyle="--", color="grey", linewidth=1)
1845 if model_colors_map: 1845 ↛ 1849line 1845 didn't jump to line 1849 because the condition on line 1845 was always true
1846 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1848 # Add colorbar if hexbin output
1849 if hexbin:
1850 cb = plt.colorbar(
1851 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1852 )
1853 cb.set_label("Number of data points", size=12)
1855 # Save plot.
1856 _save_close_figure(fig, "scatter", filename)
1859def _spatial_plot(
1860 method: Literal["contourf", "pcolormesh", "scatter"],
1861 cube: iris.cube.Cube,
1862 filename: str | None,
1863 sequence_coordinate: str,
1864 stamp_coordinate: str,
1865 overlay_cube: iris.cube.Cube | None = None,
1866 contour_cube: iris.cube.Cube | None = None,
1867 point_cube: iris.cube.Cube | None = None,
1868 **kwargs,
1869):
1870 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1872 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1873 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1874 is present then postage stamp plots will be produced.
1876 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1877 be overplotted on the same figure.
1879 Parameters
1880 ----------
1881 method: "contourf" | "pcolormesh" | "scatter"
1882 The plotting method to use.
1883 Select choice of "contourf" or "pcolormesh" for gridded data.
1884 Use "scatter" for point-based data.
1885 cube: Cube
1886 Iris cube of the data to plot. It should have two spatial dimensions,
1887 such as lat and lon, and may also have a another two dimension to be
1888 plotted sequentially and/or as postage stamp plots.
1889 filename: str | None
1890 Name of the plot to write, used as a prefix for plot sequences. If None
1891 uses the recipe name.
1892 sequence_coordinate: str
1893 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1894 This coordinate must exist in the cube.
1895 stamp_coordinate: str
1896 Coordinate about which to plot postage stamp plots. Defaults to
1897 ``"realization"``.
1898 overlay_cube: Cube | None, optional
1899 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1900 contour_cube: Cube | None, optional
1901 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1902 point_cube: Cube | None, optional
1903 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
1905 Raises
1906 ------
1907 ValueError
1908 If the cube doesn't have the right dimensions.
1909 TypeError
1910 If the cube isn't a single cube.
1911 """
1912 # Ensure we've got a single cube.
1913 cube = check_single_cube(cube)
1915 # Set title based on recipe metadata or use cube name
1916 recipe_title = get_recipe_metadata().get("title", cube.name())
1918 # Check if there is a valid stamp coordinate in cube dimensions.
1919 if stamp_coordinate == "realization": 1919 ↛ 1924line 1919 didn't jump to line 1924 because the condition on line 1919 was always true
1920 stamp_coordinate = check_stamp_coordinate(cube)
1922 # Make postage stamp plots if stamp_coordinate exists and has more than a
1923 # single point.
1924 plotting_func = _plot_and_save_spatial_plot
1925 try:
1926 if cube.coord(stamp_coordinate).shape[0] > 1:
1927 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1928 except iris.exceptions.CoordinateNotFoundError:
1929 pass
1931 # Produce a geographical scatter plot if the data have a
1932 # dimension called observation or model_obs_error
1933 if any(
1934 crd.var_name == "station"
1935 or crd.var_name == "Station_Name"
1936 or crd.var_name == "model_obs_error"
1937 for crd in cube.coords()
1938 ):
1939 plotting_func = _plot_and_save_spatial_plot
1940 method = "scatter"
1942 # Must have a sequence coordinate.
1943 try:
1944 cube.coord(sequence_coordinate)
1945 except iris.exceptions.CoordinateNotFoundError as err:
1946 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1948 # Create a plot for each value of the sequence coordinate.
1949 plot_index = []
1950 nplot = np.size(cube.coord(sequence_coordinate).points)
1952 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1953 # Set plot titles and filename
1954 seq_coord = cube_slice.coord(sequence_coordinate)
1955 plot_title, plot_filename = _set_title_and_filename(
1956 seq_coord, nplot, recipe_title, filename
1957 )
1959 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1960 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1961 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1962 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1964 # Do the actual plotting.
1965 plotting_func(
1966 cube_slice,
1967 filename=plot_filename,
1968 stamp_coordinate=stamp_coordinate,
1969 title=plot_title,
1970 method=method,
1971 overlay_cube=overlay_slice,
1972 contour_cube=contour_slice,
1973 point_cube=point_slice,
1974 **kwargs,
1975 )
1976 plot_index.append(plot_filename)
1978 # Add list of plots to plot metadata.
1979 complete_plot_index = _append_to_plot_index(plot_index)
1981 # Make a page to display the plots.
1982 _make_plot_html_page(complete_plot_index)
1985####################
1986# Public functions #
1987####################
1990def spatial_contour_plot(
1991 cube: iris.cube.Cube,
1992 filename: str | None = None,
1993 sequence_coordinate: str = "time",
1994 stamp_coordinate: str = "realization",
1995 **kwargs,
1996) -> iris.cube.Cube:
1997 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1999 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2000 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2001 is present then postage stamp plots will be produced.
2003 Parameters
2004 ----------
2005 cube: Cube
2006 Iris cube of the data to plot. It should have two spatial dimensions,
2007 such as lat and lon, and may also have a another two dimension to be
2008 plotted sequentially and/or as postage stamp plots.
2009 filename: str, optional
2010 Name of the plot to write, used as a prefix for plot sequences. Defaults
2011 to the recipe name.
2012 sequence_coordinate: str, optional
2013 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2014 This coordinate must exist in the cube.
2015 stamp_coordinate: str, optional
2016 Coordinate about which to plot postage stamp plots. Defaults to
2017 ``"realization"``.
2019 Returns
2020 -------
2021 Cube
2022 The original cube (so further operations can be applied).
2024 Raises
2025 ------
2026 ValueError
2027 If the cube doesn't have the right dimensions.
2028 TypeError
2029 If the cube isn't a single cube.
2030 """
2031 _spatial_plot(
2032 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2033 )
2034 return cube
2037def spatial_pcolormesh_plot(
2038 cube: iris.cube.Cube,
2039 filename: str | None = None,
2040 sequence_coordinate: str = "time",
2041 stamp_coordinate: str = "realization",
2042 **kwargs,
2043) -> iris.cube.Cube:
2044 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2046 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2047 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2048 is present then postage stamp plots will be produced.
2050 This function is significantly faster than ``spatial_contour_plot``,
2051 especially at high resolutions, and should be preferred unless contiguous
2052 contour areas are important.
2054 Parameters
2055 ----------
2056 cube: Cube
2057 Iris cube of the data to plot. It should have two spatial dimensions,
2058 such as lat and lon, and may also have a another two dimension to be
2059 plotted sequentially and/or as postage stamp plots.
2060 filename: str, optional
2061 Name of the plot to write, used as a prefix for plot sequences. Defaults
2062 to the recipe name.
2063 sequence_coordinate: str, optional
2064 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2065 This coordinate must exist in the cube.
2066 stamp_coordinate: str, optional
2067 Coordinate about which to plot postage stamp plots. Defaults to
2068 ``"realization"``.
2070 Returns
2071 -------
2072 Cube
2073 The original cube (so further operations can be applied).
2075 Raises
2076 ------
2077 ValueError
2078 If the cube doesn't have the right dimensions.
2079 TypeError
2080 If the cube isn't a single cube.
2081 """
2082 _spatial_plot(
2083 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2084 )
2085 return cube
2088def spatial_multi_pcolormesh_plot(
2089 cube: iris.cube.Cube,
2090 overlay_cube: iris.cube.Cube | None = None,
2091 contour_cube: iris.cube.Cube | None = None,
2092 point_cube: iris.cube.Cube | None = None,
2093 filename: str | None = None,
2094 sequence_coordinate: str = "time",
2095 stamp_coordinate: str = "realization",
2096 **kwargs,
2097) -> iris.cube.Cube:
2098 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2100 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2101 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2102 is present then postage stamp plots will be produced.
2104 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2106 If specified, contours of a contour_cube can be overplotted on top of those.
2108 If specified, a spatial scatter map of point_cube can be overplotted.
2110 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2112 This function is significantly faster than ``spatial_contour_plot``,
2113 especially at high resolutions, and should be preferred unless contiguous
2114 contour areas are important.
2116 Parameters
2117 ----------
2118 cube: Cube
2119 Iris cube of the data to plot. It should have two spatial dimensions,
2120 such as lat and lon, and may also have two additional dimensions to be
2121 plotted sequentially and/or as postage stamp plots.
2122 overlay_cube: Cube, optional
2123 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2124 such as lat and lon, and may also have two additional dimensions to be
2125 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.
2126 If not provided, output plot generated without overlay cube.
2127 contour_cube: Cube, optional
2128 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,
2129 such as lat and lon, and may also have two additional dimensions to be
2130 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2131 point_cube: Cube, optional
2132 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
2133 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2134 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2135 postage stamp plots. If not provided, output plot generated without point-based layer.
2136 filename: str, optional
2137 Name of the plot to write, used as a prefix for plot sequences. Defaults
2138 to the recipe name.
2139 sequence_coordinate: str, optional
2140 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2141 This coordinate must exist in the cube.
2142 stamp_coordinate: str, optional
2143 Coordinate about which to plot postage stamp plots. Defaults to
2144 ``"realization"``.
2146 Returns
2147 -------
2148 Cube
2149 The original cube (so further operations can be applied).
2151 Raises
2152 ------
2153 ValueError
2154 If the cube doesn't have the right dimensions.
2155 TypeError
2156 If the cube isn't a single cube.
2157 """
2158 _spatial_plot(
2159 "pcolormesh",
2160 cube,
2161 filename,
2162 sequence_coordinate,
2163 stamp_coordinate,
2164 overlay_cube=overlay_cube,
2165 contour_cube=contour_cube,
2166 point_cube=point_cube,
2167 )
2168 return cube, overlay_cube, contour_cube, point_cube
2171# TODO: Expand function to handle ensemble data.
2172# line_coordinate: str, optional
2173# Coordinate about which to plot multiple lines. Defaults to
2174# ``"realization"``.
2175def plot_line_series(
2176 cube: iris.cube.Cube | iris.cube.CubeList,
2177 filename: str | None = None,
2178 series_coordinate: str = "time",
2179 sequence_coordinate: str = "time",
2180 # add the following for ensembles
2181 stamp_coordinate: str = "realization",
2182 single_plot: bool = False,
2183 **kwargs,
2184) -> iris.cube.Cube | iris.cube.CubeList:
2185 """Plot a line plot for the specified coordinate.
2187 The Cube or CubeList must be 1D.
2189 Parameters
2190 ----------
2191 iris.cube | iris.cube.CubeList
2192 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2193 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2194 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2195 filename: str, optional
2196 Name of the plot to write, used as a prefix for plot sequences. Defaults
2197 to the recipe name.
2198 series_coordinate: str, optional
2199 Coordinate about which to make a series. Defaults to ``"time"``. This
2200 coordinate must exist in the cube.
2202 Returns
2203 -------
2204 iris.cube.Cube | iris.cube.CubeList
2205 The original Cube or CubeList (so further operations can be applied).
2207 Raises
2208 ------
2209 ValueError
2210 If the cubes don't have the right dimensions.
2211 TypeError
2212 If the cube isn't a Cube or CubeList.
2213 """
2214 # Ensure we have a name for the plot file.
2215 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2217 num_models = get_num_models(cube)
2219 validate_cube_shape(cube, num_models)
2221 # Iterate over all cubes and extract coordinate to plot.
2222 cubes = iris.cube.CubeList(iter_maybe(cube))
2223 coords = []
2224 for model_cube in cubes:
2225 try:
2226 coords.append(model_cube.coord(series_coordinate))
2227 except iris.exceptions.CoordinateNotFoundError as err:
2228 raise ValueError(
2229 f"Cube must have a {series_coordinate} coordinate."
2230 ) from err
2231 if model_cube.coords("realization") and model_cube.ndim > 2:
2232 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2234 plot_index = []
2236 # Check if this is a spectral plot by looking for spectral coordinates
2237 is_spectral_plot = series_coordinate in [
2238 "frequency",
2239 "physical_wavenumber",
2240 "wavelength",
2241 ]
2243 if is_spectral_plot:
2244 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2245 # coordinate frequency/wavenumber.
2246 # If several power spectra are plotted with time as sequence_coordinate for the
2247 # time slider option.
2249 # Internal plotting function.
2250 plotting_func = _plot_and_save_line_power_spectrum_series
2252 for model_cube in cubes:
2253 try:
2254 model_cube.coord(sequence_coordinate)
2255 except iris.exceptions.CoordinateNotFoundError as err:
2256 raise ValueError(
2257 f"Cube must have a {sequence_coordinate} coordinate."
2258 ) from err
2260 if num_models == 1: 2260 ↛ 2275line 2260 didn't jump to line 2275 because the condition on line 2260 was always true
2261 # check for ensembles
2262 if ( 2262 ↛ 2266line 2262 didn't jump to line 2266 because the condition on line 2262 was never true
2263 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2264 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2265 ):
2266 if single_plot:
2267 # Plot spectra, mean and ensemble spread on 1 plot
2268 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2269 else:
2270 # Plot postage stamps
2271 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2272 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2273 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2274 else:
2275 all_points = sorted(
2276 set(
2277 itertools.chain.from_iterable(
2278 cb.coord(sequence_coordinate).points for cb in cubes
2279 )
2280 )
2281 )
2282 all_slices = list(
2283 itertools.chain.from_iterable(
2284 cb.slices_over(sequence_coordinate) for cb in cubes
2285 )
2286 )
2287 # Matched slices (matched by seq coord point; it may happen that
2288 # evaluated models do not cover the same seq coord range, hence matching
2289 # necessary)
2290 cube_iterables = [
2291 iris.cube.CubeList(
2292 s
2293 for s in all_slices
2294 if s.coord(sequence_coordinate).points[0] == point
2295 )
2296 for point in all_points
2297 ]
2298 nplot = len(all_points)
2300 # Create a plot for each value of the sequence coordinate. Allowing for
2301 # multiple cubes in a CubeList to be plotted in the same plot for similar
2302 # sequence values. Passing a CubeList into the internal plotting function
2303 # for similar values of the sequence coordinate. cube_slice can be an
2304 # iris.cube.Cube or an iris.cube.CubeList.
2306 for cube_slice in cube_iterables:
2307 # Normalize cube_slice to a list of cubes
2308 if isinstance(cube_slice, iris.cube.CubeList): 2308 ↛ 2309line 2308 didn't jump to line 2309 because the condition on line 2308 was never true
2309 cubes = list(cube_slice)
2310 elif isinstance(cube_slice, iris.cube.Cube): 2310 ↛ 2313line 2310 didn't jump to line 2313 because the condition on line 2310 was always true
2311 cubes = [cube_slice]
2312 else:
2313 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2315 # Use sequence value so multiple sequences can merge.
2316 seq_coord = cube_slice[0].coord(sequence_coordinate)
2317 plot_title, plot_filename = _set_title_and_filename(
2318 seq_coord, nplot, recipe_title, filename
2319 )
2321 # Format the coordinate value in a unit appropriate way.
2322 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2324 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2325 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2325 ↛ 2326line 2325 didn't jump to line 2326 because the condition on line 2325 was never true
2326 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2328 # Do the actual plotting.
2329 plotting_func(
2330 cube_slice,
2331 coords,
2332 stamp_coordinate,
2333 plot_filename,
2334 title,
2335 series_coordinate,
2336 )
2338 plot_index.append(plot_filename)
2339 else:
2340 # Format the title and filename using plotted series coordinate
2341 nplot = 1
2342 seq_coord = coords[0]
2343 plot_title, plot_filename = _set_title_and_filename(
2344 seq_coord, nplot, recipe_title, filename
2345 )
2347 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2348 if (
2349 "station" in [c.name() for c in cubes[0].coords()]
2350 and len(cubes[0].coord("station").points) > 1
2351 ):
2352 for station in cubes[0].coord("station").points:
2353 station_cubes = cubes.extract(iris.Constraint(station=station))
2354 station_name = station_cubes[0].coord("Station_Name").points[0]
2355 station_plotname = plot_filename.replace(
2356 ".png", "_" + station_name + ".png"
2357 )
2358 _plot_and_save_line_series(
2359 station_cubes,
2360 coords,
2361 "realization",
2362 station_plotname,
2363 f"{plot_title} {station_name}",
2364 )
2365 plot_index.append(station_plotname)
2367 else:
2368 # Do the actual plotting for all other series coordinate options.
2369 _plot_and_save_line_series(
2370 cubes, coords, stamp_coordinate, plot_filename, plot_title
2371 )
2373 plot_index.append(plot_filename)
2375 # append plot to list of plots
2376 complete_plot_index = _append_to_plot_index(plot_index)
2378 # Make a page to display the plots.
2379 _make_plot_html_page(complete_plot_index)
2381 return cube
2384def plot_vertical_line_series(
2385 cubes: iris.cube.Cube | iris.cube.CubeList,
2386 filename: str | None = None,
2387 series_coordinate: str = "model_level_number",
2388 sequence_coordinate: str = "time",
2389 # line_coordinate: str = "realization",
2390 **kwargs,
2391) -> iris.cube.Cube | iris.cube.CubeList:
2392 """Plot a line plot against a type of vertical coordinate.
2394 The Cube or CubeList must be 1D.
2396 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2397 then a sequence of plots will be produced.
2399 Parameters
2400 ----------
2401 iris.cube | iris.cube.CubeList
2402 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2403 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2404 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2405 filename: str, optional
2406 Name of the plot to write, used as a prefix for plot sequences. Defaults
2407 to the recipe name.
2408 series_coordinate: str, optional
2409 Coordinate to plot on the y-axis. Can be ``pressure`` or
2410 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2411 for LFRic. Defaults to ``model_level_number``.
2412 This coordinate must exist in the cube.
2413 sequence_coordinate: str, optional
2414 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2415 This coordinate must exist in the cube.
2417 Returns
2418 -------
2419 iris.cube.Cube | iris.cube.CubeList
2420 The original Cube or CubeList (so further operations can be applied).
2421 Plotted data.
2423 Raises
2424 ------
2425 ValueError
2426 If the cubes doesn't have the right dimensions.
2427 TypeError
2428 If the cube isn't a Cube or CubeList.
2429 """
2430 # Ensure we have a name for the plot file.
2431 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2433 cubes = iter_maybe(cubes)
2434 # Initialise empty list to hold all data from all cubes in a CubeList
2435 all_data = []
2437 # Store min/max ranges for x range.
2438 x_levels = []
2440 num_models = get_num_models(cubes)
2442 validate_cube_shape(cubes, num_models)
2444 # Iterate over all cubes in cube or CubeList and plot.
2445 coords = []
2446 for cube in cubes:
2447 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2448 try:
2449 coords.append(cube.coord(series_coordinate))
2450 except iris.exceptions.CoordinateNotFoundError as err:
2451 raise ValueError(
2452 f"Cube must have a {series_coordinate} coordinate."
2453 ) from err
2455 try:
2456 if cube.ndim > 1 or not cube.coords("realization"): 2456 ↛ 2464line 2456 didn't jump to line 2464 because the condition on line 2456 was always true
2457 cube.coord(sequence_coordinate)
2458 except iris.exceptions.CoordinateNotFoundError as err:
2459 raise ValueError(
2460 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2461 ) from err
2463 # Get minimum and maximum from levels information.
2464 _, levels, _ = colorbar_map_levels(cube, axis="x")
2465 if levels is not None: 2465 ↛ 2469line 2465 didn't jump to line 2469 because the condition on line 2465 was always true
2466 x_levels.append(min(levels))
2467 x_levels.append(max(levels))
2468 else:
2469 all_data.append(cube.data)
2471 if len(x_levels) == 0: 2471 ↛ 2473line 2471 didn't jump to line 2473 because the condition on line 2471 was never true
2472 # Combine all data into a single NumPy array
2473 combined_data = np.concatenate(all_data)
2475 # Set the lower and upper limit for the x-axis to ensure all plots have
2476 # same range. This needs to read the whole cube over the range of the
2477 # sequence and if applicable postage stamp coordinate.
2478 vmin = np.floor(combined_data.min())
2479 vmax = np.ceil(combined_data.max())
2480 else:
2481 vmin = min(x_levels)
2482 vmax = max(x_levels)
2484 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2485 # a single profile directly without iterating over a sequence.
2486 sequence_coords = [
2487 cube.coord(sequence_coordinate)
2488 for cube in cubes
2489 if cube.coords(sequence_coordinate)
2490 ]
2491 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2492 np.size(coord.points) > 1 for coord in sequence_coords
2493 )
2494 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2495 np.size(coord.points) == 1 for coord in sequence_coords
2496 )
2498 plot_index = []
2499 if has_sequence_coord: 2499 ↛ 2524line 2499 didn't jump to line 2524 because the condition on line 2499 was always true
2500 # Matching the slices (matching by seq coord point; it may happen that
2501 # evaluated models do not cover the same seq coord range, hence matching
2502 # necessary)
2503 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2504 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2505 for cubes_slice in cube_iterables:
2506 # Format the coordinate value in a unit appropriate way.
2507 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2508 plot_title, plot_filename = _set_title_and_filename(
2509 seq_coord, nplot, recipe_title, filename
2510 )
2512 # Do the actual plotting.
2513 _plot_and_save_vertical_line_series(
2514 cubes_slice,
2515 coords,
2516 "realization",
2517 plot_filename,
2518 series_coordinate,
2519 title=plot_title,
2520 vmin=vmin,
2521 vmax=vmax,
2522 )
2523 plot_index.append(plot_filename)
2524 elif has_scalar_sequence_coord:
2525 # Scalar sequence coordinate (typically aggregated time bounds):
2526 # make one plot and include sequence period in title/filename.
2527 plot_title, plot_filename = _set_title_and_filename(
2528 sequence_coords[0], 1, recipe_title, filename
2529 )
2531 _plot_and_save_vertical_line_series(
2532 cubes,
2533 coords,
2534 "realization",
2535 plot_filename,
2536 series_coordinate,
2537 title=plot_title,
2538 vmin=vmin,
2539 vmax=vmax,
2540 )
2541 plot_index.append(plot_filename)
2542 else:
2543 # 1D case: no sequence coordinate, plot a single profile.
2544 plot_title = recipe_title
2545 if filename:
2546 plot_filename = filename
2547 else:
2548 plot_filename = f"{slugify(plot_title)}.png"
2550 _plot_and_save_vertical_line_series(
2551 cubes,
2552 coords,
2553 "realization",
2554 plot_filename,
2555 series_coordinate,
2556 title=plot_title,
2557 vmin=vmin,
2558 vmax=vmax,
2559 )
2560 plot_index.append(plot_filename)
2562 # Add list of plots to plot metadata.
2563 complete_plot_index = _append_to_plot_index(plot_index)
2565 # Make a page to display the plots.
2566 _make_plot_html_page(complete_plot_index)
2568 return cubes
2571def qq_plot(
2572 cubes: iris.cube.CubeList,
2573 coordinates: list[str],
2574 percentiles: list[float],
2575 model_names: list[str],
2576 filename: str | None = None,
2577 one_to_one: bool = True,
2578 **kwargs,
2579) -> iris.cube.CubeList:
2580 """Plot a Quantile-Quantile plot between two models for common time points.
2582 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2583 collapsed within the operator over all specified coordinates such as
2584 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2585 ensemble members to ensure a 1D cube (array).
2587 Parameters
2588 ----------
2589 cubes: iris.cube.CubeList
2590 Two cubes of the same variable with different models.
2591 coordinate: list[str]
2592 The list of coordinates to collapse over. This list should be
2593 every coordinate within the cube to result in a 1D cube around
2594 the percentile coordinate.
2595 percent: list[float]
2596 A list of percentiles to appear in the plot.
2597 model_names: list[str]
2598 A list of model names to appear on the axis of the plot.
2599 filename: str, optional
2600 Filename of the plot to write.
2601 one_to_one: bool, optional
2602 If True a 1:1 line is plotted; if False it is not. Default is True.
2604 Raises
2605 ------
2606 ValueError
2607 When the cubes are not compatible.
2609 Notes
2610 -----
2611 The quantile-quantile plot is a variant on the scatter plot representing
2612 two datasets by their quantiles (percentiles) for common time points.
2613 This plot does not use a theoretical distribution to compare against, but
2614 compares percentiles of two datasets. This plot does
2615 not use all raw data points, but plots the selected percentiles (quantiles) of
2616 each variable instead for the two datasets, thereby normalising the data for a
2617 direct comparison between the selected percentiles of the two dataset distributions.
2619 Quantile-quantile plots are valuable for comparing against
2620 observations and other models. Identical percentiles between the variables
2621 will lie on the one-to-one line implying the values correspond well to each
2622 other. Where there is a deviation from the one-to-one line a range of
2623 possibilities exist depending on how and where the data is shifted (e.g.,
2624 Wilks 2011 [Wilks2011]_).
2626 For distributions above the one-to-one line the distribution is left-skewed;
2627 below is right-skewed. A distinct break implies a bimodal distribution, and
2628 closer values/values further apart at the tails imply poor representation of
2629 the extremes.
2631 """
2632 # Check cubes using same functionality as the difference operator.
2633 if len(cubes) != 2:
2634 raise ValueError("cubes should contain exactly 2 cubes.")
2635 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2636 other: Cube = cubes.extract_cube(
2637 iris.Constraint(
2638 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2639 )
2640 )
2642 # Get spatial coord names.
2643 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2644 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2646 # Ensure cubes to compare are on common differencing grid.
2647 # This is triggered if either
2648 # i) latitude and longitude shapes are not the same. Note grid points
2649 # are not compared directly as these can differ through rounding
2650 # errors.
2651 # ii) or variables are known to often sit on different grid staggering
2652 # in different models (e.g. cell center vs cell edge), as is the case
2653 # for UM and LFRic comparisons.
2654 # In future greater choice of regridding method might be applied depending
2655 # on variable type. Linear regridding can in general be appropriate for smooth
2656 # variables. Care should be taken with interpretation of differences
2657 # given this dependency on regridding.
2658 if (
2659 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2660 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2661 ) or (
2662 base.long_name
2663 in [
2664 "eastward_wind_at_10m",
2665 "northward_wind_at_10m",
2666 "northward_wind_at_cell_centres",
2667 "eastward_wind_at_cell_centres",
2668 "zonal_wind_at_pressure_levels",
2669 "meridional_wind_at_pressure_levels",
2670 "potential_vorticity_at_pressure_levels",
2671 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2672 ]
2673 ):
2674 logger.debug("Linear regridding base cube to other grid to compute differences")
2675 base = regrid_onto_cube(base, other, method="Linear")
2677 # Extract just common time points.
2678 base, other = _extract_common_time_points(base, other)
2680 # Equalise attributes so we can merge.
2681 fully_equalise_attributes([base, other])
2682 logger.debug("Base: %s\nOther: %s", base, other)
2684 # Collapse cubes.
2685 base = collapse(
2686 base,
2687 coordinate=coordinates,
2688 method="PERCENTILE",
2689 additional_percent=percentiles,
2690 )
2691 other = collapse(
2692 other,
2693 coordinate=coordinates,
2694 method="PERCENTILE",
2695 additional_percent=percentiles,
2696 )
2698 # Ensure we have a name for the plot file.
2699 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2700 title = f"{recipe_title}"
2702 if filename is None:
2703 filename = slugify(recipe_title)
2705 # Add file extension.
2706 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2708 # Do the actual plotting on a scatter plot
2709 _plot_and_save_scatter_plot(
2710 base, other, plot_filename, title, one_to_one, model_names
2711 )
2713 # Add list of plots to plot metadata.
2714 plot_index = _append_to_plot_index([plot_filename])
2716 # Make a page to display the plots.
2717 _make_plot_html_page(plot_index)
2719 return iris.cube.CubeList([base, other])
2722def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2723 """
2724 Plot a Hinton style triangle/scorecard plot.
2726 This plot type can be useful for summarising high level information, such as comparing
2727 how 'skillful' two models are when verified against observations for a variety of metrics,
2728 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2729 than customisable by the user as input arguments; many have been designed to automatically
2730 scale the plot depending on the number of x and y components.
2732 Parameters
2733 ----------
2734 change: np.ndarray
2735 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2736 size/direction.
2737 signif: np.ndarray
2738 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2739 xaxis_labels: list
2740 List of labels for the xaxis (must match the second dimension length of signif and change,
2741 along with magnitude if not None).
2742 yaxis_labels: list
2743 List of labels for the yaxis (must match the first dimension length of signif and change,
2744 along with magnitude if not None).
2745 magnitude: np.ndarray | None
2746 Optional 2D array, matching the shape of change, signif, which contains numerical values
2747 the user wishes to display under each respective triangle.
2749 Returns
2750 -------
2751 matplotlib axes object to either display or do further modifications to.
2752 """
2753 # Setup colors of triangles
2754 color_pos = "#7CAE00"
2755 color_neg = "#7B68EE"
2757 # Setup cell/text size ratios
2758 figsize = None
2759 cell_size_in = 0.35
2760 text_row_ratio = 0.25
2762 # Ensure arrays, and change to bool for sig.
2763 change = np.asarray(change)
2764 signif = np.asarray(signif).astype(bool)
2765 if magnitude is not None: 2765 ↛ 2766line 2765 didn't jump to line 2766 because the condition on line 2765 was never true
2766 magnitude = np.asarray(magnitude)
2768 # Get the number of x and y elements
2769 ny, nx = change.shape
2771 # Build non-uniform y coordinates
2772 tri_height = 1.0
2773 txt_height = text_row_ratio
2775 tri_y = []
2776 txt_y = []
2777 y_edges = [0.0]
2779 y = 0.0
2780 for _j in range(ny):
2781 tri_y.append(y + tri_height / 2)
2782 y += tri_height
2783 y_edges.append(y)
2785 if magnitude is not None: 2785 ↛ 2786line 2785 didn't jump to line 2786 because the condition on line 2785 was never true
2786 txt_y.append(y + txt_height / 2)
2787 y += txt_height
2788 y_edges.append(y)
2790 total_height = y
2792 # Dynamic figure size
2793 if figsize is None: 2793 ↛ 2798line 2793 didn't jump to line 2798 because the condition on line 2793 was always true
2794 width = nx * cell_size_in
2795 height = total_height * cell_size_in + 2
2796 figsize = (width, height)
2798 fig, ax = plt.subplots(figsize=figsize)
2800 # Setup axes and grid.
2801 ax.set_aspect("equal", adjustable="box")
2802 ax.set_xlim(-0.5, nx - 0.5)
2803 ax.set_ylim(0, total_height)
2805 ax.set_xticks(np.arange(nx))
2806 ax.set_xticklabels(xaxis_labels, rotation=90)
2808 ax.set_yticks(tri_y)
2809 ax.set_yticklabels(yaxis_labels)
2811 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2812 ax.set_yticks(y_edges, minor=True)
2814 ax.set_axisbelow(True)
2815 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2816 ax.grid(False, which="major")
2817 ax.tick_params(which="minor", length=0)
2819 ax.invert_yaxis()
2821 # Compute marker scaling (fixed overlap)
2822 fig.canvas.draw()
2824 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2825 width_in, height_in = bbox.width, bbox.height
2827 cell_w = (width_in * fig.dpi) / nx
2828 cell_h = (height_in * fig.dpi) / total_height
2829 cell_pixels = min(cell_w, cell_h)
2831 max_marker_size = (0.6 * cell_pixels) ** 2
2833 text_fontsize = cell_pixels * 0.15
2835 # Plot triangles + text
2836 for j in range(ny):
2837 for i in range(nx):
2838 val = change[j, i]
2839 if np.isnan(val): 2839 ↛ 2840line 2839 didn't jump to line 2840 because the condition on line 2839 was never true
2840 continue
2842 if abs(val) < 0.01: 2842 ↛ 2843line 2842 didn't jump to line 2843 because the condition on line 2842 was never true
2843 continue
2845 sig = signif[j, i]
2846 size = max_marker_size * abs(val)
2848 # Triangle style
2849 if val >= 0:
2850 marker = "^"
2851 color = color_pos
2852 else:
2853 marker = "v"
2854 color = color_neg
2856 if sig:
2857 edgecolor = "black"
2858 linewidth = 0.6
2859 else:
2860 edgecolor = "none"
2861 linewidth = 0.0
2863 # Triangle
2864 ax.scatter(
2865 i,
2866 tri_y[j],
2867 s=size,
2868 marker=marker,
2869 c=color,
2870 edgecolors=edgecolor,
2871 linewidths=linewidth,
2872 zorder=3,
2873 clip_on=True, # ensures no rendering bleed
2874 )
2876 # Text row
2877 if magnitude is not None: 2877 ↛ 2878line 2877 didn't jump to line 2878 because the condition on line 2877 was never true
2878 mag_val = magnitude[j, i]
2880 if not np.isnan(mag_val):
2881 ax.text(
2882 i,
2883 txt_y[j],
2884 f"{mag_val:.1f}",
2885 ha="center",
2886 va="center",
2887 fontsize=text_fontsize,
2888 color="black",
2889 zorder=4,
2890 )
2892 plt.tight_layout()
2893 return fig, ax
2896def scatter_plot(
2897 cube_x: iris.cube.Cube | iris.cube.CubeList,
2898 cube_y: iris.cube.Cube | iris.cube.CubeList,
2899 filename: str | None = None,
2900 one_to_one: bool = True,
2901 **kwargs,
2902) -> iris.cube.CubeList:
2903 """Plot a scatter plot between two variables.
2905 Both cubes must be 1D.
2907 Parameters
2908 ----------
2909 cube_x: Cube | CubeList
2910 1 dimensional Cube of the data to plot on y-axis.
2911 cube_y: Cube | CubeList
2912 1 dimensional Cube of the data to plot on x-axis.
2913 filename: str, optional
2914 Filename of the plot to write.
2915 one_to_one: bool, optional
2916 If True a 1:1 line is plotted; if False it is not. Default is True.
2918 Returns
2919 -------
2920 cubes: CubeList
2921 CubeList of the original x and y cubes for further processing.
2923 Raises
2924 ------
2925 ValueError
2926 If the cube doesn't have the right dimensions and cubes not the same
2927 size.
2928 TypeError
2929 If the cube isn't a single cube.
2931 Notes
2932 -----
2933 Scatter plots are used for determining if there is a relationship between
2934 two variables. Positive relations have a slope going from bottom left to top
2935 right; Negative relations have a slope going from top left to bottom right.
2936 """
2937 # Iterate over all cubes in cube or CubeList and plot.
2938 for cube_iter in iter_maybe(cube_x):
2939 # Check cubes are correct shape.
2940 cube_iter = check_single_cube(cube_iter)
2941 if cube_iter.ndim > 1:
2942 raise ValueError("cube_x must be 1D.")
2944 # Iterate over all cubes in cube or CubeList and plot.
2945 for cube_iter in iter_maybe(cube_y):
2946 # Check cubes are correct shape.
2947 cube_iter = check_single_cube(cube_iter)
2948 if cube_iter.ndim > 1:
2949 raise ValueError("cube_y must be 1D.")
2951 # Ensure we have a name for the plot file.
2952 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2953 title = f"{recipe_title}"
2955 if filename is None:
2956 filename = slugify(recipe_title)
2958 # Add file extension.
2959 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2961 # Do the actual plotting.
2962 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2964 # Add list of plots to plot metadata.
2965 plot_index = _append_to_plot_index([plot_filename])
2967 # Make a page to display the plots.
2968 _make_plot_html_page(plot_index)
2970 return iris.cube.CubeList([cube_x, cube_y])
2973def vector_plot(
2974 cube_u: iris.cube.Cube,
2975 cube_v: iris.cube.Cube,
2976 filename: str | None = None,
2977 sequence_coordinate: str = "time",
2978 **kwargs,
2979) -> iris.cube.CubeList:
2980 """Plot a vector plot based on the input u and v components."""
2981 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2983 # Cubes must have a matching sequence coordinate.
2984 try:
2985 # Check that the u and v cubes have the same sequence coordinate.
2986 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2986 ↛ anywhereline 2986 didn't jump anywhere: it always raised an exception.
2987 raise ValueError("Coordinates do not match.")
2988 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2989 raise ValueError(
2990 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2991 ) from err
2993 # Create a plot for each value of the sequence coordinate.
2994 plot_index = []
2995 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2996 for cube_u_slice, cube_v_slice in zip(
2997 cube_u.slices_over(sequence_coordinate),
2998 cube_v.slices_over(sequence_coordinate),
2999 strict=True,
3000 ):
3001 # Format the coordinate value in a unit appropriate way.
3002 seq_coord = cube_u_slice.coord(sequence_coordinate)
3003 plot_title, plot_filename = _set_title_and_filename(
3004 seq_coord, nplot, recipe_title, filename
3005 )
3007 # Do the actual plotting.
3008 _plot_and_save_vector_plot(
3009 cube_u_slice,
3010 cube_v_slice,
3011 filename=plot_filename,
3012 title=plot_title,
3013 method="pcolormesh",
3014 )
3015 plot_index.append(plot_filename)
3017 # Add list of plots to plot metadata.
3018 complete_plot_index = _append_to_plot_index(plot_index)
3020 # Make a page to display the plots.
3021 _make_plot_html_page(complete_plot_index)
3023 return iris.cube.CubeList([cube_u, cube_v])
3026def plot_histogram_series(
3027 cubes: iris.cube.Cube | iris.cube.CubeList,
3028 filename: str | None = None,
3029 sequence_coordinate: str = "time",
3030 stamp_coordinate: str = "realization",
3031 single_plot: bool = False,
3032 **kwargs,
3033) -> iris.cube.Cube | iris.cube.CubeList:
3034 """Plot a histogram plot for each vertical level provided.
3036 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3037 is present then a sequence of plots will be produced using the time slider
3038 functionality to scroll through histograms against time. If a
3039 stamp_coordinate is present then postage stamp plots will be produced. If
3040 stamp_coordinate and single_plot is True, all postage stamp plots will be
3041 plotted in a single plot instead of separate postage stamp plots.
3043 Parameters
3044 ----------
3045 cubes: Cube | iris.cube.CubeList
3046 Iris cube or CubeList of the data to plot. It should have a single dimension other
3047 than the stamp coordinate.
3048 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3049 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3050 filename: str, optional
3051 Name of the plot to write, used as a prefix for plot sequences. Defaults
3052 to the recipe name.
3053 sequence_coordinate: str, optional
3054 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3055 This coordinate must exist in the cube and will be used for the time
3056 slider.
3057 stamp_coordinate: str, optional
3058 Coordinate about which to plot postage stamp plots. Defaults to
3059 ``"realization"``.
3060 single_plot: bool, optional
3061 If True, all postage stamp plots will be plotted in a single plot. If
3062 False, each postage stamp plot will be plotted separately. Is only valid
3063 if stamp_coordinate exists and has more than a single point.
3065 Returns
3066 -------
3067 iris.cube.Cube | iris.cube.CubeList
3068 The original Cube or CubeList (so further operations can be applied).
3069 Plotted data.
3071 Raises
3072 ------
3073 ValueError
3074 If the cube doesn't have the right dimensions.
3075 TypeError
3076 If the cube isn't a Cube or CubeList.
3077 """
3078 recipe_title = get_recipe_metadata().get("title", "Histogram")
3080 cubes = iter_maybe(cubes)
3082 # Internal plotting function.
3083 plotting_func = _plot_and_save_histogram_series
3085 num_models = get_num_models(cubes)
3087 validate_cube_shape(cubes, num_models)
3089 # If several histograms are plotted, check sequence_coordinate
3090 check_sequence_coordinate(cubes, sequence_coordinate)
3092 # Get axis minimum and maximum from levels information.
3093 # If no levels set, derive minima and maxima from data in CubeList.
3094 vmin, vmax = _set_axis_range(cubes)
3096 # Make postage stamp plots if stamp_coordinate exists and has more than a
3097 # single point. If single_plot is True:
3098 # -- all postage stamp plots will be plotted in a single plot instead of
3099 # separate postage stamp plots.
3100 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3101 # produced per single model only
3102 if num_models == 1:
3103 if ( 3103 ↛ 3107line 3103 didn't jump to line 3107 because the condition on line 3103 was never true
3104 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3105 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3106 ):
3107 if single_plot:
3108 plotting_func = (
3109 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3110 )
3111 else:
3112 plotting_func = _plot_and_save_postage_stamp_histogram_series
3113 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3114 else:
3115 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3117 plot_index = []
3118 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3119 # Create a plot for each value of the sequence coordinate. Allowing for
3120 # multiple cubes in a CubeList to be plotted in the same plot for similar
3121 # sequence values. Passing a CubeList into the internal plotting function
3122 # for similar values of the sequence coordinate. cube_slice can be an
3123 # iris.cube.Cube or an iris.cube.CubeList.
3124 for cube_slice in cube_iterables:
3125 single_cube = cube_slice
3126 if isinstance(cube_slice, iris.cube.CubeList):
3127 single_cube = cube_slice[0]
3129 # Ensure valid stamp coordinate in cube dimensions
3130 if stamp_coordinate == "realization": 3130 ↛ 3133line 3130 didn't jump to line 3133 because the condition on line 3130 was always true
3131 stamp_coordinate = check_stamp_coordinate(single_cube)
3132 # Set plot titles and filename, based on sequence coordinate
3133 seq_coord = single_cube.coord(sequence_coordinate)
3134 # Use time coordinate in title and filename if single histogram output.
3135 if sequence_coordinate == "realization" and nplot == 1: 3135 ↛ 3136line 3135 didn't jump to line 3136 because the condition on line 3135 was never true
3136 seq_coord = single_cube.coord("time")
3137 # Use station name in title and filename if model vs obs comparison
3138 if sequence_coordinate == "station": 3138 ↛ 3139line 3138 didn't jump to line 3139 because the condition on line 3138 was never true
3139 seq_coord = single_cube.coord("Station_Name")
3141 plot_title, plot_filename = _set_title_and_filename(
3142 seq_coord, nplot, recipe_title, filename
3143 )
3145 # Do the actual plotting.
3146 plotting_func(
3147 cube_slice,
3148 filename=plot_filename,
3149 stamp_coordinate=stamp_coordinate,
3150 title=plot_title,
3151 vmin=vmin,
3152 vmax=vmax,
3153 )
3154 plot_index.append(plot_filename)
3156 # Add list of plots to plot metadata.
3157 complete_plot_index = _append_to_plot_index(plot_index)
3159 # Make a page to display the plots.
3160 _make_plot_html_page(complete_plot_index)
3162 return cubes
3165def plot_scatter_series(
3166 cubes: iris.cube.Cube | iris.cube.CubeList,
3167 filename: str | None = None,
3168 sequence_coordinate: str = "time",
3169 stamp_coordinate: str = "realization",
3170 hexbin: bool = False,
3171 **kwargs,
3172) -> iris.cube.Cube | iris.cube.CubeList:
3173 """Plot a scatter plot for each sequence coordinate provided.
3175 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3176 is present then a sequence of plots will be produced using the time slider
3177 functionality to scroll through scatter against time. If a
3178 stamp_coordinate is present then postage stamp plots will be produced. If
3179 stamp_coordinate and single_plot is True, all postage stamp plots will be
3180 plotted in a single plot instead of separate postage stamp plots.
3182 Parameters
3183 ----------
3184 cubes: Cube | iris.cube.CubeList
3185 Iris cube or CubeList of the data to plot. It should have a single dimension other
3186 than the stamp coordinate.
3187 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3188 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3189 filename: str, optional
3190 Name of the plot to write, used as a prefix for plot sequences. Defaults
3191 to the recipe name.
3192 sequence_coordinate: str, optional
3193 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3194 This coordinate must exist in the cube and will be used for the time
3195 slider.
3196 stamp_coordinate: str, optional
3197 Coordinate about which to plot postage stamp plots. Defaults to
3198 ``"realization"``.
3199 hexbin: bool, optional
3200 If True, generate hexbin comparison plot.
3201 If False, generate point-by-point scatter plot.
3203 Returns
3204 -------
3205 iris.cube.Cube | iris.cube.CubeList
3206 The original Cube or CubeList (so further operations can be applied).
3207 Plotted data.
3209 Raises
3210 ------
3211 ValueError
3212 If the cube doesn't have the right dimensions.
3213 TypeError
3214 If the cube isn't a Cube or CubeList.
3215 """
3216 recipe_title = get_recipe_metadata().get("title", "Scatter")
3218 cubes = iter_maybe(cubes)
3220 # Internal plotting function.
3221 plotting_func = _plot_and_save_scatter_series
3223 num_models = get_num_models(cubes)
3225 validate_cube_shape(cubes, num_models)
3227 check_sequence_coordinate(cubes, sequence_coordinate)
3229 vmin, vmax = _set_axis_range(cubes)
3231 # Require >1 models to compare on scatter plot
3232 if num_models > 1:
3233 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3234 else:
3235 raise ValueError(
3236 "Scatter plot series requires multiple number of models in input data."
3237 )
3239 plot_index = []
3240 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3241 # Create a plot for each value of the sequence coordinate. Allowing for
3242 # multiple cubes in a CubeList to be plotted in the same plot for similar
3243 # sequence values. Passing a CubeList into the internal plotting function
3244 # for similar values of the sequence coordinate. cube_slice can be an
3245 # iris.cube.Cube or an iris.cube.CubeList.
3246 for cube_slice in cube_iterables:
3247 single_cube = cube_slice
3248 if isinstance(cube_slice, iris.cube.CubeList): 3248 ↛ 3252line 3248 didn't jump to line 3252 because the condition on line 3248 was always true
3249 single_cube = cube_slice[0]
3251 # Ensure valid stamp coordinate in cube dimensions
3252 if stamp_coordinate == "realization": 3252 ↛ 3255line 3252 didn't jump to line 3255 because the condition on line 3252 was always true
3253 stamp_coordinate = check_stamp_coordinate(single_cube)
3254 # Set plot titles and filename, based on sequence coordinate
3255 seq_coord = single_cube.coord(sequence_coordinate)
3256 # Use time coordinate in title and filename if single histogram output.
3257 if sequence_coordinate == "realization" and nplot == 1:
3258 seq_coord = single_cube.coord("time")
3259 # Use station name in title and filename if model vs obs comparison
3260 if sequence_coordinate == "station":
3261 seq_coord = single_cube.coord("Station_Name")
3263 plot_title, plot_filename = _set_title_and_filename(
3264 seq_coord, nplot, recipe_title, filename
3265 )
3267 # Do the actual plotting.
3268 plotting_func(
3269 cube_slice,
3270 filename=plot_filename,
3271 stamp_coordinate=stamp_coordinate,
3272 title=plot_title,
3273 vmin=vmin,
3274 vmax=vmax,
3275 hexbin=hexbin,
3276 )
3277 plot_index.append(plot_filename)
3279 # Add list of plots to plot metadata.
3280 complete_plot_index = _append_to_plot_index(plot_index)
3282 # Make a page to display the plots.
3283 _make_plot_html_page(complete_plot_index)
3285 return cubes
3288def _plot_and_save_postage_stamp_power_spectrum_series(
3289 cubes: iris.cube.Cube,
3290 coords: list[iris.coords.Coord],
3291 stamp_coordinate: str,
3292 filename: str,
3293 title: str,
3294 series_coordinate: str | None = None,
3295 **kwargs,
3296):
3297 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3299 Parameters
3300 ----------
3301 cubes: Cube or CubeList
3302 Cube or Cubelist of the power spectrum data.
3303 coords: list[Coord]
3304 Coordinates to plot on the x-axis, one per cube.
3305 stamp_coordinate: str
3306 Coordinate that becomes different plots.
3307 filename: str
3308 Filename of the plot to write.
3309 title: str
3310 Plot title.
3311 series_coordinate: str, optional
3312 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3314 """
3315 # Use the smallest square grid that will fit the members.
3316 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3318 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3319 model_colors_map = get_model_colors_map(cubes)
3320 # ax = plt.gca()
3321 # Make a subplot for each member.
3322 for member, subplot in zip(
3323 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3324 ):
3325 ax = plt.subplot(grid_size, grid_size, subplot)
3327 # Store min/max ranges.
3328 y_levels = []
3330 line_marker = None
3331 line_width = 1
3333 for cube in iter_maybe(member):
3334 xcoord = _select_series_coord(cube, series_coordinate)
3335 xname = xcoord.points
3337 yfield = cube.data # power spectrum
3338 label = None
3339 color = "black"
3340 if model_colors_map: 3340 ↛ 3341line 3340 didn't jump to line 3341 because the condition on line 3340 was never true
3341 label = cube.attributes.get("model_name")
3342 color = model_colors_map.get(label)
3344 if member.coord(stamp_coordinate).points == [0]:
3345 ax.plot(
3346 xname,
3347 yfield,
3348 color=color,
3349 marker=line_marker,
3350 ls="-",
3351 lw=line_width,
3352 label=f"{label} (control)"
3353 if len(cube.coord(stamp_coordinate).points) > 1
3354 else label,
3355 )
3356 # Label with member if part of an ensemble and not the control.
3357 else:
3358 ax.plot(
3359 xname,
3360 yfield,
3361 color=color,
3362 ls="-",
3363 lw=1.5,
3364 alpha=0.75,
3365 label=f"{label} (member)",
3366 )
3368 # Calculate the global min/max if multiple cubes are given.
3369 _, levels, _ = colorbar_map_levels(cube, axis="y")
3370 if levels is not None: 3370 ↛ 3371line 3370 didn't jump to line 3371 because the condition on line 3370 was never true
3371 y_levels.append(min(levels))
3372 y_levels.append(max(levels))
3374 # Add some labels and tweak the style.
3375 title = f"{title}"
3376 ax.set_title(title, fontsize=16)
3378 # Set appropriate x-axis label based on coordinate
3379 if series_coordinate == "wavelength" or ( 3379 ↛ 3382line 3379 didn't jump to line 3382 because the condition on line 3379 was never true
3380 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3381 ):
3382 ax.set_xlabel("Wavelength (km)", fontsize=14)
3383 elif series_coordinate == "physical_wavenumber" or ( 3383 ↛ 3388line 3383 didn't jump to line 3388 because the condition on line 3383 was always true
3384 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3385 ):
3386 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3387 else: # frequency or check units
3388 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3389 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3390 else:
3391 ax.set_xlabel("Wavenumber", fontsize=14)
3393 ax.set_ylabel("Power Spectral Density", fontsize=14)
3394 ax.tick_params(axis="both", labelsize=12)
3396 # Set log-log scale
3397 ax.set_xscale("log")
3398 ax.set_yscale("log")
3400 # Add gridlines
3401 ax.grid(linestyle="--", color="grey", linewidth=1)
3402 # Ientify unique labels for legend
3403 handles = list(
3404 {
3405 label: handle
3406 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3407 }.values()
3408 )
3409 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3411 ax = plt.gca()
3412 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3414 # Save plot.
3415 _save_close_figure(fig, "histogram postage stamp", filename)
3418def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3419 cubes: iris.cube.Cube,
3420 coords: list[iris.coords.Coord],
3421 stamp_coordinate: str,
3422 filename: str,
3423 title: str,
3424 series_coordinate: str | None = None,
3425 **kwargs,
3426):
3427 """Plot and save power spectra for ensemble members in single plot.
3429 Parameters
3430 ----------
3431 cubes: Cube or CubeList
3432 Cube or Cubelist of the power spectrum data.
3433 coords: list[Coord]
3434 Coordinates to plot on the x-axis, one per cube.
3435 stamp_coordinate: str
3436 Coordinate that becomes different plots.
3437 filename: str
3438 Filename of the plot to write.
3439 title: str
3440 Plot title.
3441 series_coordinate: str, optional
3442 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3444 """
3445 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3446 model_colors_map = get_model_colors_map(cubes)
3448 line_marker = None
3449 line_width = 1
3451 # Compute ensemble statistics to show spread
3452 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3453 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3454 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3456 xcoord_global = mean_cube.coord(series_coordinate)
3457 x_global = xcoord_global.points
3459 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3460 xcoord = _select_series_coord(member, series_coordinate)
3461 xname = xcoord.points
3463 yfield = member.data # power spectrum
3464 color = "black"
3465 if model_colors_map: 3465 ↛ 3469line 3465 didn't jump to line 3469 because the condition on line 3465 was always true
3466 label = member.attributes.get("model_name") if i == 0 else None
3467 color = model_colors_map.get(label)
3469 if member.coord(stamp_coordinate).points == [0]:
3470 ax.plot(
3471 xname,
3472 yfield,
3473 color=color,
3474 marker=line_marker,
3475 ls="-",
3476 lw=line_width,
3477 label=f"{label} (control)"
3478 if len(member.coord(stamp_coordinate).points) > 1
3479 else label,
3480 )
3481 # Label with member number if part of an ensemble and not the control.
3482 else:
3483 ax.plot(
3484 xname,
3485 yfield,
3486 color=color,
3487 ls="-",
3488 lw=1.5,
3489 alpha=0.75,
3490 label=label,
3491 )
3493 # Set appropriate x-axis label based on coordinate
3494 if series_coordinate == "wavelength" or ( 3494 ↛ 3497line 3494 didn't jump to line 3497 because the condition on line 3494 was never true
3495 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3496 ):
3497 ax.set_xlabel("Wavelength (km)", fontsize=14)
3498 elif series_coordinate == "physical_wavenumber" or ( 3498 ↛ 3503line 3498 didn't jump to line 3503 because the condition on line 3498 was always true
3499 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3500 ):
3501 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3502 else: # frequency or check units
3503 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3504 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3505 else:
3506 ax.set_xlabel("Wavenumber", fontsize=14)
3508 # Add ensemble spread shading
3509 ax.fill_between(
3510 x_global,
3511 min_cube.data,
3512 max_cube.data,
3513 color="grey",
3514 alpha=0.3,
3515 label="Ensemble spread",
3516 )
3518 # Add ensemble mean line
3519 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3521 ax.set_ylabel("Power Spectral Density", fontsize=14)
3522 ax.tick_params(axis="both", labelsize=12)
3524 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3525 # Set log-log scale
3526 ax.set_xscale("log")
3527 ax.set_yscale("log")
3529 # Add gridlines
3530 ax.grid(linestyle="--", color="grey", linewidth=1)
3531 # Identify unique labels for legend
3532 handles = list(
3533 {
3534 label: handle
3535 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3536 }.values()
3537 )
3538 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3540 # Figure title.
3541 ax.set_title(title, fontsize=16)
3543 # Save plot.
3544 _save_close_figure(fig, "power spectra postage stamp", filename)