Coverage for src/CSET/operators/plot.py: 83%
1107 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-12 16:03 +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 specific fixed ancillary spatial plots
256 if any(name in cube.name() for name in ("land_", "orography", "altitude")):
257 pass
258 else:
259 if cmap.name in ["viridis", "Greys"]:
260 coastcol = "magenta"
261 else:
262 coastcol = "black"
263 logger.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
264 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
265 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
267 # Add gridlines.
268 gl = axes.gridlines(
269 alpha=0.3,
270 draw_labels=True,
271 dms=False,
272 x_inline=False,
273 y_inline=False,
274 )
275 gl.top_labels = False
276 gl.right_labels = False
277 if subplot:
278 gl.bottom_labels = False
279 gl.left_labels = False
280 if subplot % grid_size[1] == 1:
281 gl.left_labels = True
282 if subplot > ((grid_size[0] - 1) * grid_size[1]): 282 ↛ 287line 282 didn't jump to line 287 because the condition on line 282 was always true
283 gl.bottom_labels = True
285 # If is lat/lon spatial map, fix extent to keep plot tight.
286 # Specifying crs within set_extent helps ensure only data region is shown.
287 if isinstance(
288 coord_system, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
289 ):
290 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
292 except ValueError:
293 # Skip if not both x and y map coordinates.
294 axes = figure.gca()
296 return axes
299def _get_plot_resolution() -> int:
300 """Get resolution of rasterised plots in pixels per inch."""
301 return get_recipe_metadata().get("plot_resolution", 100)
304def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
305 """Return title and filename based on start and end points or bounds."""
306 if use_bounds and seq_coord.has_bounds():
307 vals = seq_coord.bounds.flatten()
308 else:
309 vals = seq_coord.points
310 start = seq_coord.units.title(vals[0])
311 end = seq_coord.units.title(vals[-1])
313 if start == end:
314 sequence_title = f"\n [{start}]"
315 sequence_fname = f"_{filename_slugify(start)}"
316 else:
317 sequence_title = f"\n [{start} to {end}]"
318 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
320 # Do not include time if coord set to zero.
321 if (
322 seq_coord.units == "hours since 0001-01-01 00:00:00"
323 and vals[0] == 0
324 and vals[-1] == 0
325 ):
326 sequence_title = ""
327 sequence_fname = ""
329 return sequence_title, sequence_fname
332def _set_title_and_filename(
333 seq_coord: iris.coords.Coord,
334 nplot: int,
335 recipe_title: str,
336 filename: str,
337):
338 """Set plot title and filename based on cube coordinate.
340 Parameters
341 ----------
342 sequence_coordinate: iris.coords.Coord
343 Coordinate about which to make a plot sequence.
344 nplot: int
345 Number of output plots to generate - controls title/naming.
346 recipe_title: str
347 Default plot title, potentially to update.
348 filename: str
349 Input plot filename, potentially to update.
351 Returns
352 -------
353 plot_title: str
354 Output formatted plot title string, based on plotted data.
355 plot_filename: str
356 Output formatted plot filename string.
357 """
358 ndim = seq_coord.ndim
359 npoints = np.size(seq_coord.points)
360 sequence_title = ""
361 sequence_fname = ""
363 # Case 1: Multiple dimension sequence input - list number of aggregated cases
364 # (e.g. aggregation histogram plots)
365 if ndim > 1:
366 ncase = np.shape(seq_coord)[0]
367 sequence_title = f"\n [{ncase} cases]"
368 sequence_fname = f"_{ncase}cases"
370 # Case 2: Single dimension input
371 else:
372 # Single sequence point
373 if npoints == 1:
374 if nplot > 1:
375 # Default labels for sequence inputs
376 sequence_value = seq_coord.units.title(seq_coord.points[0])
377 sequence_value = sequence_value.replace(" unknown", "")
378 sequence_title = f"\n [{sequence_value}]"
379 sequence_fname = f"_{filename_slugify(sequence_value)}"
380 else:
381 # Aggregated attribute available where input collapsed over aggregation
382 try:
383 ncase = seq_coord.attributes["number_reference_times"]
384 sequence_title = f"\n [{ncase} cases]"
385 sequence_fname = f"_{ncase}cases"
386 except KeyError:
387 sequence_title, sequence_fname = _get_start_end_strings(
388 seq_coord, use_bounds=seq_coord.has_bounds()
389 )
390 # Multiple sequence (e.g. time) points
391 else:
392 sequence_title, sequence_fname = _get_start_end_strings(
393 seq_coord, use_bounds=False
394 )
396 # Set plot title and filename
397 plot_title = f"{recipe_title}{sequence_title}"
399 # Set plot filename, defaulting to user input if provided.
400 if filename is None:
401 filename = slugify(recipe_title)
402 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
403 else:
404 if nplot > 1:
405 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
406 else:
407 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
409 return plot_title, plot_filename
412def _select_series_coord(cube, series_coordinate):
413 """Determine the grid coordinates to use to calculate grid spacing."""
414 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
415 if series_coordinate in spacing_coordinates: 415 ↛ 421line 415 didn't jump to line 421 because the condition on line 415 was always true
416 # Try the requested coordinate first then the fallbacks in order.
417 fallbacks = [series_coordinate] + [
418 c for c in spacing_coordinates if c != series_coordinate
419 ]
420 else:
421 fallbacks = {series_coordinate}
423 # Try each possible coordinate.
424 for coord in fallbacks:
425 try:
426 return cube.coord(coord)
427 except iris.exceptions.CoordinateNotFoundError:
428 logger.debug("Coordinate %s not found.", coord)
430 # If we get here, none of the fallback options were found.
431 raise iris.exceptions.CoordinateNotFoundError(
432 f"No valid coordinate found for '{series_coordinate}' "
433 f"or fallback options {fallbacks}"
434 )
437def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
438 """Control postage stamp plot output titles based on stamp coordinate."""
439 if stamp_coord.name() == "realization":
440 mtitle = "Member"
441 else:
442 mtitle = stamp_coord.name().capitalize()
444 if stamp_coord.name() == "time":
445 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
446 else:
447 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
449 return mtitle
452def _set_axis_range(cubes):
453 """Get minimum and maximum from levels information."""
454 levels = None
455 for cube in cubes: 455 ↛ 471line 455 didn't jump to line 471 because the loop on line 455 didn't complete
456 # First check if user-specified "auto" range variable.
457 # This maintains the value of levels as None, so proceed.
458 _, levels, _ = colorbar_map_levels(cube, axis="y")
459 if levels is None:
460 break
461 # If levels is changed, recheck to use the vmin,vmax or
462 # levels-based ranges for histogram plots.
463 _, levels, _ = colorbar_map_levels(cube)
464 logger.debug("levels: %s", levels)
465 if levels is not None: 465 ↛ 455line 465 didn't jump to line 455 because the condition on line 465 was always true
466 vmin = min(levels)
467 vmax = max(levels)
468 logger.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
469 break
471 if levels is None:
472 vmin = min(cb.data.min() for cb in cubes)
473 vmax = max(cb.data.max() for cb in cubes)
475 return vmin, vmax
478def _find_matched_slices(cubes, sequence_coordinate):
479 """Identify matched cubes in CubeList by sequence_coordinate values.
481 Ensures common points are compared for multiple cube inputs.
482 """
483 all_points = sorted(
484 set(
485 itertools.chain.from_iterable(
486 cb.coord(sequence_coordinate).points for cb in cubes
487 )
488 )
489 )
490 all_slices = list(
491 itertools.chain.from_iterable(
492 cb.slices_over(sequence_coordinate) for cb in cubes
493 )
494 )
495 # Matched slices (matched by seq coord point; it may happen that
496 # evaluated models do not cover the same seq coord range, hence matching
497 # necessary)
498 cube_iterables = [
499 iris.cube.CubeList(
500 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
501 )
502 for point in all_points
503 ]
505 return cube_iterables
508def _plot_and_save_spatial_plot(
509 cube: iris.cube.Cube,
510 filename: str,
511 title: str,
512 method: Literal["contourf", "pcolormesh", "scatter"],
513 overlay_cube: iris.cube.Cube | None = None,
514 contour_cube: iris.cube.Cube | None = None,
515 point_cube: iris.cube.Cube | None = None,
516 **kwargs,
517):
518 """Plot and save a spatial plot.
520 Parameters
521 ----------
522 cube: Cube
523 2 dimensional (lat and lon) Cube of the data to plot.
524 filename: str
525 Filename of the plot to write.
526 title: str
527 Plot title.
528 method: "contourf" | "pcolormesh" | "scatter"
529 The plotting method to use
530 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
531 overlay_cube: Cube, optional
532 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
533 contour_cube: Cube, optional
534 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
535 point_cube: Cube, optional
536 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
537 """
538 # Setup plot details, size, resolution, etc.
539 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
541 # Specify the color bar
542 cmap, levels, norm = colorbar_map_levels(cube)
544 # If overplotting, set required colorbars
545 if overlay_cube:
546 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
547 if contour_cube:
548 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
550 # Setup plot map projection, extent and coastlines and borderlines.
551 axes = _setup_spatial_map(cube, fig, cmap)
553 # Set colorscale bounds
554 try:
555 vmin = min(levels)
556 vmax = max(levels)
557 except TypeError:
558 vmin, vmax = None, None
559 # Ensure to use norm and not vmin/vmax if levels are defined.
560 if norm is not None:
561 vmin = None
562 vmax = None
563 logger.debug("Plotting using defined levels.")
565 # Plot the field.
566 if method == "contourf":
567 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
568 elif method == "pcolormesh":
569 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
570 elif method == "scatter":
571 # Scatter plot of the field. The marker size is chosen to give
572 # symbols that decrease in size as the number of data points
573 # increases, although the fraction of the figure covered by
574 # symbols increases roughly as N^(1/2), disregarding overlaps,
575 # and has been selected for the default figure size of (10, 10).
576 # Should this be changed, the marker size should be adjusted in
577 # proportion to the area of the figure.
578 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
579 lat_axis, lon_axis = get_cube_yxcoordname(cube)
580 plot = iplt.scatter(
581 cube.coord(lon_axis),
582 cube.coord(lat_axis),
583 c=cube.data[:],
584 s=mrk_size,
585 cmap=cmap,
586 edgecolors="k",
587 norm=norm,
588 vmin=vmin,
589 vmax=vmax,
590 )
591 else:
592 raise ValueError(f"Unknown plotting method: {method}")
594 # Overplot overlay field, if required
595 if overlay_cube:
596 try:
597 over_vmin = min(over_levels)
598 over_vmax = max(over_levels)
599 except TypeError:
600 over_vmin, over_vmax = None, None
601 if over_norm is not None: 601 ↛ 602line 601 didn't jump to line 602 because the condition on line 601 was never true
602 over_vmin = None
603 over_vmax = None
604 overlay = iplt.pcolormesh(
605 overlay_cube,
606 cmap=over_cmap,
607 norm=over_norm,
608 alpha=0.8,
609 vmin=over_vmin,
610 vmax=over_vmax,
611 )
612 # Overplot contour field, if required, with contour labelling.
613 if contour_cube:
614 contour = iplt.contour(
615 contour_cube,
616 colors="darkgray",
617 levels=cntr_levels,
618 norm=cntr_norm,
619 alpha=0.5,
620 linestyles="--",
621 linewidths=1,
622 )
623 plt.clabel(contour)
624 # Overplot valid elements of point-based field, if required.
625 # Check for non-masked points only to avoid plotting missing data.
626 if point_cube:
627 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
628 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
629 lon_coord = point_cube.coord(lon_axis)
630 lat_coord = point_cube.coord(lat_axis)
631 valid = ~point_cube.data.mask
632 valid_lon = iris.coords.AuxCoord(
633 lon_coord.points[valid],
634 standard_name=lon_coord.standard_name,
635 units=lon_coord.units,
636 coord_system=lon_coord.coord_system,
637 )
638 valid_lat = iris.coords.AuxCoord(
639 lat_coord.points[valid],
640 standard_name=lat_coord.standard_name,
641 units=lat_coord.units,
642 coord_system=lat_coord.coord_system,
643 )
644 iplt.scatter(
645 valid_lon,
646 valid_lat,
647 c=point_cube.data[valid],
648 s=mrk_size,
649 cmap=cmap,
650 edgecolors="k",
651 norm=norm,
652 vmin=vmin,
653 vmax=vmax,
654 )
656 # Check to see if transect, and if so, adjust y axis.
657 if is_transect(cube):
658 if "pressure" in [coord.name() for coord in cube.coords()]:
659 axes.invert_yaxis()
660 axes.set_yscale("log")
661 axes.set_ylim(1100, 100)
662 # If both model_level_number and level_height exists, iplt can construct
663 # plot as a function of height above orography (NOT sea level).
664 elif {"model_level_number", "level_height"}.issubset( 664 ↛ 669line 664 didn't jump to line 669 because the condition on line 664 was always true
665 {coord.name() for coord in cube.coords()}
666 ):
667 axes.set_yscale("log")
669 axes.set_title(
670 f"{title}\n"
671 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
672 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
673 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
674 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
675 fontsize=16,
676 )
678 # Inset code
679 axins = inset_axes(
680 axes,
681 width="20%",
682 height="20%",
683 loc="upper right",
684 axes_class=GeoAxes,
685 axes_kwargs={"map_projection": ccrs.PlateCarree()},
686 )
688 # Slightly transparent to reduce plot blocking.
689 axins.patch.set_alpha(0.4)
691 axins.coastlines(resolution="50m")
692 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
694 SLat, SLon, ELat, ELon = (
695 float(coord) for coord in cube.attributes["transect_coords"].split("_")
696 )
698 # Draw line between them
699 axins.plot(
700 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
701 )
703 # Plot points (note: lon, lat order for Cartopy)
704 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
705 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
707 lon_min, lon_max = sorted([SLon, ELon])
708 lat_min, lat_max = sorted([SLat, ELat])
710 # Midpoints
711 lon_mid = (lon_min + lon_max) / 2
712 lat_mid = (lat_min + lat_max) / 2
714 # Maximum half-range
715 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
716 if half_range == 0: # points identical → provide small default 716 ↛ 720line 716 didn't jump to line 720 because the condition on line 716 was always true
717 half_range = 1
719 # Set square extent
720 axins.set_extent(
721 [
722 lon_mid - half_range,
723 lon_mid + half_range,
724 lat_mid - half_range,
725 lat_mid + half_range,
726 ],
727 crs=ccrs.PlateCarree(),
728 )
730 # Ensure square aspect
731 axins.set_aspect("equal")
733 else:
734 # Add title.
735 axes.set_title(title, fontsize=16)
737 # Adjust padding if spatial plot or transect
738 if is_transect(cube):
739 yinfopad = -0.1
740 ycbarpad = 0.1
741 else:
742 yinfopad = 0.01
743 ycbarpad = 0.042
745 # Add watermark with min/max/mean. Currently not user togglable.
746 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
747 axes.annotate(
748 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
749 xy=(0.025, yinfopad),
750 xycoords="axes fraction",
751 xytext=(-5, 5),
752 textcoords="offset points",
753 ha="left",
754 va="bottom",
755 size=11,
756 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
757 )
759 # Add secondary colour bar for overlay_cube field if required.
760 if overlay_cube:
761 cbarB = fig.colorbar(
762 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
763 )
764 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
765 # add ticks and tick_labels for every levels if less than 20 levels exist
766 if over_levels is not None and len(over_levels) < 20: 766 ↛ 767line 766 didn't jump to line 767 because the condition on line 766 was never true
767 cbarB.set_ticks(over_levels)
768 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
769 if any(
770 var in overlay_cube.name()
771 for var in ("rainfall", "snowfall", "visibility")
772 ):
773 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
774 logger.debug("Set secondary colorbar ticks and labels.")
776 # Add main colour bar.
777 cbar = fig.colorbar(
778 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
779 )
781 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
782 # add ticks and tick_labels for every levels if less than 20 levels exist
783 if levels is not None and len(levels) < 20:
784 cbar.set_ticks(levels)
785 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
786 if any(var in cube.name() for var in ("rainfall", "snowfall", "visibility")): 786 ↛ 789line 786 didn't jump to line 789 because the condition on line 786 was always true
787 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
788 # Tick labels for rainfall rates from Nimrod radar data.
789 if "rainfall rate composite" in cube.name(): 789 ↛ 790line 789 didn't jump to line 790 because the condition on line 789 was never true
790 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
791 # Tick labels for rain accumulations from Nimrod radar data.
792 if "rain accumulation" in cube.name(): 792 ↛ 793line 792 didn't jump to line 793 because the condition on line 792 was never true
793 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
794 if "wts accumulation" in cube.name(): 794 ↛ 795line 794 didn't jump to line 795 because the condition on line 794 was never true
795 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
796 cbar.minorticks_off()
797 cbar.set_ticks(tick_levels)
798 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
799 cbar.set_label(label=f"{cube.name()}", size=14)
800 # Tick labels for model rainfall data.
801 if "surface_microphysical" in cube.name(): 801 ↛ 804line 801 didn't jump to line 804 because the condition on line 801 was always true
802 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
803 # Tick labels for Nimrod weights data.
804 logger.debug("Set colorbar ticks and labels.")
806 # Save plot.
807 _save_close_figure(fig, "spatial", filename)
810def _plot_and_save_postage_stamp_spatial_plot(
811 cube: iris.cube.Cube,
812 filename: str,
813 stamp_coordinate: str,
814 title: str,
815 method: Literal["contourf", "pcolormesh"],
816 overlay_cube: iris.cube.Cube | None = None,
817 contour_cube: iris.cube.Cube | None = None,
818 **kwargs,
819):
820 """Plot postage stamp spatial plots from an ensemble.
822 Parameters
823 ----------
824 cube: Cube
825 Iris cube of data to be plotted. It must have the stamp coordinate.
826 filename: str
827 Filename of the plot to write.
828 stamp_coordinate: str
829 Coordinate that becomes different plots.
830 method: "contourf" | "pcolormesh"
831 The plotting method to use.
832 overlay_cube: Cube, optional
833 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
834 contour_cube: Cube, optional
835 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
837 Raises
838 ------
839 ValueError
840 If the cube doesn't have the right dimensions.
841 """
842 # Use the smallest square grid that will fit the members.
843 nmember = len(cube.coord(stamp_coordinate).points)
844 grid_rows = int(math.sqrt(nmember))
845 grid_size = math.ceil(nmember / grid_rows)
847 fig = plt.figure(
848 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
849 )
851 # Specify the color bar
852 cmap, levels, norm = colorbar_map_levels(cube)
853 # If overplotting, set required colorbars
854 if overlay_cube: 854 ↛ 855line 854 didn't jump to line 855 because the condition on line 854 was never true
855 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
856 if contour_cube: 856 ↛ 857line 856 didn't jump to line 857 because the condition on line 856 was never true
857 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
859 # Make a subplot for each member.
860 for member, subplot in zip(
861 cube.slices_over(stamp_coordinate),
862 range(1, grid_size * grid_rows + 1),
863 strict=False,
864 ):
865 # Setup subplot map projection, extent and coastlines and borderlines.
866 axes = _setup_spatial_map(
867 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
868 )
869 if method == "contourf":
870 # Filled contour plot of the field.
871 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
872 elif method == "pcolormesh":
873 if levels is not None:
874 vmin = min(levels)
875 vmax = max(levels)
876 else:
877 raise TypeError("Unknown vmin and vmax range.")
878 vmin, vmax = None, None
879 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
880 # if levels are defined.
881 if norm is not None: 881 ↛ 882line 881 didn't jump to line 882 because the condition on line 881 was never true
882 vmin = None
883 vmax = None
884 # pcolormesh plot of the field.
885 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
886 else:
887 raise ValueError(f"Unknown plotting method: {method}")
889 # Overplot overlay field, if required
890 if overlay_cube: 890 ↛ 891line 890 didn't jump to line 891 because the condition on line 890 was never true
891 try:
892 over_vmin = min(over_levels)
893 over_vmax = max(over_levels)
894 except TypeError:
895 over_vmin, over_vmax = None, None
896 if over_norm is not None:
897 over_vmin = None
898 over_vmax = None
899 iplt.pcolormesh(
900 overlay_cube[member.coord(stamp_coordinate).points[0]],
901 cmap=over_cmap,
902 norm=over_norm,
903 alpha=0.6,
904 vmin=over_vmin,
905 vmax=over_vmax,
906 )
907 # Overplot contour field, if required
908 if contour_cube: 908 ↛ 909line 908 didn't jump to line 909 because the condition on line 908 was never true
909 iplt.contour(
910 contour_cube[member.coord(stamp_coordinate).points[0]],
911 colors="darkgray",
912 levels=cntr_levels,
913 norm=cntr_norm,
914 alpha=0.6,
915 linestyles="--",
916 linewidths=1,
917 )
918 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
919 axes.set_title(f"{mtitle}")
921 # Put the shared colorbar in its own axes.
922 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
923 colorbar = fig.colorbar(
924 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
925 )
926 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
928 # Overall figure title.
929 fig.suptitle(title, fontsize=16)
931 # Save plot.
932 _save_close_figure(fig, "contour postate stamp", filename)
935def _plot_and_save_line_series(
936 cubes: iris.cube.CubeList,
937 coords: list[iris.coords.Coord],
938 ensemble_coord: str,
939 filename: str,
940 title: str,
941 **kwargs,
942):
943 """Plot and save a 1D line series.
945 Parameters
946 ----------
947 cubes: Cube or CubeList
948 Cube or CubeList containing the cubes to plot on the y-axis.
949 coords: list[Coord]
950 Coordinates to plot on the x-axis, one per cube.
951 ensemble_coord: str
952 Ensemble coordinate in the cube.
953 filename: str
954 Filename of the plot to write.
955 title: str
956 Plot title.
957 """
958 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
960 model_colors_map = get_model_colors_map(cubes)
962 # Store min/max ranges.
963 y_levels = []
965 # Check match-up across sequence coords gives consistent sizes
966 validate_cubes_coords(cubes, coords)
968 for cube, coord in zip(cubes, coords, strict=True):
969 label = None
970 color = "black"
971 if model_colors_map:
972 label = cube.attributes.get("model_name")
973 color = model_colors_map.get(label)
974 if not cube.coords(ensemble_coord): 974 ↛ 976line 974 didn't jump to line 976 because the condition on line 974 was never true
975 # No ensemble coordinate — plot the cube directly as a single line.
976 iplt.plot(coord, cube, color=color, marker="o", ls="-", lw=3, label=label)
977 else:
978 for cube_slice in cube.slices_over(ensemble_coord):
979 # Label with (control) if part of an ensemble or not otherwise.
980 if cube_slice.coord(ensemble_coord).points == [0]:
981 iplt.plot(
982 coord,
983 cube_slice,
984 color=color,
985 marker="o",
986 ls="-",
987 lw=3,
988 label=f"{label} (control)"
989 if len(cube.coord(ensemble_coord).points) > 1
990 else label,
991 )
992 # Label with (perturbed) if part of an ensemble and not the control.
993 else:
994 iplt.plot(
995 coord,
996 cube_slice,
997 color=color,
998 ls="-",
999 lw=1.5,
1000 alpha=0.75,
1001 label=f"{label} (member)",
1002 )
1004 # Calculate the global min/max if multiple cubes are given.
1005 _, levels, _ = colorbar_map_levels(cube, axis="y")
1006 if levels is not None: 1006 ↛ 1007line 1006 didn't jump to line 1007 because the condition on line 1006 was never true
1007 y_levels.append(min(levels))
1008 y_levels.append(max(levels))
1010 # Get the current axes.
1011 ax = plt.gca()
1013 # Add some labels and tweak the style.
1014 # check if cubes[0] works for single cube if not CubeList
1015 if coords[0].name() == "time":
1016 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1017 else:
1018 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1019 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1020 ax.set_title(title, fontsize=16)
1022 ax.ticklabel_format(axis="y", useOffset=False)
1023 ax.tick_params(axis="x", labelrotation=15)
1024 ax.tick_params(axis="both", labelsize=12)
1026 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1027 if y_levels: 1027 ↛ 1028line 1027 didn't jump to line 1028 because the condition on line 1027 was never true
1028 ax.set_ylim(min(y_levels), max(y_levels))
1029 logger.debug("Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels))
1030 else:
1031 ax.autoscale()
1033 # Add gridlines
1034 ax.grid(linestyle="--", color="grey", linewidth=1)
1035 # Add zero line
1036 ymin, ymax = ax.get_ylim()
1037 if ymin < 0.0 and ymax > 0.0:
1038 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1039 # Identify unique labels for legend
1040 handles = list(
1041 {
1042 label: handle
1043 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1044 }.values()
1045 )
1046 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1048 # Save plot.
1049 _save_close_figure(fig, "line", filename)
1052def _plot_and_save_line_power_spectrum_series(
1053 cubes: iris.cube.Cube | iris.cube.CubeList,
1054 coords: list[iris.coords.Coord],
1055 ensemble_coord: str,
1056 filename: str,
1057 title: str,
1058 series_coordinate: str,
1059 **kwargs,
1060):
1061 """Plot and save a 1D line series.
1063 Parameters
1064 ----------
1065 cubes: Cube or CubeList
1066 Cube or CubeList containing the cubes to plot on the y-axis.
1067 coords: list[Coord]
1068 Coordinates to plot on the x-axis, one per cube.
1069 ensemble_coord: str
1070 Ensemble coordinate in the cube.
1071 filename: str
1072 Filename of the plot to write.
1073 title: str
1074 Plot title.
1075 series_coordinate: str
1076 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1077 """
1078 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1079 model_colors_map = get_model_colors_map(cubes)
1080 ax = plt.gca()
1082 # Store min/max ranges.
1083 y_levels = []
1085 line_marker = None
1086 line_width = 1
1088 for cube in iter_maybe(cubes):
1089 # next 2 lines replace chunk of code.
1090 xcoord = _select_series_coord(cube, series_coordinate)
1091 xname = xcoord.points
1093 yfield = cube.data # power spectrum
1095 # If data from power spectra is all np.nans (like T+0h rainfall field which
1096 # might be full of zeros), then set yfield to zeros so it doesn't crash the
1097 # plotting.
1098 if np.all(np.isnan(yfield)):
1099 yfield = np.zeros_like(yfield)
1101 label = None
1102 color = "black"
1103 if model_colors_map: 1103 ↛ 1106line 1103 didn't jump to line 1106 because the condition on line 1103 was always true
1104 label = cube.attributes.get("model_name")
1105 color = model_colors_map.get(label)
1106 for cube_slice in cube.slices_over(ensemble_coord):
1107 # Label with (control) if part of an ensemble or not otherwise.
1108 if cube_slice.coord(ensemble_coord).points == [0]: 1108 ↛ 1122line 1108 didn't jump to line 1122 because the condition on line 1108 was always true
1109 ax.plot(
1110 xname,
1111 yfield,
1112 color=color,
1113 marker=line_marker,
1114 ls="-",
1115 lw=line_width,
1116 label=f"{label} (control)"
1117 if len(cube.coord(ensemble_coord).points) > 1
1118 else label,
1119 )
1120 # Label with (perturbed) if part of an ensemble and not the control.
1121 else:
1122 ax.plot(
1123 xname,
1124 yfield,
1125 color=color,
1126 ls="-",
1127 lw=1.5,
1128 alpha=0.75,
1129 label=f"{label} (member)",
1130 )
1132 # Calculate the global min/max if multiple cubes are given.
1133 _, levels, _ = colorbar_map_levels(cube, axis="y")
1134 if levels is not None: 1134 ↛ 1135line 1134 didn't jump to line 1135 because the condition on line 1134 was never true
1135 y_levels.append(min(levels))
1136 y_levels.append(max(levels))
1138 # Add some labels and tweak the style.
1140 title = f"{title}"
1141 ax.set_title(title, fontsize=16)
1143 # Set appropriate x-axis label based on coordinate
1144 if series_coordinate == "wavelength" or ( 1144 ↛ 1147line 1144 didn't jump to line 1147 because the condition on line 1144 was never true
1145 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1146 ):
1147 ax.set_xlabel("Wavelength (km)", fontsize=14)
1148 elif series_coordinate == "physical_wavenumber" or ( 1148 ↛ 1151line 1148 didn't jump to line 1151 because the condition on line 1148 was never true
1149 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1150 ):
1151 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1152 else: # frequency or check units
1153 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1153 ↛ 1154line 1153 didn't jump to line 1154 because the condition on line 1153 was never true
1154 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1155 else:
1156 ax.set_xlabel("Wavenumber", fontsize=14)
1158 ax.set_ylabel("Power Spectral Density", fontsize=14)
1159 ax.tick_params(axis="both", labelsize=12)
1161 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1163 # Set log-log scale
1164 ax.set_xscale("log")
1165 ax.set_yscale("log")
1167 # Add gridlines
1168 ax.grid(linestyle="--", color="grey", linewidth=1)
1169 # Ientify unique labels for legend
1170 handles = list(
1171 {
1172 label: handle
1173 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1174 }.values()
1175 )
1176 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1178 # Save plot.
1179 _save_close_figure(fig, "line power spectrum", filename)
1182def _plot_and_save_vertical_line_series(
1183 cubes: iris.cube.CubeList,
1184 coords: list[iris.coords.Coord],
1185 ensemble_coord: str,
1186 filename: str,
1187 series_coordinate: str,
1188 title: str,
1189 vmin: float,
1190 vmax: float,
1191 **kwargs,
1192):
1193 """Plot and save a 1D line series in vertical.
1195 Parameters
1196 ----------
1197 cubes: CubeList
1198 1 dimensional Cube or CubeList of the data to plot on x-axis.
1199 coord: list[Coord]
1200 Coordinates to plot on the y-axis, one per cube.
1201 ensemble_coord: str
1202 Ensemble coordinate in the cube.
1203 filename: str
1204 Filename of the plot to write.
1205 series_coordinate: str
1206 Coordinate to use as vertical axis.
1207 title: str
1208 Plot title.
1209 vmin: float
1210 Minimum value for the x-axis.
1211 vmax: float
1212 Maximum value for the x-axis.
1213 """
1214 # plot the vertical pressure axis using log scale
1215 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1217 model_colors_map = get_model_colors_map(cubes)
1219 # Check match-up across sequence coords gives consistent sizes
1220 validate_cubes_coords(cubes, coords)
1222 for cube, coord in zip(cubes, coords, strict=True):
1223 label = None
1224 color = "black"
1225 if model_colors_map: 1225 ↛ 1226line 1225 didn't jump to line 1226 because the condition on line 1225 was never true
1226 label = cube.attributes.get("model_name")
1227 color = model_colors_map.get(label)
1229 for cube_slice in cube.slices_over(ensemble_coord):
1230 # If ensemble data given plot control member with (control)
1231 # unless single forecast.
1232 if cube_slice.coord(ensemble_coord).points == [0]:
1233 iplt.plot(
1234 cube_slice,
1235 coord,
1236 color=color,
1237 marker="o",
1238 ls="-",
1239 lw=3,
1240 label=f"{label} (control)"
1241 if len(cube.coord(ensemble_coord).points) > 1
1242 else label,
1243 )
1244 # If ensemble data given plot perturbed members with (perturbed).
1245 else:
1246 iplt.plot(
1247 cube_slice,
1248 coord,
1249 color=color,
1250 ls="-",
1251 lw=1.5,
1252 alpha=0.75,
1253 label=f"{label} (member)",
1254 )
1256 # Get the current axis
1257 ax = plt.gca()
1259 # Special handling for pressure level data.
1260 if series_coordinate == "pressure": 1260 ↛ 1282line 1260 didn't jump to line 1282 because the condition on line 1260 was always true
1261 # Invert y-axis and set to log scale.
1262 ax.invert_yaxis()
1263 ax.set_yscale("log")
1265 # Define y-ticks and labels for pressure log axis.
1266 y_tick_labels = [
1267 "1000",
1268 "850",
1269 "700",
1270 "500",
1271 "300",
1272 "200",
1273 "100",
1274 ]
1275 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1277 # Set y-axis limits and ticks.
1278 ax.set_ylim(1100, 100)
1280 # Test if series_coordinate is model level data. The UM data uses
1281 # model_level_number and lfric uses full_levels as coordinate.
1282 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1283 # Define y-ticks and labels for vertical axis.
1284 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1285 y_tick_labels = [str(int(i)) for i in y_ticks]
1286 ax.set_ylim(min(y_ticks), max(y_ticks))
1288 ax.set_yticks(y_ticks)
1289 ax.set_yticklabels(y_tick_labels)
1291 # Set x-axis limits.
1292 ax.set_xlim(vmin, vmax)
1293 # Mark y=0 if present in plot.
1294 if vmin < 0.0 and vmax > 0.0: 1294 ↛ 1295line 1294 didn't jump to line 1295 because the condition on line 1294 was never true
1295 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1297 # Add some labels and tweak the style.
1298 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1299 ax.set_xlabel(
1300 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1301 )
1302 ax.set_title(title, fontsize=16)
1303 ax.ticklabel_format(axis="x")
1304 ax.tick_params(axis="y")
1305 ax.tick_params(axis="both", labelsize=12)
1307 # Add gridlines
1308 ax.grid(linestyle="--", color="grey", linewidth=1)
1309 # Ientify unique labels for legend
1310 handles = list(
1311 {
1312 label: handle
1313 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1314 }.values()
1315 )
1316 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1318 # Save plot.
1319 _save_close_figure(fig, "vertical line", filename)
1322def _plot_and_save_scatter_plot(
1323 cube_x: iris.cube.Cube | iris.cube.CubeList,
1324 cube_y: iris.cube.Cube | iris.cube.CubeList,
1325 filename: str,
1326 title: str,
1327 one_to_one: bool,
1328 model_names: list[str] | None = None,
1329 **kwargs,
1330):
1331 """Plot and save a 2D scatter plot.
1333 Parameters
1334 ----------
1335 cube_x: Cube | CubeList
1336 1 dimensional Cube or CubeList of the data to plot on x-axis.
1337 cube_y: Cube | CubeList
1338 1 dimensional Cube or CubeList of the data to plot on y-axis.
1339 filename: str
1340 Filename of the plot to write.
1341 title: str
1342 Plot title.
1343 one_to_one: bool
1344 Whether a 1:1 line is plotted.
1345 """
1346 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1347 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1348 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1349 # over the pairs simultaneously.
1351 # Ensure cube_x and cube_y are iterable
1352 cube_x_iterable = iter_maybe(cube_x)
1353 cube_y_iterable = iter_maybe(cube_y)
1355 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1356 iplt.scatter(cube_x_iter, cube_y_iter)
1357 if one_to_one is True:
1358 plt.plot(
1359 [
1360 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1361 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1362 ],
1363 [
1364 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1365 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1366 ],
1367 "k",
1368 linestyle="--",
1369 )
1370 ax = plt.gca()
1372 # Add some labels and tweak the style.
1373 if model_names is None:
1374 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1375 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1376 else:
1377 # Add the model names, these should be order of base (x) and other (y).
1378 ax.set_xlabel(
1379 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1380 )
1381 ax.set_ylabel(
1382 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1383 )
1384 ax.set_title(title, fontsize=16)
1385 ax.ticklabel_format(axis="y", useOffset=False)
1386 ax.tick_params(axis="x", labelrotation=15)
1387 ax.tick_params(axis="both", labelsize=12)
1388 ax.autoscale()
1390 # Save plot.
1391 _save_close_figure(fig, "scatter", filename)
1394def _plot_and_save_vector_plot(
1395 cube_u: iris.cube.Cube,
1396 cube_v: iris.cube.Cube,
1397 filename: str,
1398 title: str,
1399 method: Literal["contourf", "pcolormesh"],
1400 **kwargs,
1401):
1402 """Plot and save a 2D vector plot.
1404 Parameters
1405 ----------
1406 cube_u: Cube
1407 2 dimensional Cube of u component of the data.
1408 cube_v: Cube
1409 2 dimensional Cube of v component of the data.
1410 filename: str
1411 Filename of the plot to write.
1412 title: str
1413 Plot title.
1414 """
1415 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1416 # Create a cube containing the magnitude of the vector field.
1417 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1418 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1419 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1420 cube_vec_mag.rename(
1421 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1422 )
1424 # Specify the color bar
1425 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1427 # Setup plot map projection, extent and coastlines and borderlines.
1428 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1430 if method == "contourf":
1431 # Filled contour plot of the field.
1432 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1433 elif method == "pcolormesh":
1434 try:
1435 vmin = min(levels)
1436 vmax = max(levels)
1437 except TypeError:
1438 vmin, vmax = None, None
1439 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1440 # if levels are defined.
1441 if norm is not None:
1442 vmin = None
1443 vmax = None
1444 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1445 else:
1446 raise ValueError(f"Unknown plotting method: {method}")
1448 # Check to see if transect, and if so, adjust y axis.
1449 if is_transect(cube_vec_mag):
1450 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1451 axes.invert_yaxis()
1452 axes.set_yscale("log")
1453 axes.set_ylim(1100, 100)
1454 # If both model_level_number and level_height exists, iplt can construct
1455 # plot as a function of height above orography (NOT sea level).
1456 elif {"model_level_number", "level_height"}.issubset(
1457 {coord.name() for coord in cube_vec_mag.coords()}
1458 ):
1459 axes.set_yscale("log")
1461 axes.set_title(
1462 f"{title}\n"
1463 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1464 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1465 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1466 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1467 fontsize=16,
1468 )
1470 else:
1471 # Add title.
1472 axes.set_title(title, fontsize=16)
1474 # Add watermark with min/max/mean. Currently not user togglable.
1475 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1476 axes.annotate(
1477 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}",
1478 xy=(0.05, -0.05),
1479 xycoords="axes fraction",
1480 xytext=(-5, 5),
1481 textcoords="offset points",
1482 ha="right",
1483 va="bottom",
1484 size=11,
1485 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1486 )
1488 # Add colour bar.
1489 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1490 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1491 # add ticks and tick_labels for every levels if less than 20 levels exist
1492 if levels is not None and len(levels) < 20:
1493 cbar.set_ticks(levels)
1494 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1496 # 30 barbs along the longest axis of the plot, or a barb per point for data
1497 # with less than 30 points.
1498 step = max(max(cube_u.shape) // 30, 1)
1499 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1501 # Save plot.
1502 _save_close_figure(fig, "vector", filename)
1505def _plot_and_save_histogram_series(
1506 cubes: iris.cube.Cube | iris.cube.CubeList,
1507 filename: str,
1508 title: str,
1509 vmin: float,
1510 vmax: float,
1511 **kwargs,
1512):
1513 """Plot and save a histogram series.
1515 Parameters
1516 ----------
1517 cubes: Cube or CubeList
1518 2 dimensional Cube or CubeList of the data to plot as histogram.
1519 filename: str
1520 Filename of the plot to write.
1521 title: str
1522 Plot title.
1523 vmin: float
1524 minimum for colorbar
1525 vmax: float
1526 maximum for colorbar
1527 """
1528 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1529 ax = plt.gca()
1531 model_colors_map = get_model_colors_map(cubes)
1533 # Set default that histograms will produce probability density function
1534 # at each bin (integral over range sums to 1).
1535 density = True
1537 for cube in iter_maybe(cubes):
1538 # Easier to check title (where var name originates)
1539 # than seeing if long names exist etc.
1540 # Exception case, where distribution better fits log scales/bins.
1541 if (
1542 ("surface_microphysical" in title)
1543 or ("rain accumulation" in title)
1544 or ("Rainfall rate Composite" in title)
1545 or ("Nimrod_5min" in title)
1546 ):
1547 if "amount" in title:
1548 # Compute histogram following Klingaman et al. (2017): ASoP
1549 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1550 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1551 density = False
1552 else:
1553 bins = 10.0 ** (
1554 np.arange(-10, 27, 1) / 10.0
1555 ) # Suggestion from RMED toolbox.
1556 bins = np.insert(bins, 0, 0)
1557 ax.set_yscale("log")
1558 vmin = bins[1]
1559 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1560 ax.set_xscale("log")
1561 elif "lightning" in title:
1562 bins = [0, 1, 2, 3, 4, 5]
1563 else:
1564 bins = np.linspace(vmin, vmax, 51)
1565 logger.debug(
1566 "Plotting histogram with %s bins %s - %s.",
1567 np.size(bins),
1568 np.min(bins),
1569 np.max(bins),
1570 )
1572 # Reshape cube data into a single array to allow for a single histogram.
1573 # Otherwise we plot xdim histograms stacked.
1574 cube_data_1d = (cube.data).flatten()
1576 label = None
1577 color = "black"
1578 if model_colors_map:
1579 label = cube.attributes.get("model_name")
1580 color = model_colors_map[label]
1581 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1583 # Compute area under curve.
1584 if (
1585 ("surface_microphysical" in title and "amount" in title)
1586 or ("rain_accumulation" in title)
1587 or ("Rainfall rate Composite" in title)
1588 or ("Nimrod_5min" in title)
1589 ):
1590 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1591 x = x * bin_mean / x.sum()
1592 x = x[1:]
1593 y = y[1:]
1595 ax.plot(
1596 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1597 )
1599 # Add some labels and tweak the style.
1600 ax.set_title(title, fontsize=16)
1601 ax.set_xlabel(
1602 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1603 )
1604 ax.set_ylabel("Normalised probability density", fontsize=14)
1605 if (
1606 ("surface_microphysical" in title and "amount" in title)
1607 or ("rain accumulation" in title)
1608 or ("Nimrod_5min" in title)
1609 ):
1610 ax.set_ylabel(
1611 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1612 )
1613 ax.set_xlim(vmin, vmax)
1614 ax.tick_params(axis="both", labelsize=12)
1616 # Overlay grid-lines onto histogram plot.
1617 ax.grid(linestyle="--", color="grey", linewidth=1)
1618 if model_colors_map:
1619 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1621 # Save plot.
1622 _save_close_figure(fig, "histogram", filename)
1625def _plot_and_save_postage_stamp_histogram_series(
1626 cube: iris.cube.Cube,
1627 filename: str,
1628 title: str,
1629 stamp_coordinate: str,
1630 vmin: float,
1631 vmax: float,
1632 **kwargs,
1633):
1634 """Plot and save postage (ensemble members) stamps for a histogram series.
1636 Parameters
1637 ----------
1638 cube: Cube
1639 2 dimensional Cube of the data to plot as histogram.
1640 filename: str
1641 Filename of the plot to write.
1642 title: str
1643 Plot title.
1644 stamp_coordinate: str
1645 Coordinate that becomes different plots.
1646 vmin: float
1647 minimum for pdf x-axis
1648 vmax: float
1649 maximum for pdf x-axis
1650 """
1651 # Use the smallest square grid that will fit the members.
1652 nmember = len(cube.coord(stamp_coordinate).points)
1653 grid_rows = int(math.sqrt(nmember))
1654 grid_size = math.ceil(nmember / grid_rows)
1656 fig = plt.figure(
1657 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1658 )
1659 # Make a subplot for each member.
1660 for member, subplot in zip(
1661 cube.slices_over(stamp_coordinate),
1662 range(1, grid_size * grid_rows + 1),
1663 strict=False,
1664 ):
1665 # Implicit interface is much easier here, due to needing to have the
1666 # cartopy GeoAxes generated.
1667 plt.subplot(grid_rows, grid_size, subplot)
1668 # Reshape cube data into a single array to allow for a single histogram.
1669 # Otherwise we plot xdim histograms stacked.
1670 member_data_1d = (member.data).flatten()
1671 plt.hist(member_data_1d, density=True, stacked=True)
1672 axes = plt.gca()
1673 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1674 axes.set_title(f"{mtitle}")
1675 axes.set_xlim(vmin, vmax)
1677 # Overall figure title.
1678 fig.suptitle(title, fontsize=16)
1680 # Save plot.
1681 _save_close_figure(fig, "histogram postage stamp", filename)
1684def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1685 cube: iris.cube.Cube,
1686 filename: str,
1687 title: str,
1688 stamp_coordinate: str,
1689 vmin: float,
1690 vmax: float,
1691 **kwargs,
1692):
1693 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1694 ax.set_title(title, fontsize=16)
1695 ax.set_xlim(vmin, vmax)
1696 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1697 ax.set_ylabel("normalised probability density", fontsize=14)
1698 # Loop over all slices along the stamp_coordinate
1699 for member in cube.slices_over(stamp_coordinate):
1700 # Flatten the member data to 1D
1701 member_data_1d = member.data.flatten()
1702 # Plot the histogram using plt.hist
1703 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1704 plt.hist(
1705 member_data_1d,
1706 density=True,
1707 stacked=True,
1708 label=f"{mtitle}",
1709 )
1711 # Add a legend
1712 ax.legend(fontsize=16)
1714 # Save plot.
1715 _save_close_figure(fig, "histogram postage stamp", filename)
1718def _plot_and_save_scatter_series(
1719 cubes: iris.cube.Cube | iris.cube.CubeList,
1720 filename: str,
1721 title: str,
1722 vmin: float,
1723 vmax: float,
1724 hexbin: bool,
1725 **kwargs,
1726):
1727 """Plot and save a scatter plot series.
1729 Parameters
1730 ----------
1731 cubes: Cube or CubeList
1732 2 dimensional Cube or CubeList of the data to plot as scatter.
1733 filename: str
1734 Filename of the plot to write.
1735 title: str
1736 Plot title.
1737 vmin: float
1738 minimum for colorbar
1739 vmax: float
1740 maximum for colorbar
1741 hexbin: bool
1742 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1743 Else scatter of all points, with potential to overplot many comparisons on same plot.
1744 """
1745 if hexbin:
1746 # Check cubes using same functionality as the difference operator.
1747 if len(cubes) != 2:
1748 raise ValueError(
1749 "Cubes should contain exactly 2 cubes for hexbin plotting."
1750 )
1751 title = title.replace("scatter", "hexbin")
1752 filename = filename.replace("scatter", "hexbin")
1754 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1755 ax = plt.gca()
1757 model_colors_map = get_model_colors_map(cubes)
1759 percentiles = np.arange(0, 100, 5)
1760 percentiles[0] = 1
1761 percentiles[-1] = 99
1762 quantiles = iris.cube.CubeList()
1764 # Loop through all output cubes for both data points and overplotting quantiles.
1765 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1766 for plottype in ["points", "quantiles"]:
1767 nplot = 0
1768 for cube in iter_maybe(cubes):
1769 label = None
1770 color = "black"
1771 if model_colors_map: 1771 ↛ 1776line 1771 didn't jump to line 1776 because the condition on line 1771 was always true
1772 label = cube.attributes.get("model_name")
1773 color = model_colors_map[label]
1775 # Plot all data points
1776 if plottype == "points":
1777 if nplot > 0:
1778 if hexbin:
1779 hb = plt.hexbin(
1780 cubes[0].data.flatten(),
1781 cube.data.flatten(),
1782 alpha=0.3,
1783 gridsize=100,
1784 mincnt=1,
1785 )
1786 else:
1787 plt.scatter(
1788 cubes[0].data.flatten(),
1789 cube.data.flatten(),
1790 color=color,
1791 marker="+",
1792 label=None,
1793 alpha=0.3,
1794 )
1796 elif plottype == "quantiles": 1796 ↛ 1815line 1796 didn't jump to line 1815 because the condition on line 1796 was always true
1797 # Construct Q-Q plot
1798 quantiles.append(
1799 cube.collapsed(
1800 cube.coords(dim_coords=True),
1801 iris.analysis.PERCENTILE,
1802 percent=percentiles,
1803 )
1804 )
1805 if nplot > 0:
1806 iplt.scatter(
1807 quantiles[0],
1808 quantiles[-1],
1809 color=color,
1810 marker="o",
1811 label=label,
1812 edgecolors="black",
1813 )
1815 nplot = nplot + 1
1817 # Add some labels and tweak the style.
1818 ax.set_title(title, fontsize=16)
1819 ax.set_xlabel(
1820 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1821 )
1822 ax.set_ylabel(
1823 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1824 )
1825 ax.tick_params(axis="both", labelsize=12)
1826 ax.autoscale()
1828 # Set 1:1 line and equal axes if scatter plot of common cube names
1829 nameA = iter_maybe(cubes)[0].name()
1830 nameB = iter_maybe(cubes)[1].name()
1831 if any(part in nameB.split("_") for part in nameA.split("_")): 1831 ↛ 1842line 1831 didn't jump to line 1842 because the condition on line 1831 was always true
1832 lims = [
1833 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1834 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1835 ]
1836 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1837 ax.set_aspect("equal")
1838 ax.set_xlim(lims)
1839 ax.set_ylim(lims)
1841 # Overlay grid-lines onto scatter plot.
1842 ax.grid(linestyle="--", color="grey", linewidth=1)
1843 if model_colors_map: 1843 ↛ 1847line 1843 didn't jump to line 1847 because the condition on line 1843 was always true
1844 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1846 # Add colorbar if hexbin output
1847 if hexbin:
1848 cb = plt.colorbar(
1849 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1850 )
1851 cb.set_label("Number of data points", size=12)
1853 # Save plot.
1854 _save_close_figure(fig, "scatter", filename)
1857def _spatial_plot(
1858 method: Literal["contourf", "pcolormesh", "scatter"],
1859 cube: iris.cube.Cube,
1860 filename: str | None,
1861 sequence_coordinate: str,
1862 stamp_coordinate: str,
1863 overlay_cube: iris.cube.Cube | None = None,
1864 contour_cube: iris.cube.Cube | None = None,
1865 point_cube: iris.cube.Cube | None = None,
1866 **kwargs,
1867):
1868 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1870 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1871 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1872 is present then postage stamp plots will be produced.
1874 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1875 be overplotted on the same figure.
1877 Parameters
1878 ----------
1879 method: "contourf" | "pcolormesh" | "scatter"
1880 The plotting method to use.
1881 Select choice of "contourf" or "pcolormesh" for gridded data.
1882 Use "scatter" for point-based data.
1883 cube: Cube
1884 Iris cube of the data to plot. It should have two spatial dimensions,
1885 such as lat and lon, and may also have a another two dimension to be
1886 plotted sequentially and/or as postage stamp plots.
1887 filename: str | None
1888 Name of the plot to write, used as a prefix for plot sequences. If None
1889 uses the recipe name.
1890 sequence_coordinate: str
1891 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1892 This coordinate must exist in the cube.
1893 stamp_coordinate: str
1894 Coordinate about which to plot postage stamp plots. Defaults to
1895 ``"realization"``.
1896 overlay_cube: Cube | None, optional
1897 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1898 contour_cube: Cube | None, optional
1899 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1900 point_cube: Cube | None, optional
1901 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
1903 Raises
1904 ------
1905 ValueError
1906 If the cube doesn't have the right dimensions.
1907 TypeError
1908 If the cube isn't a single cube.
1909 """
1910 # Ensure we've got a single cube.
1911 cube = check_single_cube(cube)
1913 # Set title based on recipe metadata or use cube name
1914 recipe_title = get_recipe_metadata().get("title", cube.name())
1916 # Check if there is a valid stamp coordinate in cube dimensions.
1917 if stamp_coordinate == "realization": 1917 ↛ 1922line 1917 didn't jump to line 1922 because the condition on line 1917 was always true
1918 stamp_coordinate = check_stamp_coordinate(cube)
1920 # Make postage stamp plots if stamp_coordinate exists and has more than a
1921 # single point.
1922 plotting_func = _plot_and_save_spatial_plot
1923 try:
1924 if cube.coord(stamp_coordinate).shape[0] > 1:
1925 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1926 except iris.exceptions.CoordinateNotFoundError:
1927 pass
1929 # Produce a geographical scatter plot if the data have a
1930 # dimension called observation or model_obs_error
1931 if any(
1932 crd.var_name == "station"
1933 or crd.var_name == "Station_Name"
1934 or crd.var_name == "model_obs_error"
1935 for crd in cube.coords()
1936 ):
1937 plotting_func = _plot_and_save_spatial_plot
1938 method = "scatter"
1940 # Must have a sequence coordinate.
1941 try:
1942 cube.coord(sequence_coordinate)
1943 except iris.exceptions.CoordinateNotFoundError as err:
1944 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1946 # Create a plot for each value of the sequence coordinate.
1947 plot_index = []
1948 nplot = np.size(cube.coord(sequence_coordinate).points)
1950 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1951 # Set plot titles and filename
1952 seq_coord = cube_slice.coord(sequence_coordinate)
1953 plot_title, plot_filename = _set_title_and_filename(
1954 seq_coord, nplot, recipe_title, filename
1955 )
1957 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1958 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1959 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1960 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1962 # Do the actual plotting.
1963 plotting_func(
1964 cube_slice,
1965 filename=plot_filename,
1966 stamp_coordinate=stamp_coordinate,
1967 title=plot_title,
1968 method=method,
1969 overlay_cube=overlay_slice,
1970 contour_cube=contour_slice,
1971 point_cube=point_slice,
1972 **kwargs,
1973 )
1974 plot_index.append(plot_filename)
1976 # Add list of plots to plot metadata.
1977 complete_plot_index = _append_to_plot_index(plot_index)
1979 # Make a page to display the plots.
1980 _make_plot_html_page(complete_plot_index)
1983####################
1984# Public functions #
1985####################
1988def spatial_contour_plot(
1989 cube: iris.cube.Cube,
1990 filename: str | None = None,
1991 sequence_coordinate: str = "time",
1992 stamp_coordinate: str = "realization",
1993 **kwargs,
1994) -> iris.cube.Cube:
1995 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1997 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1998 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1999 is present then postage stamp plots will be produced.
2001 Parameters
2002 ----------
2003 cube: Cube
2004 Iris cube of the data to plot. It should have two spatial dimensions,
2005 such as lat and lon, and may also have a another two dimension to be
2006 plotted sequentially and/or as postage stamp plots.
2007 filename: str, optional
2008 Name of the plot to write, used as a prefix for plot sequences. Defaults
2009 to the recipe name.
2010 sequence_coordinate: str, optional
2011 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2012 This coordinate must exist in the cube.
2013 stamp_coordinate: str, optional
2014 Coordinate about which to plot postage stamp plots. Defaults to
2015 ``"realization"``.
2017 Returns
2018 -------
2019 Cube
2020 The original cube (so further operations can be applied).
2022 Raises
2023 ------
2024 ValueError
2025 If the cube doesn't have the right dimensions.
2026 TypeError
2027 If the cube isn't a single cube.
2028 """
2029 _spatial_plot(
2030 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2031 )
2032 return cube
2035def spatial_pcolormesh_plot(
2036 cube: iris.cube.Cube,
2037 filename: str | None = None,
2038 sequence_coordinate: str = "time",
2039 stamp_coordinate: str = "realization",
2040 **kwargs,
2041) -> iris.cube.Cube:
2042 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2044 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2045 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2046 is present then postage stamp plots will be produced.
2048 This function is significantly faster than ``spatial_contour_plot``,
2049 especially at high resolutions, and should be preferred unless contiguous
2050 contour areas are important.
2052 Parameters
2053 ----------
2054 cube: Cube
2055 Iris cube of the data to plot. It should have two spatial dimensions,
2056 such as lat and lon, and may also have a another two dimension to be
2057 plotted sequentially and/or as postage stamp plots.
2058 filename: str, optional
2059 Name of the plot to write, used as a prefix for plot sequences. Defaults
2060 to the recipe name.
2061 sequence_coordinate: str, optional
2062 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2063 This coordinate must exist in the cube.
2064 stamp_coordinate: str, optional
2065 Coordinate about which to plot postage stamp plots. Defaults to
2066 ``"realization"``.
2068 Returns
2069 -------
2070 Cube
2071 The original cube (so further operations can be applied).
2073 Raises
2074 ------
2075 ValueError
2076 If the cube doesn't have the right dimensions.
2077 TypeError
2078 If the cube isn't a single cube.
2079 """
2080 _spatial_plot(
2081 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2082 )
2083 return cube
2086def spatial_multi_pcolormesh_plot(
2087 cube: iris.cube.Cube,
2088 overlay_cube: iris.cube.Cube | None = None,
2089 contour_cube: iris.cube.Cube | None = None,
2090 point_cube: iris.cube.Cube | None = None,
2091 filename: str | None = None,
2092 sequence_coordinate: str = "time",
2093 stamp_coordinate: str = "realization",
2094 **kwargs,
2095) -> iris.cube.Cube:
2096 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2098 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2099 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2100 is present then postage stamp plots will be produced.
2102 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2104 If specified, contours of a contour_cube can be overplotted on top of those.
2106 If specified, a spatial scatter map of point_cube can be overplotted.
2108 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2110 This function is significantly faster than ``spatial_contour_plot``,
2111 especially at high resolutions, and should be preferred unless contiguous
2112 contour areas are important.
2114 Parameters
2115 ----------
2116 cube: Cube
2117 Iris cube of the data to plot. It should have two spatial dimensions,
2118 such as lat and lon, and may also have two additional dimensions to be
2119 plotted sequentially and/or as postage stamp plots.
2120 overlay_cube: Cube, optional
2121 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2122 such as lat and lon, and may also have two additional dimensions to be
2123 plotted sequentially and/or as postage stamp plots. This is likely to be a masked cube in order not to hide the underlying basis cube.
2124 If not provided, output plot generated without overlay cube.
2125 contour_cube: Cube, optional
2126 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,
2127 such as lat and lon, and may also have two additional dimensions to be
2128 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2129 point_cube: Cube, optional
2130 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
2131 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2132 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2133 postage stamp plots. If not provided, output plot generated without point-based layer.
2134 filename: str, optional
2135 Name of the plot to write, used as a prefix for plot sequences. Defaults
2136 to the recipe name.
2137 sequence_coordinate: str, optional
2138 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2139 This coordinate must exist in the cube.
2140 stamp_coordinate: str, optional
2141 Coordinate about which to plot postage stamp plots. Defaults to
2142 ``"realization"``.
2144 Returns
2145 -------
2146 Cube
2147 The original cube (so further operations can be applied).
2149 Raises
2150 ------
2151 ValueError
2152 If the cube doesn't have the right dimensions.
2153 TypeError
2154 If the cube isn't a single cube.
2155 """
2156 _spatial_plot(
2157 "pcolormesh",
2158 cube,
2159 filename,
2160 sequence_coordinate,
2161 stamp_coordinate,
2162 overlay_cube=overlay_cube,
2163 contour_cube=contour_cube,
2164 point_cube=point_cube,
2165 )
2166 return cube, overlay_cube, contour_cube, point_cube
2169# TODO: Expand function to handle ensemble data.
2170# line_coordinate: str, optional
2171# Coordinate about which to plot multiple lines. Defaults to
2172# ``"realization"``.
2173def plot_line_series(
2174 cube: iris.cube.Cube | iris.cube.CubeList,
2175 filename: str | None = None,
2176 series_coordinate: str = "time",
2177 sequence_coordinate: str = "time",
2178 # add the following for ensembles
2179 stamp_coordinate: str = "realization",
2180 single_plot: bool = False,
2181 **kwargs,
2182) -> iris.cube.Cube | iris.cube.CubeList:
2183 """Plot a line plot for the specified coordinate.
2185 The Cube or CubeList must be 1D.
2187 Parameters
2188 ----------
2189 iris.cube | iris.cube.CubeList
2190 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2191 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2192 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2193 filename: str, optional
2194 Name of the plot to write, used as a prefix for plot sequences. Defaults
2195 to the recipe name.
2196 series_coordinate: str, optional
2197 Coordinate about which to make a series. Defaults to ``"time"``. This
2198 coordinate must exist in the cube.
2200 Returns
2201 -------
2202 iris.cube.Cube | iris.cube.CubeList
2203 The original Cube or CubeList (so further operations can be applied).
2205 Raises
2206 ------
2207 ValueError
2208 If the cubes don't have the right dimensions.
2209 TypeError
2210 If the cube isn't a Cube or CubeList.
2211 """
2212 # Ensure we have a name for the plot file.
2213 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2215 num_models = get_num_models(cube)
2217 validate_cube_shape(cube, num_models)
2219 # Iterate over all cubes and extract coordinate to plot.
2220 cubes = iris.cube.CubeList(iter_maybe(cube))
2221 coords = []
2222 for model_cube in cubes:
2223 try:
2224 coords.append(model_cube.coord(series_coordinate))
2225 except iris.exceptions.CoordinateNotFoundError as err:
2226 raise ValueError(
2227 f"Cube must have a {series_coordinate} coordinate."
2228 ) from err
2229 # Count dimensions excluding realization
2230 ndim = model_cube.ndim
2232 if model_cube.coords("realization"): 2232 ↛ 2239line 2232 didn't jump to line 2239 because the condition on line 2232 was always true
2233 realization_dims = model_cube.coord_dims("realization")
2235 # Only subtract if realization is a dimension coordinate
2236 if realization_dims:
2237 ndim -= len(realization_dims)
2239 if ndim > 2:
2240 raise ValueError(
2241 "Cube must be 1D or 2D (excluding any realization dimension)."
2242 )
2244 plot_index = []
2246 # Check if this is a spectral plot by looking for spectral coordinates
2247 is_spectral_plot = series_coordinate in [
2248 "frequency",
2249 "physical_wavenumber",
2250 "wavelength",
2251 ]
2253 if is_spectral_plot:
2254 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2255 # coordinate frequency/wavenumber.
2256 # If several power spectra are plotted with time as sequence_coordinate for the
2257 # time slider option.
2259 # Internal plotting function.
2260 plotting_func = _plot_and_save_line_power_spectrum_series
2262 for model_cube in cubes:
2263 try:
2264 model_cube.coord(sequence_coordinate)
2265 except iris.exceptions.CoordinateNotFoundError as err:
2266 raise ValueError(
2267 f"Cube must have a {sequence_coordinate} coordinate."
2268 ) from err
2270 if num_models == 1: 2270 ↛ 2285line 2270 didn't jump to line 2285 because the condition on line 2270 was always true
2271 # check for ensembles
2272 if ( 2272 ↛ 2276line 2272 didn't jump to line 2276 because the condition on line 2272 was never true
2273 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2274 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2275 ):
2276 if single_plot:
2277 # Plot spectra, mean and ensemble spread on 1 plot
2278 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2279 else:
2280 # Plot postage stamps
2281 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2282 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2283 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2284 else:
2285 all_points = sorted(
2286 set(
2287 itertools.chain.from_iterable(
2288 cb.coord(sequence_coordinate).points for cb in cubes
2289 )
2290 )
2291 )
2292 all_slices = list(
2293 itertools.chain.from_iterable(
2294 cb.slices_over(sequence_coordinate) for cb in cubes
2295 )
2296 )
2297 # Matched slices (matched by seq coord point; it may happen that
2298 # evaluated models do not cover the same seq coord range, hence matching
2299 # necessary)
2300 cube_iterables = [
2301 iris.cube.CubeList(
2302 s
2303 for s in all_slices
2304 if s.coord(sequence_coordinate).points[0] == point
2305 )
2306 for point in all_points
2307 ]
2308 nplot = len(all_points)
2310 # Create a plot for each value of the sequence coordinate. Allowing for
2311 # multiple cubes in a CubeList to be plotted in the same plot for similar
2312 # sequence values. Passing a CubeList into the internal plotting function
2313 # for similar values of the sequence coordinate. cube_slice can be an
2314 # iris.cube.Cube or an iris.cube.CubeList.
2316 for cube_slice in cube_iterables:
2317 # Normalize cube_slice to a list of cubes
2318 if isinstance(cube_slice, iris.cube.CubeList): 2318 ↛ 2319line 2318 didn't jump to line 2319 because the condition on line 2318 was never true
2319 cubes = list(cube_slice)
2320 elif isinstance(cube_slice, iris.cube.Cube): 2320 ↛ 2323line 2320 didn't jump to line 2323 because the condition on line 2320 was always true
2321 cubes = [cube_slice]
2322 else:
2323 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2325 # Use sequence value so multiple sequences can merge.
2326 seq_coord = cube_slice[0].coord(sequence_coordinate)
2327 plot_title, plot_filename = _set_title_and_filename(
2328 seq_coord, nplot, recipe_title, filename
2329 )
2331 # Format the coordinate value in a unit appropriate way.
2332 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2334 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2335 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2335 ↛ 2336line 2335 didn't jump to line 2336 because the condition on line 2335 was never true
2336 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2338 # Do the actual plotting.
2339 plotting_func(
2340 cube_slice,
2341 coords,
2342 stamp_coordinate,
2343 plot_filename,
2344 title,
2345 series_coordinate,
2346 )
2348 plot_index.append(plot_filename)
2349 else:
2350 # Format the title and filename using plotted series coordinate
2351 nplot = 1
2352 seq_coord = coords[0]
2353 plot_title, plot_filename = _set_title_and_filename(
2354 seq_coord, nplot, recipe_title, filename
2355 )
2357 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2358 if (
2359 "station" in [c.name() for c in cubes[0].coords()]
2360 and len(cubes[0].coord("station").points) > 1
2361 ):
2362 for station in cubes[0].coord("station").points:
2363 station_cubes = cubes.extract(iris.Constraint(station=station))
2364 station_name = station_cubes[0].coord("Station_Name").points[0]
2365 station_plotname = plot_filename.replace(
2366 ".png", "_" + station_name + ".png"
2367 )
2368 _plot_and_save_line_series(
2369 station_cubes,
2370 coords,
2371 "realization",
2372 station_plotname,
2373 f"{plot_title} {station_name}",
2374 )
2375 plot_index.append(station_plotname)
2377 else:
2378 # Do the actual plotting for all other series coordinate options.
2379 _plot_and_save_line_series(
2380 cubes, coords, stamp_coordinate, plot_filename, plot_title
2381 )
2383 plot_index.append(plot_filename)
2385 # append plot to list of plots
2386 complete_plot_index = _append_to_plot_index(plot_index)
2388 # Make a page to display the plots.
2389 _make_plot_html_page(complete_plot_index)
2391 return cube
2394def plot_vertical_line_series(
2395 cubes: iris.cube.Cube | iris.cube.CubeList,
2396 filename: str | None = None,
2397 series_coordinate: str = "model_level_number",
2398 sequence_coordinate: str = "time",
2399 # line_coordinate: str = "realization",
2400 **kwargs,
2401) -> iris.cube.Cube | iris.cube.CubeList:
2402 """Plot a line plot against a type of vertical coordinate.
2404 The Cube or CubeList must be 1D.
2406 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2407 then a sequence of plots will be produced.
2409 Parameters
2410 ----------
2411 iris.cube | iris.cube.CubeList
2412 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2413 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2414 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2415 filename: str, optional
2416 Name of the plot to write, used as a prefix for plot sequences. Defaults
2417 to the recipe name.
2418 series_coordinate: str, optional
2419 Coordinate to plot on the y-axis. Can be ``pressure`` or
2420 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2421 for LFRic. Defaults to ``model_level_number``.
2422 This coordinate must exist in the cube.
2423 sequence_coordinate: str, optional
2424 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2425 This coordinate must exist in the cube.
2427 Returns
2428 -------
2429 iris.cube.Cube | iris.cube.CubeList
2430 The original Cube or CubeList (so further operations can be applied).
2431 Plotted data.
2433 Raises
2434 ------
2435 ValueError
2436 If the cubes doesn't have the right dimensions.
2437 TypeError
2438 If the cube isn't a Cube or CubeList.
2439 """
2440 # Ensure we have a name for the plot file.
2441 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2443 cubes = iter_maybe(cubes)
2444 # Initialise empty list to hold all data from all cubes in a CubeList
2445 all_data = []
2447 # Store min/max ranges for x range.
2448 x_levels = []
2450 num_models = get_num_models(cubes)
2452 validate_cube_shape(cubes, num_models)
2454 # Iterate over all cubes in cube or CubeList and plot.
2455 coords = []
2456 for cube in cubes:
2457 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2458 try:
2459 coords.append(cube.coord(series_coordinate))
2460 except iris.exceptions.CoordinateNotFoundError as err:
2461 raise ValueError(
2462 f"Cube must have a {series_coordinate} coordinate."
2463 ) from err
2465 try:
2466 if cube.ndim > 1 or not cube.coords("realization"): 2466 ↛ 2474line 2466 didn't jump to line 2474 because the condition on line 2466 was always true
2467 cube.coord(sequence_coordinate)
2468 except iris.exceptions.CoordinateNotFoundError as err:
2469 raise ValueError(
2470 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2471 ) from err
2473 # Get minimum and maximum from levels information.
2474 _, levels, _ = colorbar_map_levels(cube, axis="x")
2475 if levels is not None: 2475 ↛ 2479line 2475 didn't jump to line 2479 because the condition on line 2475 was always true
2476 x_levels.append(min(levels))
2477 x_levels.append(max(levels))
2478 else:
2479 all_data.append(cube.data)
2481 if len(x_levels) == 0: 2481 ↛ 2483line 2481 didn't jump to line 2483 because the condition on line 2481 was never true
2482 # Combine all data into a single NumPy array
2483 combined_data = np.concatenate(all_data)
2485 # Set the lower and upper limit for the x-axis to ensure all plots have
2486 # same range. This needs to read the whole cube over the range of the
2487 # sequence and if applicable postage stamp coordinate.
2488 vmin = np.floor(combined_data.min())
2489 vmax = np.ceil(combined_data.max())
2490 else:
2491 vmin = min(x_levels)
2492 vmax = max(x_levels)
2494 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2495 # a single profile directly without iterating over a sequence.
2496 sequence_coords = [
2497 cube.coord(sequence_coordinate)
2498 for cube in cubes
2499 if cube.coords(sequence_coordinate)
2500 ]
2501 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2502 np.size(coord.points) > 1 for coord in sequence_coords
2503 )
2504 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2505 np.size(coord.points) == 1 for coord in sequence_coords
2506 )
2508 plot_index = []
2509 if has_sequence_coord: 2509 ↛ 2534line 2509 didn't jump to line 2534 because the condition on line 2509 was always true
2510 # Matching the slices (matching by seq coord point; it may happen that
2511 # evaluated models do not cover the same seq coord range, hence matching
2512 # necessary)
2513 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2514 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2515 for cubes_slice in cube_iterables:
2516 # Format the coordinate value in a unit appropriate way.
2517 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2518 plot_title, plot_filename = _set_title_and_filename(
2519 seq_coord, nplot, recipe_title, filename
2520 )
2522 # Do the actual plotting.
2523 _plot_and_save_vertical_line_series(
2524 cubes_slice,
2525 coords,
2526 "realization",
2527 plot_filename,
2528 series_coordinate,
2529 title=plot_title,
2530 vmin=vmin,
2531 vmax=vmax,
2532 )
2533 plot_index.append(plot_filename)
2534 elif has_scalar_sequence_coord:
2535 # Scalar sequence coordinate (typically aggregated time bounds):
2536 # make one plot and include sequence period in title/filename.
2537 plot_title, plot_filename = _set_title_and_filename(
2538 sequence_coords[0], 1, recipe_title, filename
2539 )
2541 _plot_and_save_vertical_line_series(
2542 cubes,
2543 coords,
2544 "realization",
2545 plot_filename,
2546 series_coordinate,
2547 title=plot_title,
2548 vmin=vmin,
2549 vmax=vmax,
2550 )
2551 plot_index.append(plot_filename)
2552 else:
2553 # 1D case: no sequence coordinate, plot a single profile.
2554 plot_title = recipe_title
2555 if filename:
2556 plot_filename = filename
2557 else:
2558 plot_filename = f"{slugify(plot_title)}.png"
2560 _plot_and_save_vertical_line_series(
2561 cubes,
2562 coords,
2563 "realization",
2564 plot_filename,
2565 series_coordinate,
2566 title=plot_title,
2567 vmin=vmin,
2568 vmax=vmax,
2569 )
2570 plot_index.append(plot_filename)
2572 # Add list of plots to plot metadata.
2573 complete_plot_index = _append_to_plot_index(plot_index)
2575 # Make a page to display the plots.
2576 _make_plot_html_page(complete_plot_index)
2578 return cubes
2581def qq_plot(
2582 cubes: iris.cube.CubeList,
2583 coordinates: list[str],
2584 percentiles: list[float],
2585 model_names: list[str],
2586 filename: str | None = None,
2587 one_to_one: bool = True,
2588 **kwargs,
2589) -> iris.cube.CubeList:
2590 """Plot a Quantile-Quantile plot between two models for common time points.
2592 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2593 collapsed within the operator over all specified coordinates such as
2594 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2595 ensemble members to ensure a 1D cube (array).
2597 Parameters
2598 ----------
2599 cubes: iris.cube.CubeList
2600 Two cubes of the same variable with different models.
2601 coordinate: list[str]
2602 The list of coordinates to collapse over. This list should be
2603 every coordinate within the cube to result in a 1D cube around
2604 the percentile coordinate.
2605 percent: list[float]
2606 A list of percentiles to appear in the plot.
2607 model_names: list[str]
2608 A list of model names to appear on the axis of the plot.
2609 filename: str, optional
2610 Filename of the plot to write.
2611 one_to_one: bool, optional
2612 If True a 1:1 line is plotted; if False it is not. Default is True.
2614 Raises
2615 ------
2616 ValueError
2617 When the cubes are not compatible.
2619 Notes
2620 -----
2621 The quantile-quantile plot is a variant on the scatter plot representing
2622 two datasets by their quantiles (percentiles) for common time points.
2623 This plot does not use a theoretical distribution to compare against, but
2624 compares percentiles of two datasets. This plot does
2625 not use all raw data points, but plots the selected percentiles (quantiles) of
2626 each variable instead for the two datasets, thereby normalising the data for a
2627 direct comparison between the selected percentiles of the two dataset distributions.
2629 Quantile-quantile plots are valuable for comparing against
2630 observations and other models. Identical percentiles between the variables
2631 will lie on the one-to-one line implying the values correspond well to each
2632 other. Where there is a deviation from the one-to-one line a range of
2633 possibilities exist depending on how and where the data is shifted (e.g.,
2634 Wilks 2011 [Wilks2011]_).
2636 For distributions above the one-to-one line the distribution is left-skewed;
2637 below is right-skewed. A distinct break implies a bimodal distribution, and
2638 closer values/values further apart at the tails imply poor representation of
2639 the extremes.
2641 """
2642 # Check cubes using same functionality as the difference operator.
2643 if len(cubes) != 2:
2644 raise ValueError("cubes should contain exactly 2 cubes.")
2645 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2646 other: Cube = cubes.extract_cube(
2647 iris.Constraint(
2648 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2649 )
2650 )
2652 # Get spatial coord names.
2653 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2654 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2656 # Ensure cubes to compare are on common differencing grid.
2657 # This is triggered if either
2658 # i) latitude and longitude shapes are not the same. Note grid points
2659 # are not compared directly as these can differ through rounding
2660 # errors.
2661 # ii) or variables are known to often sit on different grid staggering
2662 # in different models (e.g. cell center vs cell edge), as is the case
2663 # for UM and LFRic comparisons.
2664 # In future greater choice of regridding method might be applied depending
2665 # on variable type. Linear regridding can in general be appropriate for smooth
2666 # variables. Care should be taken with interpretation of differences
2667 # given this dependency on regridding.
2668 if (
2669 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2670 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2671 ) or (
2672 base.long_name
2673 in [
2674 "eastward_wind_at_10m",
2675 "northward_wind_at_10m",
2676 "northward_wind_at_cell_centres",
2677 "eastward_wind_at_cell_centres",
2678 "zonal_wind_at_pressure_levels",
2679 "meridional_wind_at_pressure_levels",
2680 "potential_vorticity_at_pressure_levels",
2681 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2682 ]
2683 ):
2684 logger.debug("Linear regridding base cube to other grid to compute differences")
2685 base = regrid_onto_cube(base, other, method="Linear")
2687 # Extract just common time points.
2688 base, other = _extract_common_time_points(base, other)
2690 # Equalise attributes so we can merge.
2691 fully_equalise_attributes([base, other])
2692 logger.debug("Base: %s\nOther: %s", base, other)
2694 # Collapse cubes.
2695 base = collapse(
2696 base,
2697 coordinate=coordinates,
2698 method="PERCENTILE",
2699 additional_percent=percentiles,
2700 )
2701 other = collapse(
2702 other,
2703 coordinate=coordinates,
2704 method="PERCENTILE",
2705 additional_percent=percentiles,
2706 )
2708 # Ensure we have a name for the plot file.
2709 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2710 title = f"{recipe_title}"
2712 if filename is None:
2713 filename = slugify(recipe_title)
2715 # Add file extension.
2716 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2718 # Do the actual plotting on a scatter plot
2719 _plot_and_save_scatter_plot(
2720 base, other, plot_filename, title, one_to_one, model_names
2721 )
2723 # Add list of plots to plot metadata.
2724 plot_index = _append_to_plot_index([plot_filename])
2726 # Make a page to display the plots.
2727 _make_plot_html_page(plot_index)
2729 return iris.cube.CubeList([base, other])
2732def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2733 """
2734 Plot a Hinton style triangle/scorecard plot.
2736 This plot type can be useful for summarising high level information, such as comparing
2737 how 'skillful' two models are when verified against observations for a variety of metrics,
2738 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2739 than customisable by the user as input arguments; many have been designed to automatically
2740 scale the plot depending on the number of x and y components.
2742 Parameters
2743 ----------
2744 change: np.ndarray
2745 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2746 size/direction.
2747 signif: np.ndarray
2748 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2749 xaxis_labels: list
2750 List of labels for the xaxis (must match the second dimension length of signif and change,
2751 along with magnitude if not None).
2752 yaxis_labels: list
2753 List of labels for the yaxis (must match the first dimension length of signif and change,
2754 along with magnitude if not None).
2755 magnitude: np.ndarray | None
2756 Optional 2D array, matching the shape of change, signif, which contains numerical values
2757 the user wishes to display under each respective triangle.
2759 Returns
2760 -------
2761 matplotlib axes object to either display or do further modifications to.
2762 """
2763 # Setup colors of triangles
2764 color_pos = "#7CAE00"
2765 color_neg = "#7B68EE"
2767 # Setup cell/text size ratios
2768 figsize = None
2769 cell_size_in = 0.35
2770 text_row_ratio = 0.25
2772 # Ensure arrays, and change to bool for sig.
2773 change = np.asarray(change)
2774 signif = np.asarray(signif).astype(bool)
2775 if magnitude is not None: 2775 ↛ 2776line 2775 didn't jump to line 2776 because the condition on line 2775 was never true
2776 magnitude = np.asarray(magnitude)
2778 # Get the number of x and y elements
2779 ny, nx = change.shape
2781 # Build non-uniform y coordinates
2782 tri_height = 1.0
2783 txt_height = text_row_ratio
2785 tri_y = []
2786 txt_y = []
2787 y_edges = [0.0]
2789 y = 0.0
2790 for _j in range(ny):
2791 tri_y.append(y + tri_height / 2)
2792 y += tri_height
2793 y_edges.append(y)
2795 if magnitude is not None: 2795 ↛ 2796line 2795 didn't jump to line 2796 because the condition on line 2795 was never true
2796 txt_y.append(y + txt_height / 2)
2797 y += txt_height
2798 y_edges.append(y)
2800 total_height = y
2802 # Dynamic figure size
2803 if figsize is None: 2803 ↛ 2808line 2803 didn't jump to line 2808 because the condition on line 2803 was always true
2804 width = nx * cell_size_in
2805 height = total_height * cell_size_in + 2
2806 figsize = (width, height)
2808 fig, ax = plt.subplots(figsize=figsize)
2810 # Setup axes and grid.
2811 ax.set_aspect("equal", adjustable="box")
2812 ax.set_xlim(-0.5, nx - 0.5)
2813 ax.set_ylim(0, total_height)
2815 ax.set_xticks(np.arange(nx))
2816 ax.set_xticklabels(xaxis_labels, rotation=90)
2818 ax.set_yticks(tri_y)
2819 ax.set_yticklabels(yaxis_labels)
2821 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2822 ax.set_yticks(y_edges, minor=True)
2824 ax.set_axisbelow(True)
2825 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2826 ax.grid(False, which="major")
2827 ax.tick_params(which="minor", length=0)
2829 ax.invert_yaxis()
2831 # Compute marker scaling (fixed overlap)
2832 fig.canvas.draw()
2834 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2835 width_in, height_in = bbox.width, bbox.height
2837 cell_w = (width_in * fig.dpi) / nx
2838 cell_h = (height_in * fig.dpi) / total_height
2839 cell_pixels = min(cell_w, cell_h)
2841 max_marker_size = (0.6 * cell_pixels) ** 2
2843 text_fontsize = cell_pixels * 0.15
2845 # Plot triangles + text
2846 for j in range(ny):
2847 for i in range(nx):
2848 val = change[j, i]
2849 if np.isnan(val): 2849 ↛ 2850line 2849 didn't jump to line 2850 because the condition on line 2849 was never true
2850 continue
2852 if abs(val) < 0.01: 2852 ↛ 2853line 2852 didn't jump to line 2853 because the condition on line 2852 was never true
2853 continue
2855 sig = signif[j, i]
2856 size = max_marker_size * abs(val)
2858 # Triangle style
2859 if val >= 0:
2860 marker = "^"
2861 color = color_pos
2862 else:
2863 marker = "v"
2864 color = color_neg
2866 if sig:
2867 edgecolor = "black"
2868 linewidth = 0.6
2869 else:
2870 edgecolor = "none"
2871 linewidth = 0.0
2873 # Triangle
2874 ax.scatter(
2875 i,
2876 tri_y[j],
2877 s=size,
2878 marker=marker,
2879 c=color,
2880 edgecolors=edgecolor,
2881 linewidths=linewidth,
2882 zorder=3,
2883 clip_on=True, # ensures no rendering bleed
2884 )
2886 # Text row
2887 if magnitude is not None: 2887 ↛ 2888line 2887 didn't jump to line 2888 because the condition on line 2887 was never true
2888 mag_val = magnitude[j, i]
2890 if not np.isnan(mag_val):
2891 ax.text(
2892 i,
2893 txt_y[j],
2894 f"{mag_val:.1f}",
2895 ha="center",
2896 va="center",
2897 fontsize=text_fontsize,
2898 color="black",
2899 zorder=4,
2900 )
2902 plt.tight_layout()
2903 return fig, ax
2906def scatter_plot(
2907 cube_x: iris.cube.Cube | iris.cube.CubeList,
2908 cube_y: iris.cube.Cube | iris.cube.CubeList,
2909 filename: str | None = None,
2910 one_to_one: bool = True,
2911 **kwargs,
2912) -> iris.cube.CubeList:
2913 """Plot a scatter plot between two variables.
2915 Both cubes must be 1D.
2917 Parameters
2918 ----------
2919 cube_x: Cube | CubeList
2920 1 dimensional Cube of the data to plot on y-axis.
2921 cube_y: Cube | CubeList
2922 1 dimensional Cube of the data to plot on x-axis.
2923 filename: str, optional
2924 Filename of the plot to write.
2925 one_to_one: bool, optional
2926 If True a 1:1 line is plotted; if False it is not. Default is True.
2928 Returns
2929 -------
2930 cubes: CubeList
2931 CubeList of the original x and y cubes for further processing.
2933 Raises
2934 ------
2935 ValueError
2936 If the cube doesn't have the right dimensions and cubes not the same
2937 size.
2938 TypeError
2939 If the cube isn't a single cube.
2941 Notes
2942 -----
2943 Scatter plots are used for determining if there is a relationship between
2944 two variables. Positive relations have a slope going from bottom left to top
2945 right; Negative relations have a slope going from top left to bottom right.
2946 """
2947 # Iterate over all cubes in cube or CubeList and plot.
2948 for cube_iter in iter_maybe(cube_x):
2949 # Check cubes are correct shape.
2950 cube_iter = check_single_cube(cube_iter)
2951 if cube_iter.ndim > 1:
2952 raise ValueError("cube_x must be 1D.")
2954 # Iterate over all cubes in cube or CubeList and plot.
2955 for cube_iter in iter_maybe(cube_y):
2956 # Check cubes are correct shape.
2957 cube_iter = check_single_cube(cube_iter)
2958 if cube_iter.ndim > 1:
2959 raise ValueError("cube_y must be 1D.")
2961 # Ensure we have a name for the plot file.
2962 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2963 title = f"{recipe_title}"
2965 if filename is None:
2966 filename = slugify(recipe_title)
2968 # Add file extension.
2969 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2971 # Do the actual plotting.
2972 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2974 # Add list of plots to plot metadata.
2975 plot_index = _append_to_plot_index([plot_filename])
2977 # Make a page to display the plots.
2978 _make_plot_html_page(plot_index)
2980 return iris.cube.CubeList([cube_x, cube_y])
2983def vector_plot(
2984 cube_u: iris.cube.Cube,
2985 cube_v: iris.cube.Cube,
2986 filename: str | None = None,
2987 sequence_coordinate: str = "time",
2988 **kwargs,
2989) -> iris.cube.CubeList:
2990 """Plot a vector plot based on the input u and v components."""
2991 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2993 # Cubes must have a matching sequence coordinate.
2994 try:
2995 # Check that the u and v cubes have the same sequence coordinate.
2996 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2996 ↛ anywhereline 2996 didn't jump anywhere: it always raised an exception.
2997 raise ValueError("Coordinates do not match.")
2998 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2999 raise ValueError(
3000 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3001 ) from err
3003 # Create a plot for each value of the sequence coordinate.
3004 plot_index = []
3005 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3006 for cube_u_slice, cube_v_slice in zip(
3007 cube_u.slices_over(sequence_coordinate),
3008 cube_v.slices_over(sequence_coordinate),
3009 strict=True,
3010 ):
3011 # Format the coordinate value in a unit appropriate way.
3012 seq_coord = cube_u_slice.coord(sequence_coordinate)
3013 plot_title, plot_filename = _set_title_and_filename(
3014 seq_coord, nplot, recipe_title, filename
3015 )
3017 # Do the actual plotting.
3018 _plot_and_save_vector_plot(
3019 cube_u_slice,
3020 cube_v_slice,
3021 filename=plot_filename,
3022 title=plot_title,
3023 method="pcolormesh",
3024 )
3025 plot_index.append(plot_filename)
3027 # Add list of plots to plot metadata.
3028 complete_plot_index = _append_to_plot_index(plot_index)
3030 # Make a page to display the plots.
3031 _make_plot_html_page(complete_plot_index)
3033 return iris.cube.CubeList([cube_u, cube_v])
3036def plot_histogram_series(
3037 cubes: iris.cube.Cube | iris.cube.CubeList,
3038 filename: str | None = None,
3039 sequence_coordinate: str = "time",
3040 stamp_coordinate: str = "realization",
3041 single_plot: bool = False,
3042 **kwargs,
3043) -> iris.cube.Cube | iris.cube.CubeList:
3044 """Plot a histogram plot for each vertical level provided.
3046 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3047 is present then a sequence of plots will be produced using the time slider
3048 functionality to scroll through histograms against time. If a
3049 stamp_coordinate is present then postage stamp plots will be produced. If
3050 stamp_coordinate and single_plot is True, all postage stamp plots will be
3051 plotted in a single plot instead of separate postage stamp plots.
3053 Parameters
3054 ----------
3055 cubes: Cube | iris.cube.CubeList
3056 Iris cube or CubeList of the data to plot. It should have a single dimension other
3057 than the stamp coordinate.
3058 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3059 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3060 filename: str, optional
3061 Name of the plot to write, used as a prefix for plot sequences. Defaults
3062 to the recipe name.
3063 sequence_coordinate: str, optional
3064 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3065 This coordinate must exist in the cube and will be used for the time
3066 slider.
3067 stamp_coordinate: str, optional
3068 Coordinate about which to plot postage stamp plots. Defaults to
3069 ``"realization"``.
3070 single_plot: bool, optional
3071 If True, all postage stamp plots will be plotted in a single plot. If
3072 False, each postage stamp plot will be plotted separately. Is only valid
3073 if stamp_coordinate exists and has more than a single point.
3075 Returns
3076 -------
3077 iris.cube.Cube | iris.cube.CubeList
3078 The original Cube or CubeList (so further operations can be applied).
3079 Plotted data.
3081 Raises
3082 ------
3083 ValueError
3084 If the cube doesn't have the right dimensions.
3085 TypeError
3086 If the cube isn't a Cube or CubeList.
3087 """
3088 recipe_title = get_recipe_metadata().get("title", "Histogram")
3090 cubes = iter_maybe(cubes)
3092 # Internal plotting function.
3093 plotting_func = _plot_and_save_histogram_series
3095 num_models = get_num_models(cubes)
3097 validate_cube_shape(cubes, num_models)
3099 # If several histograms are plotted, check sequence_coordinate
3100 check_sequence_coordinate(cubes, sequence_coordinate)
3102 # Get axis minimum and maximum from levels information.
3103 # If no levels set, derive minima and maxima from data in CubeList.
3104 vmin, vmax = _set_axis_range(cubes)
3106 # Make postage stamp plots if stamp_coordinate exists and has more than a
3107 # single point. If single_plot is True:
3108 # -- all postage stamp plots will be plotted in a single plot instead of
3109 # separate postage stamp plots.
3110 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3111 # produced per single model only
3112 if num_models == 1:
3113 if ( 3113 ↛ 3117line 3113 didn't jump to line 3117 because the condition on line 3113 was never true
3114 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3115 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3116 ):
3117 if single_plot:
3118 plotting_func = (
3119 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3120 )
3121 else:
3122 plotting_func = _plot_and_save_postage_stamp_histogram_series
3123 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3124 else:
3125 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3127 plot_index = []
3128 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3129 # Create a plot for each value of the sequence coordinate. Allowing for
3130 # multiple cubes in a CubeList to be plotted in the same plot for similar
3131 # sequence values. Passing a CubeList into the internal plotting function
3132 # for similar values of the sequence coordinate. cube_slice can be an
3133 # iris.cube.Cube or an iris.cube.CubeList.
3134 for cube_slice in cube_iterables:
3135 single_cube = cube_slice
3136 if isinstance(cube_slice, iris.cube.CubeList):
3137 single_cube = cube_slice[0]
3139 # Ensure valid stamp coordinate in cube dimensions
3140 if stamp_coordinate == "realization": 3140 ↛ 3143line 3140 didn't jump to line 3143 because the condition on line 3140 was always true
3141 stamp_coordinate = check_stamp_coordinate(single_cube)
3142 # Set plot titles and filename, based on sequence coordinate
3143 seq_coord = single_cube.coord(sequence_coordinate)
3144 # Use time coordinate in title and filename if single histogram output.
3145 if sequence_coordinate == "realization" and nplot == 1: 3145 ↛ 3146line 3145 didn't jump to line 3146 because the condition on line 3145 was never true
3146 seq_coord = single_cube.coord("time")
3147 # Use station name in title and filename if model vs obs comparison
3148 if sequence_coordinate == "station": 3148 ↛ 3149line 3148 didn't jump to line 3149 because the condition on line 3148 was never true
3149 seq_coord = single_cube.coord("Station_Name")
3151 plot_title, plot_filename = _set_title_and_filename(
3152 seq_coord, nplot, recipe_title, filename
3153 )
3155 # Do the actual plotting.
3156 plotting_func(
3157 cube_slice,
3158 filename=plot_filename,
3159 stamp_coordinate=stamp_coordinate,
3160 title=plot_title,
3161 vmin=vmin,
3162 vmax=vmax,
3163 )
3164 plot_index.append(plot_filename)
3166 # Add list of plots to plot metadata.
3167 complete_plot_index = _append_to_plot_index(plot_index)
3169 # Make a page to display the plots.
3170 _make_plot_html_page(complete_plot_index)
3172 return cubes
3175def plot_scatter_series(
3176 cubes: iris.cube.Cube | iris.cube.CubeList,
3177 filename: str | None = None,
3178 sequence_coordinate: str = "time",
3179 stamp_coordinate: str = "realization",
3180 hexbin: bool = False,
3181 **kwargs,
3182) -> iris.cube.Cube | iris.cube.CubeList:
3183 """Plot a scatter plot for each sequence coordinate provided.
3185 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3186 is present then a sequence of plots will be produced using the time slider
3187 functionality to scroll through scatter against time. If a
3188 stamp_coordinate is present then postage stamp plots will be produced. If
3189 stamp_coordinate and single_plot is True, all postage stamp plots will be
3190 plotted in a single plot instead of separate postage stamp plots.
3192 Parameters
3193 ----------
3194 cubes: Cube | iris.cube.CubeList
3195 Iris cube or CubeList of the data to plot. It should have a single dimension other
3196 than the stamp coordinate.
3197 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3198 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3199 filename: str, optional
3200 Name of the plot to write, used as a prefix for plot sequences. Defaults
3201 to the recipe name.
3202 sequence_coordinate: str, optional
3203 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3204 This coordinate must exist in the cube and will be used for the time
3205 slider.
3206 stamp_coordinate: str, optional
3207 Coordinate about which to plot postage stamp plots. Defaults to
3208 ``"realization"``.
3209 hexbin: bool, optional
3210 If True, generate hexbin comparison plot.
3211 If False, generate point-by-point scatter plot.
3213 Returns
3214 -------
3215 iris.cube.Cube | iris.cube.CubeList
3216 The original Cube or CubeList (so further operations can be applied).
3217 Plotted data.
3219 Raises
3220 ------
3221 ValueError
3222 If the cube doesn't have the right dimensions.
3223 TypeError
3224 If the cube isn't a Cube or CubeList.
3225 """
3226 recipe_title = get_recipe_metadata().get("title", "Scatter")
3228 cubes = iter_maybe(cubes)
3230 # Internal plotting function.
3231 plotting_func = _plot_and_save_scatter_series
3233 num_models = get_num_models(cubes)
3235 validate_cube_shape(cubes, num_models)
3237 check_sequence_coordinate(cubes, sequence_coordinate)
3239 vmin, vmax = _set_axis_range(cubes)
3241 # Require >1 models to compare on scatter plot
3242 if num_models > 1:
3243 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3244 else:
3245 raise ValueError(
3246 "Scatter plot series requires multiple number of models in input data."
3247 )
3249 plot_index = []
3250 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3251 # Create a plot for each value of the sequence coordinate. Allowing for
3252 # multiple cubes in a CubeList to be plotted in the same plot for similar
3253 # sequence values. Passing a CubeList into the internal plotting function
3254 # for similar values of the sequence coordinate. cube_slice can be an
3255 # iris.cube.Cube or an iris.cube.CubeList.
3256 for cube_slice in cube_iterables:
3257 single_cube = cube_slice
3258 if isinstance(cube_slice, iris.cube.CubeList): 3258 ↛ 3262line 3258 didn't jump to line 3262 because the condition on line 3258 was always true
3259 single_cube = cube_slice[0]
3261 # Ensure valid stamp coordinate in cube dimensions
3262 if stamp_coordinate == "realization": 3262 ↛ 3265line 3262 didn't jump to line 3265 because the condition on line 3262 was always true
3263 stamp_coordinate = check_stamp_coordinate(single_cube)
3264 # Set plot titles and filename, based on sequence coordinate
3265 seq_coord = single_cube.coord(sequence_coordinate)
3266 # Use time coordinate in title and filename if single histogram output.
3267 if sequence_coordinate == "realization" and nplot == 1:
3268 seq_coord = single_cube.coord("time")
3269 # Use station name in title and filename if model vs obs comparison
3270 if sequence_coordinate == "station":
3271 seq_coord = single_cube.coord("Station_Name")
3273 plot_title, plot_filename = _set_title_and_filename(
3274 seq_coord, nplot, recipe_title, filename
3275 )
3277 # Do the actual plotting.
3278 plotting_func(
3279 cube_slice,
3280 filename=plot_filename,
3281 stamp_coordinate=stamp_coordinate,
3282 title=plot_title,
3283 vmin=vmin,
3284 vmax=vmax,
3285 hexbin=hexbin,
3286 )
3287 plot_index.append(plot_filename)
3289 # Add list of plots to plot metadata.
3290 complete_plot_index = _append_to_plot_index(plot_index)
3292 # Make a page to display the plots.
3293 _make_plot_html_page(complete_plot_index)
3295 return cubes
3298def _plot_and_save_postage_stamp_power_spectrum_series(
3299 cubes: iris.cube.Cube,
3300 coords: list[iris.coords.Coord],
3301 stamp_coordinate: str,
3302 filename: str,
3303 title: str,
3304 series_coordinate: str | None = None,
3305 **kwargs,
3306):
3307 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3309 Parameters
3310 ----------
3311 cubes: Cube or CubeList
3312 Cube or Cubelist of the power spectrum data.
3313 coords: list[Coord]
3314 Coordinates to plot on the x-axis, one per cube.
3315 stamp_coordinate: str
3316 Coordinate that becomes different plots.
3317 filename: str
3318 Filename of the plot to write.
3319 title: str
3320 Plot title.
3321 series_coordinate: str, optional
3322 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3324 """
3325 # Use the smallest square grid that will fit the members.
3326 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3328 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3329 model_colors_map = get_model_colors_map(cubes)
3330 # ax = plt.gca()
3331 # Make a subplot for each member.
3332 for member, subplot in zip(
3333 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3334 ):
3335 ax = plt.subplot(grid_size, grid_size, subplot)
3337 # Store min/max ranges.
3338 y_levels = []
3340 line_marker = None
3341 line_width = 1
3343 for cube in iter_maybe(member):
3344 xcoord = _select_series_coord(cube, series_coordinate)
3345 xname = xcoord.points
3347 yfield = cube.data # power spectrum
3348 label = None
3349 color = "black"
3350 if model_colors_map: 3350 ↛ 3351line 3350 didn't jump to line 3351 because the condition on line 3350 was never true
3351 label = cube.attributes.get("model_name")
3352 color = model_colors_map.get(label)
3354 if member.coord(stamp_coordinate).points == [0]:
3355 ax.plot(
3356 xname,
3357 yfield,
3358 color=color,
3359 marker=line_marker,
3360 ls="-",
3361 lw=line_width,
3362 label=f"{label} (control)"
3363 if len(cube.coord(stamp_coordinate).points) > 1
3364 else label,
3365 )
3366 # Label with member if part of an ensemble and not the control.
3367 else:
3368 ax.plot(
3369 xname,
3370 yfield,
3371 color=color,
3372 ls="-",
3373 lw=1.5,
3374 alpha=0.75,
3375 label=f"{label} (member)",
3376 )
3378 # Calculate the global min/max if multiple cubes are given.
3379 _, levels, _ = colorbar_map_levels(cube, axis="y")
3380 if levels is not None: 3380 ↛ 3381line 3380 didn't jump to line 3381 because the condition on line 3380 was never true
3381 y_levels.append(min(levels))
3382 y_levels.append(max(levels))
3384 # Add some labels and tweak the style.
3385 title = f"{title}"
3386 ax.set_title(title, fontsize=16)
3388 # Set appropriate x-axis label based on coordinate
3389 if series_coordinate == "wavelength" or ( 3389 ↛ 3392line 3389 didn't jump to line 3392 because the condition on line 3389 was never true
3390 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3391 ):
3392 ax.set_xlabel("Wavelength (km)", fontsize=14)
3393 elif series_coordinate == "physical_wavenumber" or ( 3393 ↛ 3398line 3393 didn't jump to line 3398 because the condition on line 3393 was always true
3394 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3395 ):
3396 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3397 else: # frequency or check units
3398 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3399 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3400 else:
3401 ax.set_xlabel("Wavenumber", fontsize=14)
3403 ax.set_ylabel("Power Spectral Density", fontsize=14)
3404 ax.tick_params(axis="both", labelsize=12)
3406 # Set log-log scale
3407 ax.set_xscale("log")
3408 ax.set_yscale("log")
3410 # Add gridlines
3411 ax.grid(linestyle="--", color="grey", linewidth=1)
3412 # Ientify unique labels for legend
3413 handles = list(
3414 {
3415 label: handle
3416 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3417 }.values()
3418 )
3419 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3421 ax = plt.gca()
3422 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3424 # Save plot.
3425 _save_close_figure(fig, "histogram postage stamp", filename)
3428def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3429 cubes: iris.cube.Cube,
3430 coords: list[iris.coords.Coord],
3431 stamp_coordinate: str,
3432 filename: str,
3433 title: str,
3434 series_coordinate: str | None = None,
3435 **kwargs,
3436):
3437 """Plot and save power spectra for ensemble members in single plot.
3439 Parameters
3440 ----------
3441 cubes: Cube or CubeList
3442 Cube or Cubelist of the power spectrum data.
3443 coords: list[Coord]
3444 Coordinates to plot on the x-axis, one per cube.
3445 stamp_coordinate: str
3446 Coordinate that becomes different plots.
3447 filename: str
3448 Filename of the plot to write.
3449 title: str
3450 Plot title.
3451 series_coordinate: str, optional
3452 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3454 """
3455 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3456 model_colors_map = get_model_colors_map(cubes)
3458 line_marker = None
3459 line_width = 1
3461 # Compute ensemble statistics to show spread
3462 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3463 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3464 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3466 xcoord_global = mean_cube.coord(series_coordinate)
3467 x_global = xcoord_global.points
3469 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3470 xcoord = _select_series_coord(member, series_coordinate)
3471 xname = xcoord.points
3473 yfield = member.data # power spectrum
3474 color = "black"
3475 if model_colors_map: 3475 ↛ 3479line 3475 didn't jump to line 3479 because the condition on line 3475 was always true
3476 label = member.attributes.get("model_name") if i == 0 else None
3477 color = model_colors_map.get(label)
3479 if member.coord(stamp_coordinate).points == [0]:
3480 ax.plot(
3481 xname,
3482 yfield,
3483 color=color,
3484 marker=line_marker,
3485 ls="-",
3486 lw=line_width,
3487 label=f"{label} (control)"
3488 if len(member.coord(stamp_coordinate).points) > 1
3489 else label,
3490 )
3491 # Label with member number if part of an ensemble and not the control.
3492 else:
3493 ax.plot(
3494 xname,
3495 yfield,
3496 color=color,
3497 ls="-",
3498 lw=1.5,
3499 alpha=0.75,
3500 label=label,
3501 )
3503 # Set appropriate x-axis label based on coordinate
3504 if series_coordinate == "wavelength" or ( 3504 ↛ 3507line 3504 didn't jump to line 3507 because the condition on line 3504 was never true
3505 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3506 ):
3507 ax.set_xlabel("Wavelength (km)", fontsize=14)
3508 elif series_coordinate == "physical_wavenumber" or ( 3508 ↛ 3513line 3508 didn't jump to line 3513 because the condition on line 3508 was always true
3509 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3510 ):
3511 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3512 else: # frequency or check units
3513 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3514 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3515 else:
3516 ax.set_xlabel("Wavenumber", fontsize=14)
3518 # Add ensemble spread shading
3519 ax.fill_between(
3520 x_global,
3521 min_cube.data,
3522 max_cube.data,
3523 color="grey",
3524 alpha=0.3,
3525 label="Ensemble spread",
3526 )
3528 # Add ensemble mean line
3529 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3531 ax.set_ylabel("Power Spectral Density", fontsize=14)
3532 ax.tick_params(axis="both", labelsize=12)
3534 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3535 # Set log-log scale
3536 ax.set_xscale("log")
3537 ax.set_yscale("log")
3539 # Add gridlines
3540 ax.grid(linestyle="--", color="grey", linewidth=1)
3541 # Identify unique labels for legend
3542 handles = list(
3543 {
3544 label: handle
3545 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3546 }.values()
3547 )
3548 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3550 # Figure title.
3551 ax.set_title(title, fontsize=16)
3553 # Save plot.
3554 _save_close_figure(fig, "power spectra postage stamp", filename)