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