Coverage for src/CSET/operators/plot.py: 83%
1102 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-10 15:54 +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 if model_cube.coords("realization") and model_cube.ndim > 2:
2230 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2232 plot_index = []
2234 # Check if this is a spectral plot by looking for spectral coordinates
2235 is_spectral_plot = series_coordinate in [
2236 "frequency",
2237 "physical_wavenumber",
2238 "wavelength",
2239 ]
2241 if is_spectral_plot:
2242 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2243 # coordinate frequency/wavenumber.
2244 # If several power spectra are plotted with time as sequence_coordinate for the
2245 # time slider option.
2247 # Internal plotting function.
2248 plotting_func = _plot_and_save_line_power_spectrum_series
2250 for model_cube in cubes:
2251 try:
2252 model_cube.coord(sequence_coordinate)
2253 except iris.exceptions.CoordinateNotFoundError as err:
2254 raise ValueError(
2255 f"Cube must have a {sequence_coordinate} coordinate."
2256 ) from err
2258 if num_models == 1: 2258 ↛ 2273line 2258 didn't jump to line 2273 because the condition on line 2258 was always true
2259 # check for ensembles
2260 if ( 2260 ↛ 2264line 2260 didn't jump to line 2264 because the condition on line 2260 was never true
2261 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2262 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2263 ):
2264 if single_plot:
2265 # Plot spectra, mean and ensemble spread on 1 plot
2266 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2267 else:
2268 # Plot postage stamps
2269 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2270 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2271 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2272 else:
2273 all_points = sorted(
2274 set(
2275 itertools.chain.from_iterable(
2276 cb.coord(sequence_coordinate).points for cb in cubes
2277 )
2278 )
2279 )
2280 all_slices = list(
2281 itertools.chain.from_iterable(
2282 cb.slices_over(sequence_coordinate) for cb in cubes
2283 )
2284 )
2285 # Matched slices (matched by seq coord point; it may happen that
2286 # evaluated models do not cover the same seq coord range, hence matching
2287 # necessary)
2288 cube_iterables = [
2289 iris.cube.CubeList(
2290 s
2291 for s in all_slices
2292 if s.coord(sequence_coordinate).points[0] == point
2293 )
2294 for point in all_points
2295 ]
2296 nplot = len(all_points)
2298 # Create a plot for each value of the sequence coordinate. Allowing for
2299 # multiple cubes in a CubeList to be plotted in the same plot for similar
2300 # sequence values. Passing a CubeList into the internal plotting function
2301 # for similar values of the sequence coordinate. cube_slice can be an
2302 # iris.cube.Cube or an iris.cube.CubeList.
2304 for cube_slice in cube_iterables:
2305 # Normalize cube_slice to a list of cubes
2306 if isinstance(cube_slice, iris.cube.CubeList): 2306 ↛ 2307line 2306 didn't jump to line 2307 because the condition on line 2306 was never true
2307 cubes = list(cube_slice)
2308 elif isinstance(cube_slice, iris.cube.Cube): 2308 ↛ 2311line 2308 didn't jump to line 2311 because the condition on line 2308 was always true
2309 cubes = [cube_slice]
2310 else:
2311 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2313 # Use sequence value so multiple sequences can merge.
2314 seq_coord = cube_slice[0].coord(sequence_coordinate)
2315 plot_title, plot_filename = _set_title_and_filename(
2316 seq_coord, nplot, recipe_title, filename
2317 )
2319 # Format the coordinate value in a unit appropriate way.
2320 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2322 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2323 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2323 ↛ 2324line 2323 didn't jump to line 2324 because the condition on line 2323 was never true
2324 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2326 # Do the actual plotting.
2327 plotting_func(
2328 cube_slice,
2329 coords,
2330 stamp_coordinate,
2331 plot_filename,
2332 title,
2333 series_coordinate,
2334 )
2336 plot_index.append(plot_filename)
2337 else:
2338 # Format the title and filename using plotted series coordinate
2339 nplot = 1
2340 seq_coord = coords[0]
2341 plot_title, plot_filename = _set_title_and_filename(
2342 seq_coord, nplot, recipe_title, filename
2343 )
2345 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2346 if (
2347 "station" in [c.name() for c in cubes[0].coords()]
2348 and len(cubes[0].coord("station").points) > 1
2349 ):
2350 for station in cubes[0].coord("station").points:
2351 station_cubes = cubes.extract(iris.Constraint(station=station))
2352 station_name = station_cubes[0].coord("Station_Name").points[0]
2353 station_plotname = plot_filename.replace(
2354 ".png", "_" + station_name + ".png"
2355 )
2356 _plot_and_save_line_series(
2357 station_cubes,
2358 coords,
2359 "realization",
2360 station_plotname,
2361 f"{plot_title} {station_name}",
2362 )
2363 plot_index.append(station_plotname)
2365 else:
2366 # Do the actual plotting for all other series coordinate options.
2367 _plot_and_save_line_series(
2368 cubes, coords, stamp_coordinate, plot_filename, plot_title
2369 )
2371 plot_index.append(plot_filename)
2373 # append plot to list of plots
2374 complete_plot_index = _append_to_plot_index(plot_index)
2376 # Make a page to display the plots.
2377 _make_plot_html_page(complete_plot_index)
2379 return cube
2382def plot_vertical_line_series(
2383 cubes: iris.cube.Cube | iris.cube.CubeList,
2384 filename: str | None = None,
2385 series_coordinate: str = "model_level_number",
2386 sequence_coordinate: str = "time",
2387 # line_coordinate: str = "realization",
2388 **kwargs,
2389) -> iris.cube.Cube | iris.cube.CubeList:
2390 """Plot a line plot against a type of vertical coordinate.
2392 The Cube or CubeList must be 1D.
2394 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2395 then a sequence of plots will be produced.
2397 Parameters
2398 ----------
2399 iris.cube | iris.cube.CubeList
2400 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2401 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2402 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2403 filename: str, optional
2404 Name of the plot to write, used as a prefix for plot sequences. Defaults
2405 to the recipe name.
2406 series_coordinate: str, optional
2407 Coordinate to plot on the y-axis. Can be ``pressure`` or
2408 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2409 for LFRic. Defaults to ``model_level_number``.
2410 This coordinate must exist in the cube.
2411 sequence_coordinate: str, optional
2412 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2413 This coordinate must exist in the cube.
2415 Returns
2416 -------
2417 iris.cube.Cube | iris.cube.CubeList
2418 The original Cube or CubeList (so further operations can be applied).
2419 Plotted data.
2421 Raises
2422 ------
2423 ValueError
2424 If the cubes doesn't have the right dimensions.
2425 TypeError
2426 If the cube isn't a Cube or CubeList.
2427 """
2428 # Ensure we have a name for the plot file.
2429 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2431 cubes = iter_maybe(cubes)
2432 # Initialise empty list to hold all data from all cubes in a CubeList
2433 all_data = []
2435 # Store min/max ranges for x range.
2436 x_levels = []
2438 num_models = get_num_models(cubes)
2440 validate_cube_shape(cubes, num_models)
2442 # Iterate over all cubes in cube or CubeList and plot.
2443 coords = []
2444 for cube in cubes:
2445 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2446 try:
2447 coords.append(cube.coord(series_coordinate))
2448 except iris.exceptions.CoordinateNotFoundError as err:
2449 raise ValueError(
2450 f"Cube must have a {series_coordinate} coordinate."
2451 ) from err
2453 try:
2454 if cube.ndim > 1 or not cube.coords("realization"): 2454 ↛ 2462line 2454 didn't jump to line 2462 because the condition on line 2454 was always true
2455 cube.coord(sequence_coordinate)
2456 except iris.exceptions.CoordinateNotFoundError as err:
2457 raise ValueError(
2458 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2459 ) from err
2461 # Get minimum and maximum from levels information.
2462 _, levels, _ = colorbar_map_levels(cube, axis="x")
2463 if levels is not None: 2463 ↛ 2467line 2463 didn't jump to line 2467 because the condition on line 2463 was always true
2464 x_levels.append(min(levels))
2465 x_levels.append(max(levels))
2466 else:
2467 all_data.append(cube.data)
2469 if len(x_levels) == 0: 2469 ↛ 2471line 2469 didn't jump to line 2471 because the condition on line 2469 was never true
2470 # Combine all data into a single NumPy array
2471 combined_data = np.concatenate(all_data)
2473 # Set the lower and upper limit for the x-axis to ensure all plots have
2474 # same range. This needs to read the whole cube over the range of the
2475 # sequence and if applicable postage stamp coordinate.
2476 vmin = np.floor(combined_data.min())
2477 vmax = np.ceil(combined_data.max())
2478 else:
2479 vmin = min(x_levels)
2480 vmax = max(x_levels)
2482 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2483 # a single profile directly without iterating over a sequence.
2484 sequence_coords = [
2485 cube.coord(sequence_coordinate)
2486 for cube in cubes
2487 if cube.coords(sequence_coordinate)
2488 ]
2489 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2490 np.size(coord.points) > 1 for coord in sequence_coords
2491 )
2492 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2493 np.size(coord.points) == 1 for coord in sequence_coords
2494 )
2496 plot_index = []
2497 if has_sequence_coord: 2497 ↛ 2522line 2497 didn't jump to line 2522 because the condition on line 2497 was always true
2498 # Matching the slices (matching by seq coord point; it may happen that
2499 # evaluated models do not cover the same seq coord range, hence matching
2500 # necessary)
2501 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2502 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2503 for cubes_slice in cube_iterables:
2504 # Format the coordinate value in a unit appropriate way.
2505 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2506 plot_title, plot_filename = _set_title_and_filename(
2507 seq_coord, nplot, recipe_title, filename
2508 )
2510 # Do the actual plotting.
2511 _plot_and_save_vertical_line_series(
2512 cubes_slice,
2513 coords,
2514 "realization",
2515 plot_filename,
2516 series_coordinate,
2517 title=plot_title,
2518 vmin=vmin,
2519 vmax=vmax,
2520 )
2521 plot_index.append(plot_filename)
2522 elif has_scalar_sequence_coord:
2523 # Scalar sequence coordinate (typically aggregated time bounds):
2524 # make one plot and include sequence period in title/filename.
2525 plot_title, plot_filename = _set_title_and_filename(
2526 sequence_coords[0], 1, recipe_title, filename
2527 )
2529 _plot_and_save_vertical_line_series(
2530 cubes,
2531 coords,
2532 "realization",
2533 plot_filename,
2534 series_coordinate,
2535 title=plot_title,
2536 vmin=vmin,
2537 vmax=vmax,
2538 )
2539 plot_index.append(plot_filename)
2540 else:
2541 # 1D case: no sequence coordinate, plot a single profile.
2542 plot_title = recipe_title
2543 if filename:
2544 plot_filename = filename
2545 else:
2546 plot_filename = f"{slugify(plot_title)}.png"
2548 _plot_and_save_vertical_line_series(
2549 cubes,
2550 coords,
2551 "realization",
2552 plot_filename,
2553 series_coordinate,
2554 title=plot_title,
2555 vmin=vmin,
2556 vmax=vmax,
2557 )
2558 plot_index.append(plot_filename)
2560 # Add list of plots to plot metadata.
2561 complete_plot_index = _append_to_plot_index(plot_index)
2563 # Make a page to display the plots.
2564 _make_plot_html_page(complete_plot_index)
2566 return cubes
2569def qq_plot(
2570 cubes: iris.cube.CubeList,
2571 coordinates: list[str],
2572 percentiles: list[float],
2573 model_names: list[str],
2574 filename: str | None = None,
2575 one_to_one: bool = True,
2576 **kwargs,
2577) -> iris.cube.CubeList:
2578 """Plot a Quantile-Quantile plot between two models for common time points.
2580 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2581 collapsed within the operator over all specified coordinates such as
2582 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2583 ensemble members to ensure a 1D cube (array).
2585 Parameters
2586 ----------
2587 cubes: iris.cube.CubeList
2588 Two cubes of the same variable with different models.
2589 coordinate: list[str]
2590 The list of coordinates to collapse over. This list should be
2591 every coordinate within the cube to result in a 1D cube around
2592 the percentile coordinate.
2593 percent: list[float]
2594 A list of percentiles to appear in the plot.
2595 model_names: list[str]
2596 A list of model names to appear on the axis of the plot.
2597 filename: str, optional
2598 Filename of the plot to write.
2599 one_to_one: bool, optional
2600 If True a 1:1 line is plotted; if False it is not. Default is True.
2602 Raises
2603 ------
2604 ValueError
2605 When the cubes are not compatible.
2607 Notes
2608 -----
2609 The quantile-quantile plot is a variant on the scatter plot representing
2610 two datasets by their quantiles (percentiles) for common time points.
2611 This plot does not use a theoretical distribution to compare against, but
2612 compares percentiles of two datasets. This plot does
2613 not use all raw data points, but plots the selected percentiles (quantiles) of
2614 each variable instead for the two datasets, thereby normalising the data for a
2615 direct comparison between the selected percentiles of the two dataset distributions.
2617 Quantile-quantile plots are valuable for comparing against
2618 observations and other models. Identical percentiles between the variables
2619 will lie on the one-to-one line implying the values correspond well to each
2620 other. Where there is a deviation from the one-to-one line a range of
2621 possibilities exist depending on how and where the data is shifted (e.g.,
2622 Wilks 2011 [Wilks2011]_).
2624 For distributions above the one-to-one line the distribution is left-skewed;
2625 below is right-skewed. A distinct break implies a bimodal distribution, and
2626 closer values/values further apart at the tails imply poor representation of
2627 the extremes.
2629 """
2630 # Check cubes using same functionality as the difference operator.
2631 if len(cubes) != 2:
2632 raise ValueError("cubes should contain exactly 2 cubes.")
2633 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2634 other: Cube = cubes.extract_cube(
2635 iris.Constraint(
2636 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2637 )
2638 )
2640 # Get spatial coord names.
2641 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2642 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2644 # Ensure cubes to compare are on common differencing grid.
2645 # This is triggered if either
2646 # i) latitude and longitude shapes are not the same. Note grid points
2647 # are not compared directly as these can differ through rounding
2648 # errors.
2649 # ii) or variables are known to often sit on different grid staggering
2650 # in different models (e.g. cell center vs cell edge), as is the case
2651 # for UM and LFRic comparisons.
2652 # In future greater choice of regridding method might be applied depending
2653 # on variable type. Linear regridding can in general be appropriate for smooth
2654 # variables. Care should be taken with interpretation of differences
2655 # given this dependency on regridding.
2656 if (
2657 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2658 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2659 ) or (
2660 base.long_name
2661 in [
2662 "eastward_wind_at_10m",
2663 "northward_wind_at_10m",
2664 "northward_wind_at_cell_centres",
2665 "eastward_wind_at_cell_centres",
2666 "zonal_wind_at_pressure_levels",
2667 "meridional_wind_at_pressure_levels",
2668 "potential_vorticity_at_pressure_levels",
2669 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2670 ]
2671 ):
2672 logger.debug("Linear regridding base cube to other grid to compute differences")
2673 base = regrid_onto_cube(base, other, method="Linear")
2675 # Extract just common time points.
2676 base, other = _extract_common_time_points(base, other)
2678 # Equalise attributes so we can merge.
2679 fully_equalise_attributes([base, other])
2680 logger.debug("Base: %s\nOther: %s", base, other)
2682 # Collapse cubes.
2683 base = collapse(
2684 base,
2685 coordinate=coordinates,
2686 method="PERCENTILE",
2687 additional_percent=percentiles,
2688 )
2689 other = collapse(
2690 other,
2691 coordinate=coordinates,
2692 method="PERCENTILE",
2693 additional_percent=percentiles,
2694 )
2696 # Ensure we have a name for the plot file.
2697 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2698 title = f"{recipe_title}"
2700 if filename is None:
2701 filename = slugify(recipe_title)
2703 # Add file extension.
2704 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2706 # Do the actual plotting on a scatter plot
2707 _plot_and_save_scatter_plot(
2708 base, other, plot_filename, title, one_to_one, model_names
2709 )
2711 # Add list of plots to plot metadata.
2712 plot_index = _append_to_plot_index([plot_filename])
2714 # Make a page to display the plots.
2715 _make_plot_html_page(plot_index)
2717 return iris.cube.CubeList([base, other])
2720def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2721 """
2722 Plot a Hinton style triangle/scorecard plot.
2724 This plot type can be useful for summarising high level information, such as comparing
2725 how 'skillful' two models are when verified against observations for a variety of metrics,
2726 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2727 than customisable by the user as input arguments; many have been designed to automatically
2728 scale the plot depending on the number of x and y components.
2730 Parameters
2731 ----------
2732 change: np.ndarray
2733 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2734 size/direction.
2735 signif: np.ndarray
2736 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2737 xaxis_labels: list
2738 List of labels for the xaxis (must match the second dimension length of signif and change,
2739 along with magnitude if not None).
2740 yaxis_labels: list
2741 List of labels for the yaxis (must match the first dimension length of signif and change,
2742 along with magnitude if not None).
2743 magnitude: np.ndarray | None
2744 Optional 2D array, matching the shape of change, signif, which contains numerical values
2745 the user wishes to display under each respective triangle.
2747 Returns
2748 -------
2749 matplotlib axes object to either display or do further modifications to.
2750 """
2751 # Setup colors of triangles
2752 color_pos = "#7CAE00"
2753 color_neg = "#7B68EE"
2755 # Setup cell/text size ratios
2756 figsize = None
2757 cell_size_in = 0.35
2758 text_row_ratio = 0.25
2760 # Ensure arrays, and change to bool for sig.
2761 change = np.asarray(change)
2762 signif = np.asarray(signif).astype(bool)
2763 if magnitude is not None: 2763 ↛ 2764line 2763 didn't jump to line 2764 because the condition on line 2763 was never true
2764 magnitude = np.asarray(magnitude)
2766 # Get the number of x and y elements
2767 ny, nx = change.shape
2769 # Build non-uniform y coordinates
2770 tri_height = 1.0
2771 txt_height = text_row_ratio
2773 tri_y = []
2774 txt_y = []
2775 y_edges = [0.0]
2777 y = 0.0
2778 for _j in range(ny):
2779 tri_y.append(y + tri_height / 2)
2780 y += tri_height
2781 y_edges.append(y)
2783 if magnitude is not None: 2783 ↛ 2784line 2783 didn't jump to line 2784 because the condition on line 2783 was never true
2784 txt_y.append(y + txt_height / 2)
2785 y += txt_height
2786 y_edges.append(y)
2788 total_height = y
2790 # Dynamic figure size
2791 if figsize is None: 2791 ↛ 2796line 2791 didn't jump to line 2796 because the condition on line 2791 was always true
2792 width = nx * cell_size_in
2793 height = total_height * cell_size_in + 2
2794 figsize = (width, height)
2796 fig, ax = plt.subplots(figsize=figsize)
2798 # Setup axes and grid.
2799 ax.set_aspect("equal", adjustable="box")
2800 ax.set_xlim(-0.5, nx - 0.5)
2801 ax.set_ylim(0, total_height)
2803 ax.set_xticks(np.arange(nx))
2804 ax.set_xticklabels(xaxis_labels, rotation=90)
2806 ax.set_yticks(tri_y)
2807 ax.set_yticklabels(yaxis_labels)
2809 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2810 ax.set_yticks(y_edges, minor=True)
2812 ax.set_axisbelow(True)
2813 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2814 ax.grid(False, which="major")
2815 ax.tick_params(which="minor", length=0)
2817 ax.invert_yaxis()
2819 # Compute marker scaling (fixed overlap)
2820 fig.canvas.draw()
2822 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2823 width_in, height_in = bbox.width, bbox.height
2825 cell_w = (width_in * fig.dpi) / nx
2826 cell_h = (height_in * fig.dpi) / total_height
2827 cell_pixels = min(cell_w, cell_h)
2829 max_marker_size = (0.6 * cell_pixels) ** 2
2831 text_fontsize = cell_pixels * 0.15
2833 # Plot triangles + text
2834 for j in range(ny):
2835 for i in range(nx):
2836 val = change[j, i]
2837 if np.isnan(val): 2837 ↛ 2838line 2837 didn't jump to line 2838 because the condition on line 2837 was never true
2838 continue
2840 if abs(val) < 0.01: 2840 ↛ 2841line 2840 didn't jump to line 2841 because the condition on line 2840 was never true
2841 continue
2843 sig = signif[j, i]
2844 size = max_marker_size * abs(val)
2846 # Triangle style
2847 if val >= 0:
2848 marker = "^"
2849 color = color_pos
2850 else:
2851 marker = "v"
2852 color = color_neg
2854 if sig:
2855 edgecolor = "black"
2856 linewidth = 0.6
2857 else:
2858 edgecolor = "none"
2859 linewidth = 0.0
2861 # Triangle
2862 ax.scatter(
2863 i,
2864 tri_y[j],
2865 s=size,
2866 marker=marker,
2867 c=color,
2868 edgecolors=edgecolor,
2869 linewidths=linewidth,
2870 zorder=3,
2871 clip_on=True, # ensures no rendering bleed
2872 )
2874 # Text row
2875 if magnitude is not None: 2875 ↛ 2876line 2875 didn't jump to line 2876 because the condition on line 2875 was never true
2876 mag_val = magnitude[j, i]
2878 if not np.isnan(mag_val):
2879 ax.text(
2880 i,
2881 txt_y[j],
2882 f"{mag_val:.1f}",
2883 ha="center",
2884 va="center",
2885 fontsize=text_fontsize,
2886 color="black",
2887 zorder=4,
2888 )
2890 plt.tight_layout()
2891 return fig, ax
2894def scatter_plot(
2895 cube_x: iris.cube.Cube | iris.cube.CubeList,
2896 cube_y: iris.cube.Cube | iris.cube.CubeList,
2897 filename: str | None = None,
2898 one_to_one: bool = True,
2899 **kwargs,
2900) -> iris.cube.CubeList:
2901 """Plot a scatter plot between two variables.
2903 Both cubes must be 1D.
2905 Parameters
2906 ----------
2907 cube_x: Cube | CubeList
2908 1 dimensional Cube of the data to plot on y-axis.
2909 cube_y: Cube | CubeList
2910 1 dimensional Cube of the data to plot on x-axis.
2911 filename: str, optional
2912 Filename of the plot to write.
2913 one_to_one: bool, optional
2914 If True a 1:1 line is plotted; if False it is not. Default is True.
2916 Returns
2917 -------
2918 cubes: CubeList
2919 CubeList of the original x and y cubes for further processing.
2921 Raises
2922 ------
2923 ValueError
2924 If the cube doesn't have the right dimensions and cubes not the same
2925 size.
2926 TypeError
2927 If the cube isn't a single cube.
2929 Notes
2930 -----
2931 Scatter plots are used for determining if there is a relationship between
2932 two variables. Positive relations have a slope going from bottom left to top
2933 right; Negative relations have a slope going from top left to bottom right.
2934 """
2935 # Iterate over all cubes in cube or CubeList and plot.
2936 for cube_iter in iter_maybe(cube_x):
2937 # Check cubes are correct shape.
2938 cube_iter = check_single_cube(cube_iter)
2939 if cube_iter.ndim > 1:
2940 raise ValueError("cube_x must be 1D.")
2942 # Iterate over all cubes in cube or CubeList and plot.
2943 for cube_iter in iter_maybe(cube_y):
2944 # Check cubes are correct shape.
2945 cube_iter = check_single_cube(cube_iter)
2946 if cube_iter.ndim > 1:
2947 raise ValueError("cube_y must be 1D.")
2949 # Ensure we have a name for the plot file.
2950 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2951 title = f"{recipe_title}"
2953 if filename is None:
2954 filename = slugify(recipe_title)
2956 # Add file extension.
2957 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2959 # Do the actual plotting.
2960 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2962 # Add list of plots to plot metadata.
2963 plot_index = _append_to_plot_index([plot_filename])
2965 # Make a page to display the plots.
2966 _make_plot_html_page(plot_index)
2968 return iris.cube.CubeList([cube_x, cube_y])
2971def vector_plot(
2972 cube_u: iris.cube.Cube,
2973 cube_v: iris.cube.Cube,
2974 filename: str | None = None,
2975 sequence_coordinate: str = "time",
2976 **kwargs,
2977) -> iris.cube.CubeList:
2978 """Plot a vector plot based on the input u and v components."""
2979 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2981 # Cubes must have a matching sequence coordinate.
2982 try:
2983 # Check that the u and v cubes have the same sequence coordinate.
2984 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2984 ↛ anywhereline 2984 didn't jump anywhere: it always raised an exception.
2985 raise ValueError("Coordinates do not match.")
2986 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2987 raise ValueError(
2988 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2989 ) from err
2991 # Create a plot for each value of the sequence coordinate.
2992 plot_index = []
2993 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2994 for cube_u_slice, cube_v_slice in zip(
2995 cube_u.slices_over(sequence_coordinate),
2996 cube_v.slices_over(sequence_coordinate),
2997 strict=True,
2998 ):
2999 # Format the coordinate value in a unit appropriate way.
3000 seq_coord = cube_u_slice.coord(sequence_coordinate)
3001 plot_title, plot_filename = _set_title_and_filename(
3002 seq_coord, nplot, recipe_title, filename
3003 )
3005 # Do the actual plotting.
3006 _plot_and_save_vector_plot(
3007 cube_u_slice,
3008 cube_v_slice,
3009 filename=plot_filename,
3010 title=plot_title,
3011 method="pcolormesh",
3012 )
3013 plot_index.append(plot_filename)
3015 # Add list of plots to plot metadata.
3016 complete_plot_index = _append_to_plot_index(plot_index)
3018 # Make a page to display the plots.
3019 _make_plot_html_page(complete_plot_index)
3021 return iris.cube.CubeList([cube_u, cube_v])
3024def plot_histogram_series(
3025 cubes: iris.cube.Cube | iris.cube.CubeList,
3026 filename: str | None = None,
3027 sequence_coordinate: str = "time",
3028 stamp_coordinate: str = "realization",
3029 single_plot: bool = False,
3030 **kwargs,
3031) -> iris.cube.Cube | iris.cube.CubeList:
3032 """Plot a histogram plot for each vertical level provided.
3034 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3035 is present then a sequence of plots will be produced using the time slider
3036 functionality to scroll through histograms against time. If a
3037 stamp_coordinate is present then postage stamp plots will be produced. If
3038 stamp_coordinate and single_plot is True, all postage stamp plots will be
3039 plotted in a single plot instead of separate postage stamp plots.
3041 Parameters
3042 ----------
3043 cubes: Cube | iris.cube.CubeList
3044 Iris cube or CubeList of the data to plot. It should have a single dimension other
3045 than the stamp coordinate.
3046 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3047 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3048 filename: str, optional
3049 Name of the plot to write, used as a prefix for plot sequences. Defaults
3050 to the recipe name.
3051 sequence_coordinate: str, optional
3052 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3053 This coordinate must exist in the cube and will be used for the time
3054 slider.
3055 stamp_coordinate: str, optional
3056 Coordinate about which to plot postage stamp plots. Defaults to
3057 ``"realization"``.
3058 single_plot: bool, optional
3059 If True, all postage stamp plots will be plotted in a single plot. If
3060 False, each postage stamp plot will be plotted separately. Is only valid
3061 if stamp_coordinate exists and has more than a single point.
3063 Returns
3064 -------
3065 iris.cube.Cube | iris.cube.CubeList
3066 The original Cube or CubeList (so further operations can be applied).
3067 Plotted data.
3069 Raises
3070 ------
3071 ValueError
3072 If the cube doesn't have the right dimensions.
3073 TypeError
3074 If the cube isn't a Cube or CubeList.
3075 """
3076 recipe_title = get_recipe_metadata().get("title", "Histogram")
3078 cubes = iter_maybe(cubes)
3080 # Internal plotting function.
3081 plotting_func = _plot_and_save_histogram_series
3083 num_models = get_num_models(cubes)
3085 validate_cube_shape(cubes, num_models)
3087 # If several histograms are plotted, check sequence_coordinate
3088 check_sequence_coordinate(cubes, sequence_coordinate)
3090 # Get axis minimum and maximum from levels information.
3091 # If no levels set, derive minima and maxima from data in CubeList.
3092 vmin, vmax = _set_axis_range(cubes)
3094 # Make postage stamp plots if stamp_coordinate exists and has more than a
3095 # single point. If single_plot is True:
3096 # -- all postage stamp plots will be plotted in a single plot instead of
3097 # separate postage stamp plots.
3098 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3099 # produced per single model only
3100 if num_models == 1:
3101 if ( 3101 ↛ 3105line 3101 didn't jump to line 3105 because the condition on line 3101 was never true
3102 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3103 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3104 ):
3105 if single_plot:
3106 plotting_func = (
3107 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3108 )
3109 else:
3110 plotting_func = _plot_and_save_postage_stamp_histogram_series
3111 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3112 else:
3113 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3115 plot_index = []
3116 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3117 # Create a plot for each value of the sequence coordinate. Allowing for
3118 # multiple cubes in a CubeList to be plotted in the same plot for similar
3119 # sequence values. Passing a CubeList into the internal plotting function
3120 # for similar values of the sequence coordinate. cube_slice can be an
3121 # iris.cube.Cube or an iris.cube.CubeList.
3122 for cube_slice in cube_iterables:
3123 single_cube = cube_slice
3124 if isinstance(cube_slice, iris.cube.CubeList):
3125 single_cube = cube_slice[0]
3127 # Ensure valid stamp coordinate in cube dimensions
3128 if stamp_coordinate == "realization": 3128 ↛ 3131line 3128 didn't jump to line 3131 because the condition on line 3128 was always true
3129 stamp_coordinate = check_stamp_coordinate(single_cube)
3130 # Set plot titles and filename, based on sequence coordinate
3131 seq_coord = single_cube.coord(sequence_coordinate)
3132 # Use time coordinate in title and filename if single histogram output.
3133 if sequence_coordinate == "realization" and nplot == 1: 3133 ↛ 3134line 3133 didn't jump to line 3134 because the condition on line 3133 was never true
3134 seq_coord = single_cube.coord("time")
3135 # Use station name in title and filename if model vs obs comparison
3136 if sequence_coordinate == "station": 3136 ↛ 3137line 3136 didn't jump to line 3137 because the condition on line 3136 was never true
3137 seq_coord = single_cube.coord("Station_Name")
3139 plot_title, plot_filename = _set_title_and_filename(
3140 seq_coord, nplot, recipe_title, filename
3141 )
3143 # Do the actual plotting.
3144 plotting_func(
3145 cube_slice,
3146 filename=plot_filename,
3147 stamp_coordinate=stamp_coordinate,
3148 title=plot_title,
3149 vmin=vmin,
3150 vmax=vmax,
3151 )
3152 plot_index.append(plot_filename)
3154 # Add list of plots to plot metadata.
3155 complete_plot_index = _append_to_plot_index(plot_index)
3157 # Make a page to display the plots.
3158 _make_plot_html_page(complete_plot_index)
3160 return cubes
3163def plot_scatter_series(
3164 cubes: iris.cube.Cube | iris.cube.CubeList,
3165 filename: str | None = None,
3166 sequence_coordinate: str = "time",
3167 stamp_coordinate: str = "realization",
3168 hexbin: bool = False,
3169 **kwargs,
3170) -> iris.cube.Cube | iris.cube.CubeList:
3171 """Plot a scatter plot for each sequence coordinate provided.
3173 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3174 is present then a sequence of plots will be produced using the time slider
3175 functionality to scroll through scatter against time. If a
3176 stamp_coordinate is present then postage stamp plots will be produced. If
3177 stamp_coordinate and single_plot is True, all postage stamp plots will be
3178 plotted in a single plot instead of separate postage stamp plots.
3180 Parameters
3181 ----------
3182 cubes: Cube | iris.cube.CubeList
3183 Iris cube or CubeList of the data to plot. It should have a single dimension other
3184 than the stamp coordinate.
3185 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3186 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3187 filename: str, optional
3188 Name of the plot to write, used as a prefix for plot sequences. Defaults
3189 to the recipe name.
3190 sequence_coordinate: str, optional
3191 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3192 This coordinate must exist in the cube and will be used for the time
3193 slider.
3194 stamp_coordinate: str, optional
3195 Coordinate about which to plot postage stamp plots. Defaults to
3196 ``"realization"``.
3197 hexbin: bool, optional
3198 If True, generate hexbin comparison plot.
3199 If False, generate point-by-point scatter plot.
3201 Returns
3202 -------
3203 iris.cube.Cube | iris.cube.CubeList
3204 The original Cube or CubeList (so further operations can be applied).
3205 Plotted data.
3207 Raises
3208 ------
3209 ValueError
3210 If the cube doesn't have the right dimensions.
3211 TypeError
3212 If the cube isn't a Cube or CubeList.
3213 """
3214 recipe_title = get_recipe_metadata().get("title", "Scatter")
3216 cubes = iter_maybe(cubes)
3218 # Internal plotting function.
3219 plotting_func = _plot_and_save_scatter_series
3221 num_models = get_num_models(cubes)
3223 validate_cube_shape(cubes, num_models)
3225 check_sequence_coordinate(cubes, sequence_coordinate)
3227 vmin, vmax = _set_axis_range(cubes)
3229 # Require >1 models to compare on scatter plot
3230 if num_models > 1:
3231 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3232 else:
3233 raise ValueError(
3234 "Scatter plot series requires multiple number of models in input data."
3235 )
3237 plot_index = []
3238 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3239 # Create a plot for each value of the sequence coordinate. Allowing for
3240 # multiple cubes in a CubeList to be plotted in the same plot for similar
3241 # sequence values. Passing a CubeList into the internal plotting function
3242 # for similar values of the sequence coordinate. cube_slice can be an
3243 # iris.cube.Cube or an iris.cube.CubeList.
3244 for cube_slice in cube_iterables:
3245 single_cube = cube_slice
3246 if isinstance(cube_slice, iris.cube.CubeList): 3246 ↛ 3250line 3246 didn't jump to line 3250 because the condition on line 3246 was always true
3247 single_cube = cube_slice[0]
3249 # Ensure valid stamp coordinate in cube dimensions
3250 if stamp_coordinate == "realization": 3250 ↛ 3253line 3250 didn't jump to line 3253 because the condition on line 3250 was always true
3251 stamp_coordinate = check_stamp_coordinate(single_cube)
3252 # Set plot titles and filename, based on sequence coordinate
3253 seq_coord = single_cube.coord(sequence_coordinate)
3254 # Use time coordinate in title and filename if single histogram output.
3255 if sequence_coordinate == "realization" and nplot == 1:
3256 seq_coord = single_cube.coord("time")
3257 # Use station name in title and filename if model vs obs comparison
3258 if sequence_coordinate == "station":
3259 seq_coord = single_cube.coord("Station_Name")
3261 plot_title, plot_filename = _set_title_and_filename(
3262 seq_coord, nplot, recipe_title, filename
3263 )
3265 # Do the actual plotting.
3266 plotting_func(
3267 cube_slice,
3268 filename=plot_filename,
3269 stamp_coordinate=stamp_coordinate,
3270 title=plot_title,
3271 vmin=vmin,
3272 vmax=vmax,
3273 hexbin=hexbin,
3274 )
3275 plot_index.append(plot_filename)
3277 # Add list of plots to plot metadata.
3278 complete_plot_index = _append_to_plot_index(plot_index)
3280 # Make a page to display the plots.
3281 _make_plot_html_page(complete_plot_index)
3283 return cubes
3286def _plot_and_save_postage_stamp_power_spectrum_series(
3287 cubes: iris.cube.Cube,
3288 coords: list[iris.coords.Coord],
3289 stamp_coordinate: str,
3290 filename: str,
3291 title: str,
3292 series_coordinate: str | None = None,
3293 **kwargs,
3294):
3295 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3297 Parameters
3298 ----------
3299 cubes: Cube or CubeList
3300 Cube or Cubelist of the power spectrum data.
3301 coords: list[Coord]
3302 Coordinates to plot on the x-axis, one per cube.
3303 stamp_coordinate: str
3304 Coordinate that becomes different plots.
3305 filename: str
3306 Filename of the plot to write.
3307 title: str
3308 Plot title.
3309 series_coordinate: str, optional
3310 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3312 """
3313 # Use the smallest square grid that will fit the members.
3314 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3316 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3317 model_colors_map = get_model_colors_map(cubes)
3318 # ax = plt.gca()
3319 # Make a subplot for each member.
3320 for member, subplot in zip(
3321 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3322 ):
3323 ax = plt.subplot(grid_size, grid_size, subplot)
3325 # Store min/max ranges.
3326 y_levels = []
3328 line_marker = None
3329 line_width = 1
3331 for cube in iter_maybe(member):
3332 xcoord = _select_series_coord(cube, series_coordinate)
3333 xname = xcoord.points
3335 yfield = cube.data # power spectrum
3336 label = None
3337 color = "black"
3338 if model_colors_map: 3338 ↛ 3339line 3338 didn't jump to line 3339 because the condition on line 3338 was never true
3339 label = cube.attributes.get("model_name")
3340 color = model_colors_map.get(label)
3342 if member.coord(stamp_coordinate).points == [0]:
3343 ax.plot(
3344 xname,
3345 yfield,
3346 color=color,
3347 marker=line_marker,
3348 ls="-",
3349 lw=line_width,
3350 label=f"{label} (control)"
3351 if len(cube.coord(stamp_coordinate).points) > 1
3352 else label,
3353 )
3354 # Label with member if part of an ensemble and not the control.
3355 else:
3356 ax.plot(
3357 xname,
3358 yfield,
3359 color=color,
3360 ls="-",
3361 lw=1.5,
3362 alpha=0.75,
3363 label=f"{label} (member)",
3364 )
3366 # Calculate the global min/max if multiple cubes are given.
3367 _, levels, _ = colorbar_map_levels(cube, axis="y")
3368 if levels is not None: 3368 ↛ 3369line 3368 didn't jump to line 3369 because the condition on line 3368 was never true
3369 y_levels.append(min(levels))
3370 y_levels.append(max(levels))
3372 # Add some labels and tweak the style.
3373 title = f"{title}"
3374 ax.set_title(title, fontsize=16)
3376 # Set appropriate x-axis label based on coordinate
3377 if series_coordinate == "wavelength" or ( 3377 ↛ 3380line 3377 didn't jump to line 3380 because the condition on line 3377 was never true
3378 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3379 ):
3380 ax.set_xlabel("Wavelength (km)", fontsize=14)
3381 elif series_coordinate == "physical_wavenumber" or ( 3381 ↛ 3386line 3381 didn't jump to line 3386 because the condition on line 3381 was always true
3382 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3383 ):
3384 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3385 else: # frequency or check units
3386 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3387 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3388 else:
3389 ax.set_xlabel("Wavenumber", fontsize=14)
3391 ax.set_ylabel("Power Spectral Density", fontsize=14)
3392 ax.tick_params(axis="both", labelsize=12)
3394 # Set log-log scale
3395 ax.set_xscale("log")
3396 ax.set_yscale("log")
3398 # Add gridlines
3399 ax.grid(linestyle="--", color="grey", linewidth=1)
3400 # Ientify unique labels for legend
3401 handles = list(
3402 {
3403 label: handle
3404 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3405 }.values()
3406 )
3407 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3409 ax = plt.gca()
3410 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3412 # Save plot.
3413 _save_close_figure(fig, "histogram postage stamp", filename)
3416def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3417 cubes: iris.cube.Cube,
3418 coords: list[iris.coords.Coord],
3419 stamp_coordinate: str,
3420 filename: str,
3421 title: str,
3422 series_coordinate: str | None = None,
3423 **kwargs,
3424):
3425 """Plot and save power spectra for ensemble members in single plot.
3427 Parameters
3428 ----------
3429 cubes: Cube or CubeList
3430 Cube or Cubelist of the power spectrum data.
3431 coords: list[Coord]
3432 Coordinates to plot on the x-axis, one per cube.
3433 stamp_coordinate: str
3434 Coordinate that becomes different plots.
3435 filename: str
3436 Filename of the plot to write.
3437 title: str
3438 Plot title.
3439 series_coordinate: str, optional
3440 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3442 """
3443 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3444 model_colors_map = get_model_colors_map(cubes)
3446 line_marker = None
3447 line_width = 1
3449 # Compute ensemble statistics to show spread
3450 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3451 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3452 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3454 xcoord_global = mean_cube.coord(series_coordinate)
3455 x_global = xcoord_global.points
3457 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3458 xcoord = _select_series_coord(member, series_coordinate)
3459 xname = xcoord.points
3461 yfield = member.data # power spectrum
3462 color = "black"
3463 if model_colors_map: 3463 ↛ 3467line 3463 didn't jump to line 3467 because the condition on line 3463 was always true
3464 label = member.attributes.get("model_name") if i == 0 else None
3465 color = model_colors_map.get(label)
3467 if member.coord(stamp_coordinate).points == [0]:
3468 ax.plot(
3469 xname,
3470 yfield,
3471 color=color,
3472 marker=line_marker,
3473 ls="-",
3474 lw=line_width,
3475 label=f"{label} (control)"
3476 if len(member.coord(stamp_coordinate).points) > 1
3477 else label,
3478 )
3479 # Label with member number if part of an ensemble and not the control.
3480 else:
3481 ax.plot(
3482 xname,
3483 yfield,
3484 color=color,
3485 ls="-",
3486 lw=1.5,
3487 alpha=0.75,
3488 label=label,
3489 )
3491 # Set appropriate x-axis label based on coordinate
3492 if series_coordinate == "wavelength" or ( 3492 ↛ 3495line 3492 didn't jump to line 3495 because the condition on line 3492 was never true
3493 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3494 ):
3495 ax.set_xlabel("Wavelength (km)", fontsize=14)
3496 elif series_coordinate == "physical_wavenumber" or ( 3496 ↛ 3501line 3496 didn't jump to line 3501 because the condition on line 3496 was always true
3497 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3498 ):
3499 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3500 else: # frequency or check units
3501 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3502 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3503 else:
3504 ax.set_xlabel("Wavenumber", fontsize=14)
3506 # Add ensemble spread shading
3507 ax.fill_between(
3508 x_global,
3509 min_cube.data,
3510 max_cube.data,
3511 color="grey",
3512 alpha=0.3,
3513 label="Ensemble spread",
3514 )
3516 # Add ensemble mean line
3517 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3519 ax.set_ylabel("Power Spectral Density", fontsize=14)
3520 ax.tick_params(axis="both", labelsize=12)
3522 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3523 # Set log-log scale
3524 ax.set_xscale("log")
3525 ax.set_yscale("log")
3527 # Add gridlines
3528 ax.grid(linestyle="--", color="grey", linewidth=1)
3529 # Identify unique labels for legend
3530 handles = list(
3531 {
3532 label: handle
3533 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3534 }.values()
3535 )
3536 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3538 # Figure title.
3539 ax.set_title(title, fontsize=16)
3541 # Save plot.
3542 _save_close_figure(fig, "power spectra postage stamp", filename)