Coverage for src/CSET/operators/plot.py: 83%
1108 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-05 16:16 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24import sys
25from typing import Literal
27import cartopy.crs as ccrs
28import cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.pyplot as plt
36import numpy as np
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._colormaps import (
50 colorbar_map_levels,
51 get_model_colors_map,
52)
53from CSET.operators._utils import (
54 check_sequence_coordinate,
55 check_single_cube,
56 check_stamp_coordinate,
57 fully_equalise_attributes,
58 get_cube_yxcoordname,
59 get_num_models,
60 is_transect,
61 slice_over_maybe,
62 validate_cube_shape,
63 validate_cubes_coords,
64)
65from CSET.operators.collapse import collapse
66from CSET.operators.misc import _extract_common_time_points
67from CSET.operators.regrid import regrid_onto_cube
69logger = logging.getLogger(__name__)
71# Use a non-interactive plotting backend.
72mpl.use("agg")
75############################
76# Private helper functions #
77############################
80def in_sphinx_gallery():
81 """Test if running plot code in sphinx-gallery context."""
82 return "sphinx_gallery" in sys.modules
85def _append_to_plot_index(plot_index: list) -> list:
86 """Add plots into the plot index, returning the complete plot index."""
87 with open("meta.json", "r+t", encoding="UTF-8") as fp:
88 fcntl.flock(fp, fcntl.LOCK_EX)
89 fp.seek(0)
90 meta = json.load(fp)
91 complete_plot_index = meta.get("plots", [])
92 complete_plot_index = complete_plot_index + plot_index
93 meta["plots"] = complete_plot_index
94 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
95 os.getenv("DO_CASE_AGGREGATION")
96 ):
97 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
98 fp.seek(0)
99 fp.truncate()
100 json.dump(meta, fp, indent=2)
101 return complete_plot_index
104def _make_plot_html_page(plots: list):
105 """Create a HTML page to display a plot image."""
106 # Debug check that plots actually contains some strings.
107 assert isinstance(plots[0], str)
109 # Load HTML template file.
110 operator_files = importlib.resources.files()
111 template_file = operator_files.joinpath("_plot_page_template.html")
113 # Get some metadata.
114 meta = get_recipe_metadata()
115 title = meta.get("title", "Untitled")
116 description = MarkdownIt().render(meta.get("description", "*No description.*"))
118 # Prepare template variables.
119 variables = {
120 "title": title,
121 "description": description,
122 "initial_plot": plots[0],
123 "plots": plots,
124 "title_slug": slugify(title),
125 }
127 # Render template.
128 html = render_file(template_file, **variables)
130 # Save completed HTML.
131 with open("index.html", "wt", encoding="UTF-8") as fp:
132 fp.write(html)
135def _save_close_figure(figure, plot_type: str, filename: str):
136 """Save generated plot figure file and close figure.
138 If running documentation gallery generation, avoid saving to file.
140 Parameters
141 ----------
142 figure:
143 Matplotlib Figure object holding all plot elements.
144 plot_type: str
145 String identifier for plot type for logging information.
146 filename: str
147 Filename for saved figure.
148 """
149 if not in_sphinx_gallery():
150 figure.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
151 logger.info("Saved %s plot to %s", plot_type, filename)
152 plt.close(figure)
155def _setup_spatial_map(
156 cube: iris.cube.Cube,
157 figure,
158 cmap,
159 grid_size: tuple[int, int] | None = None,
160 subplot: int | None = None,
161):
162 """Define map projections, extent and add coastlines and borderlines for spatial plots.
164 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
165 is specified, and map extent defined based on the input data.
167 Parameters
168 ----------
169 cube: Cube
170 2 dimensional (lat and lon) Cube of the data to plot.
171 figure:
172 Matplotlib Figure object holding all plot elements.
173 cmap:
174 Matplotlib colormap.
175 grid_size: (int, int), optional
176 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
177 subplot: int, optional
178 Subplot index if multiple spatial subplots in figure.
180 Returns
181 -------
182 axes:
183 Matplotlib GeoAxes definition.
184 """
185 # Identify min/max plot bounds.
186 try:
187 lat_axis, lon_axis = get_cube_yxcoordname(cube)
188 xmin = np.nanmin(cube.coord(lon_axis).points)
189 xmax = np.nanmax(cube.coord(lon_axis).points)
190 ymin = np.nanmin(cube.coord(lat_axis).points)
191 ymax = np.nanmax(cube.coord(lat_axis).points)
193 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
194 if np.abs(xmax - xmin) > 180.0:
195 xmin = xmin - 180.0
196 xmax = xmax - 180.0
197 logger.debug("Adjusting plot bounds to fit global extent.")
199 # Consider map projection orientation.
200 # Adapting orientation enables plotting across international dateline.
201 # Users can adapt the default central_longitude if alternative projections views.
202 if xmax > 180.0 or xmin < -180.0:
203 central_longitude = 180.0
204 else:
205 central_longitude = 0.0
207 # Define spatial map projection.
208 coord_system = cube.coord(lat_axis).coord_system
209 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
210 # Define rotated pole map projection for rotated pole inputs.
211 projection = ccrs.RotatedPole(
212 pole_longitude=coord_system.grid_north_pole_longitude,
213 pole_latitude=coord_system.grid_north_pole_latitude,
214 central_rotated_longitude=central_longitude,
215 )
216 crs = projection
217 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 217 ↛ 219line 217 didn't jump to line 219 because the condition on line 217 was never true
218 # Define Transverse Mercator projection for TM inputs.
219 projection = ccrs.TransverseMercator(
220 central_longitude=coord_system.longitude_of_central_meridian,
221 central_latitude=coord_system.latitude_of_projection_origin,
222 false_easting=coord_system.false_easting,
223 false_northing=coord_system.false_northing,
224 scale_factor=coord_system.scale_factor_at_central_meridian,
225 )
226 crs = projection
227 else:
228 # Assume polar projection for regional grids encompassing N. Pole
229 if ymin > 20.0 and ymax > 80.0:
230 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
231 elif ymin < -80.0 and ymax < -20.0:
232 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
233 # Define regular map projection for non-rotated pole inputs.
234 # Alternatives might include e.g. for global model outputs:
235 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
236 # projection = ccrs.NearsidePerspective(
237 # central_longitude=180.0,
238 # central_latitude=0,
239 # satellite_height=35785831,
240 # )
241 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
242 else:
243 projection = ccrs.PlateCarree(central_longitude=central_longitude)
244 crs = ccrs.PlateCarree()
246 # Define axes for plot (or subplot) with required map projection.
247 if subplot is not None:
248 axes = figure.add_subplot(
249 grid_size[0], grid_size[1], subplot, projection=projection
250 )
251 else:
252 axes = figure.add_subplot(projection=projection)
254 # Add coastlines and borderlines if cube contains x and y map coordinates.
255 # Avoid adding lines for 2D masked data or specific fixed ancillary spatial plots.
256 if (cube.ndim > 1 and iris.util.is_masked(cube.data)) or any(
257 name in cube.name() for name in ["land_", "orography", "altitude"]
258 ):
259 pass
260 else:
261 if cmap.name in ["viridis", "Greys"]:
262 coastcol = "magenta"
263 else:
264 coastcol = "black"
265 logger.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
266 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
267 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
269 # Add gridlines.
270 gl = axes.gridlines(
271 alpha=0.3,
272 draw_labels=True,
273 dms=False,
274 x_inline=False,
275 y_inline=False,
276 )
277 gl.top_labels = False
278 gl.right_labels = False
279 if subplot:
280 gl.bottom_labels = False
281 gl.left_labels = False
282 if subplot % grid_size[1] == 1:
283 gl.left_labels = True
284 if subplot > ((grid_size[0] - 1) * grid_size[1]): 284 ↛ 289line 284 didn't jump to line 289 because the condition on line 284 was always true
285 gl.bottom_labels = True
287 # If is lat/lon spatial map, fix extent to keep plot tight.
288 # Specifying crs within set_extent helps ensure only data region is shown.
289 if isinstance(
290 coord_system, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
291 ):
292 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
294 except ValueError:
295 # Skip if not both x and y map coordinates.
296 axes = figure.gca()
298 return axes
301def _get_plot_resolution() -> int:
302 """Get resolution of rasterised plots in pixels per inch."""
303 return get_recipe_metadata().get("plot_resolution", 100)
306def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
307 """Return title and filename based on start and end points or bounds."""
308 if use_bounds and seq_coord.has_bounds():
309 vals = seq_coord.bounds.flatten()
310 else:
311 vals = seq_coord.points
312 start = seq_coord.units.title(vals[0])
313 end = seq_coord.units.title(vals[-1])
315 if start == end:
316 sequence_title = f"\n [{start}]"
317 sequence_fname = f"_{filename_slugify(start)}"
318 else:
319 sequence_title = f"\n [{start} to {end}]"
320 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
322 # Do not include time if coord set to zero.
323 if (
324 seq_coord.units == "hours since 0001-01-01 00:00:00"
325 and vals[0] == 0
326 and vals[-1] == 0
327 ):
328 sequence_title = ""
329 sequence_fname = ""
331 return sequence_title, sequence_fname
334def _set_title_and_filename(
335 seq_coord: iris.coords.Coord,
336 nplot: int,
337 recipe_title: str,
338 filename: str,
339):
340 """Set plot title and filename based on cube coordinate.
342 Parameters
343 ----------
344 sequence_coordinate: iris.coords.Coord
345 Coordinate about which to make a plot sequence.
346 nplot: int
347 Number of output plots to generate - controls title/naming.
348 recipe_title: str
349 Default plot title, potentially to update.
350 filename: str
351 Input plot filename, potentially to update.
353 Returns
354 -------
355 plot_title: str
356 Output formatted plot title string, based on plotted data.
357 plot_filename: str
358 Output formatted plot filename string.
359 """
360 ndim = seq_coord.ndim
361 npoints = np.size(seq_coord.points)
362 sequence_title = ""
363 sequence_fname = ""
365 # Case 1: Multiple dimension sequence input - list number of aggregated cases
366 # (e.g. aggregation histogram plots)
367 if ndim > 1:
368 ncase = np.shape(seq_coord)[0]
369 sequence_title = f"\n [{ncase} cases]"
370 sequence_fname = f"_{ncase}cases"
372 # Case 2: Single dimension input
373 else:
374 # Single sequence point
375 if npoints == 1:
376 if nplot > 1:
377 # Default labels for sequence inputs
378 sequence_value = seq_coord.units.title(seq_coord.points[0])
379 sequence_value = sequence_value.replace(" unknown", "")
380 sequence_title = f"\n [{sequence_value}]"
381 sequence_fname = f"_{filename_slugify(sequence_value)}"
382 else:
383 # Aggregated attribute available where input collapsed over aggregation
384 try:
385 ncase = seq_coord.attributes["number_reference_times"]
386 sequence_title = f"\n [{ncase} cases]"
387 sequence_fname = f"_{ncase}cases"
388 except KeyError:
389 sequence_title, sequence_fname = _get_start_end_strings(
390 seq_coord, use_bounds=seq_coord.has_bounds()
391 )
392 # Multiple sequence (e.g. time) points
393 else:
394 sequence_title, sequence_fname = _get_start_end_strings(
395 seq_coord, use_bounds=False
396 )
398 # Set plot title and filename
399 plot_title = f"{recipe_title}{sequence_title}"
401 # Set plot filename, defaulting to user input if provided.
402 if filename is None:
403 filename = slugify(recipe_title)
404 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
405 else:
406 if nplot > 1:
407 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
408 else:
409 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
411 return plot_title, plot_filename
414def _select_series_coord(cube, series_coordinate):
415 """Determine the grid coordinates to use to calculate grid spacing."""
416 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
417 if series_coordinate in spacing_coordinates: 417 ↛ 423line 417 didn't jump to line 423 because the condition on line 417 was always true
418 # Try the requested coordinate first then the fallbacks in order.
419 fallbacks = [series_coordinate] + [
420 c for c in spacing_coordinates if c != series_coordinate
421 ]
422 else:
423 fallbacks = {series_coordinate}
425 # Try each possible coordinate.
426 for coord in fallbacks:
427 try:
428 return cube.coord(coord)
429 except iris.exceptions.CoordinateNotFoundError:
430 logger.debug("Coordinate %s not found.", coord)
432 # If we get here, none of the fallback options were found.
433 raise iris.exceptions.CoordinateNotFoundError(
434 f"No valid coordinate found for '{series_coordinate}' "
435 f"or fallback options {fallbacks}"
436 )
439def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
440 """Control postage stamp plot output titles based on stamp coordinate."""
441 if stamp_coord.name() == "realization":
442 mtitle = "Member"
443 else:
444 mtitle = stamp_coord.name().capitalize()
446 if stamp_coord.name() == "time":
447 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
448 else:
449 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
451 return mtitle
454def _set_axis_range(cubes):
455 """Get minimum and maximum from levels information."""
456 levels = None
457 for cube in cubes: 457 ↛ 473line 457 didn't jump to line 473 because the loop on line 457 didn't complete
458 # First check if user-specified "auto" range variable.
459 # This maintains the value of levels as None, so proceed.
460 _, levels, _ = colorbar_map_levels(cube, axis="y")
461 if levels is None:
462 break
463 # If levels is changed, recheck to use the vmin,vmax or
464 # levels-based ranges for histogram plots.
465 _, levels, _ = colorbar_map_levels(cube)
466 logger.debug("levels: %s", levels)
467 if levels is not None: 467 ↛ 457line 467 didn't jump to line 457 because the condition on line 467 was always true
468 vmin = min(levels)
469 vmax = max(levels)
470 logger.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
471 break
473 if levels is None:
474 vmin = min(cb.data.min() for cb in cubes)
475 vmax = max(cb.data.max() for cb in cubes)
477 return vmin, vmax
480def _find_matched_slices(cubes, sequence_coordinate):
481 """Identify matched cubes in CubeList by sequence_coordinate values.
483 Ensures common points are compared for multiple cube inputs.
484 """
485 all_points = sorted(
486 set(
487 itertools.chain.from_iterable(
488 cb.coord(sequence_coordinate).points for cb in cubes
489 )
490 )
491 )
492 all_slices = list(
493 itertools.chain.from_iterable(
494 cb.slices_over(sequence_coordinate) for cb in cubes
495 )
496 )
497 # Matched slices (matched by seq coord point; it may happen that
498 # evaluated models do not cover the same seq coord range, hence matching
499 # necessary)
500 cube_iterables = [
501 iris.cube.CubeList(
502 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
503 )
504 for point in all_points
505 ]
507 return cube_iterables
510def _plot_and_save_spatial_plot(
511 cube: iris.cube.Cube,
512 filename: str,
513 title: str,
514 method: Literal["contourf", "pcolormesh", "scatter"],
515 overlay_cube: iris.cube.Cube | None = None,
516 contour_cube: iris.cube.Cube | None = None,
517 point_cube: iris.cube.Cube | None = None,
518 **kwargs,
519):
520 """Plot and save a spatial plot.
522 Parameters
523 ----------
524 cube: Cube
525 2 dimensional (lat and lon) Cube of the data to plot.
526 filename: str
527 Filename of the plot to write.
528 title: str
529 Plot title.
530 method: "contourf" | "pcolormesh" | "scatter"
531 The plotting method to use
532 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
533 overlay_cube: Cube, optional
534 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
535 contour_cube: Cube, optional
536 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
537 point_cube: Cube, optional
538 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
539 """
540 # Setup plot details, size, resolution, etc.
541 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
543 # Specify the color bar
544 cmap, levels, norm = colorbar_map_levels(cube)
546 # If overplotting, set required colorbars
547 if overlay_cube:
548 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
549 if contour_cube:
550 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
552 # Setup plot map projection, extent and coastlines and borderlines.
553 axes = _setup_spatial_map(cube, fig, cmap)
555 # Set colorscale bounds
556 try:
557 vmin = min(levels)
558 vmax = max(levels)
559 except TypeError:
560 vmin, vmax = None, None
561 # Ensure to use norm and not vmin/vmax if levels are defined.
562 if norm is not None:
563 vmin = None
564 vmax = None
565 logger.debug("Plotting using defined levels.")
567 # Plot the field.
568 if method == "contourf":
569 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
570 elif method == "pcolormesh":
571 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
572 elif method == "scatter":
573 # Scatter plot of the field. The marker size is chosen to give
574 # symbols that decrease in size as the number of data points
575 # increases, although the fraction of the figure covered by
576 # symbols increases roughly as N^(1/2), disregarding overlaps,
577 # and has been selected for the default figure size of (10, 10).
578 # Should this be changed, the marker size should be adjusted in
579 # proportion to the area of the figure.
580 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
581 lat_axis, lon_axis = get_cube_yxcoordname(cube)
582 plot = iplt.scatter(
583 cube.coord(lon_axis),
584 cube.coord(lat_axis),
585 c=cube.data[:],
586 s=mrk_size,
587 cmap=cmap,
588 edgecolors="k",
589 norm=norm,
590 vmin=vmin,
591 vmax=vmax,
592 )
593 else:
594 raise ValueError(f"Unknown plotting method: {method}")
596 # Overplot overlay field, if required
597 if overlay_cube:
598 try:
599 over_vmin = min(over_levels)
600 over_vmax = max(over_levels)
601 except TypeError:
602 over_vmin, over_vmax = None, None
603 if over_norm is not None: 603 ↛ 604line 603 didn't jump to line 604 because the condition on line 603 was never true
604 over_vmin = None
605 over_vmax = None
606 overlay = iplt.pcolormesh(
607 overlay_cube,
608 cmap=over_cmap,
609 norm=over_norm,
610 alpha=0.8,
611 vmin=over_vmin,
612 vmax=over_vmax,
613 )
614 # Overplot contour field, if required, with contour labelling.
615 if contour_cube:
616 contour = iplt.contour(
617 contour_cube,
618 colors="darkgray",
619 levels=cntr_levels,
620 norm=cntr_norm,
621 alpha=0.5,
622 linestyles="--",
623 linewidths=1,
624 )
625 plt.clabel(contour)
626 # Overplot valid elements of point-based field, if required.
627 # Check for non-masked points only to avoid plotting missing data.
628 if point_cube:
629 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
630 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
631 lon_coord = point_cube.coord(lon_axis)
632 lat_coord = point_cube.coord(lat_axis)
633 valid = ~point_cube.data.mask
634 valid_lon = iris.coords.AuxCoord(
635 lon_coord.points[valid],
636 standard_name=lon_coord.standard_name,
637 units=lon_coord.units,
638 coord_system=lon_coord.coord_system,
639 )
640 valid_lat = iris.coords.AuxCoord(
641 lat_coord.points[valid],
642 standard_name=lat_coord.standard_name,
643 units=lat_coord.units,
644 coord_system=lat_coord.coord_system,
645 )
646 iplt.scatter(
647 valid_lon,
648 valid_lat,
649 c=point_cube.data[valid],
650 s=mrk_size,
651 cmap=cmap,
652 edgecolors="k",
653 norm=norm,
654 vmin=vmin,
655 vmax=vmax,
656 )
658 # Check to see if transect, and if so, adjust y axis.
659 if is_transect(cube):
660 if "pressure" in [coord.name() for coord in cube.coords()]:
661 axes.invert_yaxis()
662 axes.set_yscale("log")
663 axes.set_ylim(1100, 100)
664 # If both model_level_number and level_height exists, iplt can construct
665 # plot as a function of height above orography (NOT sea level).
666 elif {"model_level_number", "level_height"}.issubset( 666 ↛ 671line 666 didn't jump to line 671 because the condition on line 666 was always true
667 {coord.name() for coord in cube.coords()}
668 ):
669 axes.set_yscale("log")
671 axes.set_title(
672 f"{title}\n"
673 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
674 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
675 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
676 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
677 fontsize=16,
678 )
680 # Inset code
681 axins = inset_axes(
682 axes,
683 width="20%",
684 height="20%",
685 loc="upper right",
686 axes_class=GeoAxes,
687 axes_kwargs={"map_projection": ccrs.PlateCarree()},
688 )
690 # Slightly transparent to reduce plot blocking.
691 axins.patch.set_alpha(0.4)
693 axins.coastlines(resolution="50m")
694 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
696 SLat, SLon, ELat, ELon = (
697 float(coord) for coord in cube.attributes["transect_coords"].split("_")
698 )
700 # Draw line between them
701 axins.plot(
702 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
703 )
705 # Plot points (note: lon, lat order for Cartopy)
706 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
707 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
709 lon_min, lon_max = sorted([SLon, ELon])
710 lat_min, lat_max = sorted([SLat, ELat])
712 # Midpoints
713 lon_mid = (lon_min + lon_max) / 2
714 lat_mid = (lat_min + lat_max) / 2
716 # Maximum half-range
717 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
718 if half_range == 0: # points identical → provide small default 718 ↛ 722line 718 didn't jump to line 722 because the condition on line 718 was always true
719 half_range = 1
721 # Set square extent
722 axins.set_extent(
723 [
724 lon_mid - half_range,
725 lon_mid + half_range,
726 lat_mid - half_range,
727 lat_mid + half_range,
728 ],
729 crs=ccrs.PlateCarree(),
730 )
732 # Ensure square aspect
733 axins.set_aspect("equal")
735 else:
736 # Add title.
737 axes.set_title(title, fontsize=16)
739 # Adjust padding if spatial plot or transect
740 if is_transect(cube):
741 yinfopad = -0.1
742 ycbarpad = 0.1
743 else:
744 yinfopad = 0.01
745 ycbarpad = 0.042
747 # Add watermark with min/max/mean. Currently not user togglable.
748 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
749 axes.annotate(
750 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
751 xy=(0.025, yinfopad),
752 xycoords="axes fraction",
753 xytext=(-5, 5),
754 textcoords="offset points",
755 ha="left",
756 va="bottom",
757 size=11,
758 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
759 )
761 # Add secondary colour bar for overlay_cube field if required.
762 if overlay_cube:
763 cbarB = fig.colorbar(
764 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
765 )
766 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
767 # add ticks and tick_labels for every levels if less than 20 levels exist
768 if over_levels is not None and len(over_levels) < 20: 768 ↛ 769line 768 didn't jump to line 769 because the condition on line 768 was never true
769 cbarB.set_ticks(over_levels)
770 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
771 if any(
772 var in overlay_cube.name()
773 for var in ("rainfall", "snowfall", "visibility")
774 ):
775 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
776 logger.debug("Set secondary colorbar ticks and labels.")
778 # Add main colour bar.
779 cbar = fig.colorbar(
780 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
781 )
783 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
784 # add ticks and tick_labels for every levels if less than 20 levels exist
785 if levels is not None and len(levels) < 20:
786 cbar.set_ticks(levels)
787 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
788 if any(var in cube.name() for var in ("rainfall", "snowfall", "visibility")): 788 ↛ 791line 788 didn't jump to line 791 because the condition on line 788 was always true
789 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
790 # Tick labels for rainfall rates from Nimrod radar data.
791 if "rainfall rate composite" in cube.name(): 791 ↛ 792line 791 didn't jump to line 792 because the condition on line 791 was never true
792 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
793 # Tick labels for rain accumulations from Nimrod radar data.
794 if "rain accumulation" in cube.name(): 794 ↛ 795line 794 didn't jump to line 795 because the condition on line 794 was never true
795 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
796 if "wts accumulation" in cube.name(): 796 ↛ 797line 796 didn't jump to line 797 because the condition on line 796 was never true
797 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
798 cbar.minorticks_off()
799 cbar.set_ticks(tick_levels)
800 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
801 cbar.set_label(label=f"{cube.name()}", size=14)
802 # Tick labels for model rainfall data.
803 if "surface_microphysical" in cube.name(): 803 ↛ 806line 803 didn't jump to line 806 because the condition on line 803 was always true
804 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
805 # Tick labels for Nimrod weights data.
806 logger.debug("Set colorbar ticks and labels.")
808 # Save plot.
809 _save_close_figure(fig, "spatial", filename)
812def _plot_and_save_postage_stamp_spatial_plot(
813 cube: iris.cube.Cube,
814 filename: str,
815 stamp_coordinate: str,
816 title: str,
817 method: Literal["contourf", "pcolormesh"],
818 overlay_cube: iris.cube.Cube | None = None,
819 contour_cube: iris.cube.Cube | None = None,
820 **kwargs,
821):
822 """Plot postage stamp spatial plots from an ensemble.
824 Parameters
825 ----------
826 cube: Cube
827 Iris cube of data to be plotted. It must have the stamp coordinate.
828 filename: str
829 Filename of the plot to write.
830 stamp_coordinate: str
831 Coordinate that becomes different plots.
832 method: "contourf" | "pcolormesh"
833 The plotting method to use.
834 overlay_cube: Cube, optional
835 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
836 contour_cube: Cube, optional
837 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
839 Raises
840 ------
841 ValueError
842 If the cube doesn't have the right dimensions.
843 """
844 # Use the smallest square grid that will fit the members.
845 nmember = len(cube.coord(stamp_coordinate).points)
846 grid_rows = int(math.sqrt(nmember))
847 grid_size = math.ceil(nmember / grid_rows)
849 fig = plt.figure(
850 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
851 )
853 # Specify the color bar
854 cmap, levels, norm = colorbar_map_levels(cube)
855 # If overplotting, set required colorbars
856 if overlay_cube: 856 ↛ 857line 856 didn't jump to line 857 because the condition on line 856 was never true
857 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
858 if contour_cube: 858 ↛ 859line 858 didn't jump to line 859 because the condition on line 858 was never true
859 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
861 # Make a subplot for each member.
862 for member, subplot in zip(
863 cube.slices_over(stamp_coordinate),
864 range(1, grid_size * grid_rows + 1),
865 strict=False,
866 ):
867 # Setup subplot map projection, extent and coastlines and borderlines.
868 axes = _setup_spatial_map(
869 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
870 )
871 if method == "contourf":
872 # Filled contour plot of the field.
873 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
874 elif method == "pcolormesh":
875 if levels is not None:
876 vmin = min(levels)
877 vmax = max(levels)
878 else:
879 raise TypeError("Unknown vmin and vmax range.")
880 vmin, vmax = None, None
881 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
882 # if levels are defined.
883 if norm is not None: 883 ↛ 884line 883 didn't jump to line 884 because the condition on line 883 was never true
884 vmin = None
885 vmax = None
886 # pcolormesh plot of the field.
887 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
888 else:
889 raise ValueError(f"Unknown plotting method: {method}")
891 # Overplot overlay field, if required
892 if overlay_cube: 892 ↛ 893line 892 didn't jump to line 893 because the condition on line 892 was never true
893 try:
894 over_vmin = min(over_levels)
895 over_vmax = max(over_levels)
896 except TypeError:
897 over_vmin, over_vmax = None, None
898 if over_norm is not None:
899 over_vmin = None
900 over_vmax = None
901 iplt.pcolormesh(
902 overlay_cube[member.coord(stamp_coordinate).points[0]],
903 cmap=over_cmap,
904 norm=over_norm,
905 alpha=0.6,
906 vmin=over_vmin,
907 vmax=over_vmax,
908 )
909 # Overplot contour field, if required
910 if contour_cube: 910 ↛ 911line 910 didn't jump to line 911 because the condition on line 910 was never true
911 iplt.contour(
912 contour_cube[member.coord(stamp_coordinate).points[0]],
913 colors="darkgray",
914 levels=cntr_levels,
915 norm=cntr_norm,
916 alpha=0.6,
917 linestyles="--",
918 linewidths=1,
919 )
920 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
921 axes.set_title(f"{mtitle}")
923 # Put the shared colorbar in its own axes.
924 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
925 colorbar = fig.colorbar(
926 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
927 )
928 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
930 # Overall figure title.
931 fig.suptitle(title, fontsize=16)
933 # Save plot.
934 _save_close_figure(fig, "contour postate stamp", filename)
937def _plot_and_save_line_series(
938 cubes: iris.cube.CubeList,
939 coords: list[iris.coords.Coord],
940 ensemble_coord: str,
941 filename: str,
942 title: str,
943 **kwargs,
944):
945 """Plot and save a 1D line series.
947 Parameters
948 ----------
949 cubes: Cube or CubeList
950 Cube or CubeList containing the cubes to plot on the y-axis.
951 coords: list[Coord]
952 Coordinates to plot on the x-axis, one per cube.
953 ensemble_coord: str
954 Ensemble coordinate in the cube.
955 filename: str
956 Filename of the plot to write.
957 title: str
958 Plot title.
959 """
960 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
962 title_suffix = ""
963 if cubes and cubes[0].attributes.get("cset_ensemble_mean") == "true": 963 ↛ 964line 963 didn't jump to line 964 because the condition on line 963 was never true
964 title_suffix = " (ensemble mean)"
966 model_colors_map = get_model_colors_map(cubes)
968 # Store min/max ranges.
969 y_levels = []
971 # Check match-up across sequence coords gives consistent sizes
972 validate_cubes_coords(cubes, coords)
974 for cube, coord in zip(cubes, coords, strict=True):
975 label = None
976 color = "black"
977 if model_colors_map:
978 label = cube.attributes.get("model_name")
979 color = model_colors_map.get(label)
980 if not cube.coords(ensemble_coord): 980 ↛ 983line 980 didn't jump to line 983 because the condition on line 980 was never true
981 # No ensemble coordinate (e.g. after ensemble-mean collapse): plot
982 # the cube as a single deterministic line.
983 iplt.plot(coord, cube, color=color, marker="o", ls="-", lw=3, label=label)
984 else:
985 for cube_slice in cube.slices_over(ensemble_coord):
986 # Label with (control) if part of an ensemble or not otherwise.
987 if cube_slice.coord(ensemble_coord).points == [0]:
988 iplt.plot(
989 coord,
990 cube_slice,
991 color=color,
992 marker="o",
993 ls="-",
994 lw=3,
995 label=f"{label} (control)"
996 if len(cube.coord(ensemble_coord).points) > 1
997 else label,
998 )
999 # Label with (perturbed) if part of an ensemble and not the control.
1000 else:
1001 iplt.plot(
1002 coord,
1003 cube_slice,
1004 color=color,
1005 ls="-",
1006 lw=1.5,
1007 alpha=0.75,
1008 label=f"{label} (member)",
1009 )
1011 # Calculate the global min/max if multiple cubes are given.
1012 _, levels, _ = colorbar_map_levels(cube, axis="y")
1013 if levels is not None: 1013 ↛ 1014line 1013 didn't jump to line 1014 because the condition on line 1013 was never true
1014 y_levels.append(min(levels))
1015 y_levels.append(max(levels))
1017 # Get the current axes.
1018 ax = plt.gca()
1020 # Add some labels and tweak the style.
1021 # check if cubes[0] works for single cube if not CubeList
1022 if coords[0].name() == "time":
1023 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1024 else:
1025 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1026 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1027 ax.set_title(f"{title}{title_suffix}", fontsize=16)
1029 ax.ticklabel_format(axis="y", useOffset=False)
1030 ax.tick_params(axis="x", labelrotation=15)
1031 ax.tick_params(axis="both", labelsize=12)
1033 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1034 if y_levels: 1034 ↛ 1035line 1034 didn't jump to line 1035 because the condition on line 1034 was never true
1035 ax.set_ylim(min(y_levels), max(y_levels))
1036 logger.debug("Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels))
1037 else:
1038 ax.autoscale()
1040 # Add gridlines
1041 ax.grid(linestyle="--", color="grey", linewidth=1)
1042 # Add zero line
1043 ymin, ymax = ax.get_ylim()
1044 if ymin < 0.0 and ymax > 0.0:
1045 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1046 # Identify unique labels for legend
1047 handles = list(
1048 {
1049 label: handle
1050 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1051 }.values()
1052 )
1053 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1055 # Save plot.
1056 _save_close_figure(fig, "line", filename)
1059def _plot_and_save_line_power_spectrum_series(
1060 cubes: iris.cube.Cube | iris.cube.CubeList,
1061 coords: list[iris.coords.Coord],
1062 ensemble_coord: str,
1063 filename: str,
1064 title: str,
1065 series_coordinate: str,
1066 **kwargs,
1067):
1068 """Plot and save a 1D line series.
1070 Parameters
1071 ----------
1072 cubes: Cube or CubeList
1073 Cube or CubeList containing the cubes to plot on the y-axis.
1074 coords: list[Coord]
1075 Coordinates to plot on the x-axis, one per cube.
1076 ensemble_coord: str
1077 Ensemble coordinate in the cube.
1078 filename: str
1079 Filename of the plot to write.
1080 title: str
1081 Plot title.
1082 series_coordinate: str
1083 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1084 """
1085 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1086 model_colors_map = get_model_colors_map(cubes)
1087 ax = plt.gca()
1089 # Store min/max ranges.
1090 y_levels = []
1092 line_marker = None
1093 line_width = 1
1095 for cube in iter_maybe(cubes):
1096 # next 2 lines replace chunk of code.
1097 xcoord = _select_series_coord(cube, series_coordinate)
1098 xname = xcoord.points
1100 yfield = cube.data # power spectrum
1101 label = None
1102 color = "black"
1103 if model_colors_map: 1103 ↛ 1106line 1103 didn't jump to line 1106 because the condition on line 1103 was always true
1104 label = cube.attributes.get("model_name")
1105 color = model_colors_map.get(label)
1106 for cube_slice in cube.slices_over(ensemble_coord):
1107 # Label with (control) if part of an ensemble or not otherwise.
1108 if cube_slice.coord(ensemble_coord).points == [0]: 1108 ↛ 1122line 1108 didn't jump to line 1122 because the condition on line 1108 was always true
1109 ax.plot(
1110 xname,
1111 yfield,
1112 color=color,
1113 marker=line_marker,
1114 ls="-",
1115 lw=line_width,
1116 label=f"{label} (control)"
1117 if len(cube.coord(ensemble_coord).points) > 1
1118 else label,
1119 )
1120 # Label with (perturbed) if part of an ensemble and not the control.
1121 else:
1122 ax.plot(
1123 xname,
1124 yfield,
1125 color=color,
1126 ls="-",
1127 lw=1.5,
1128 alpha=0.75,
1129 label=f"{label} (member)",
1130 )
1132 # Calculate the global min/max if multiple cubes are given.
1133 _, levels, _ = colorbar_map_levels(cube, axis="y")
1134 if levels is not None: 1134 ↛ 1135line 1134 didn't jump to line 1135 because the condition on line 1134 was never true
1135 y_levels.append(min(levels))
1136 y_levels.append(max(levels))
1138 # Add some labels and tweak the style.
1140 title = f"{title}"
1141 ax.set_title(title, fontsize=16)
1143 # Set appropriate x-axis label based on coordinate
1144 if series_coordinate == "wavelength" or ( 1144 ↛ 1147line 1144 didn't jump to line 1147 because the condition on line 1144 was never true
1145 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1146 ):
1147 ax.set_xlabel("Wavelength (km)", fontsize=14)
1148 elif series_coordinate == "physical_wavenumber" or ( 1148 ↛ 1151line 1148 didn't jump to line 1151 because the condition on line 1148 was never true
1149 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1150 ):
1151 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1152 else: # frequency or check units
1153 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1153 ↛ 1154line 1153 didn't jump to line 1154 because the condition on line 1153 was never true
1154 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1155 else:
1156 ax.set_xlabel("Wavenumber", fontsize=14)
1158 ax.set_ylabel("Power Spectral Density", fontsize=14)
1159 ax.tick_params(axis="both", labelsize=12)
1161 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1163 # Set log-log scale
1164 ax.set_xscale("log")
1165 ax.set_yscale("log")
1167 # Add gridlines
1168 ax.grid(linestyle="--", color="grey", linewidth=1)
1169 # Ientify unique labels for legend
1170 handles = list(
1171 {
1172 label: handle
1173 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1174 }.values()
1175 )
1176 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1178 # Save plot.
1179 _save_close_figure(fig, "line power spectrum", filename)
1182def _plot_and_save_vertical_line_series(
1183 cubes: iris.cube.CubeList,
1184 coords: list[iris.coords.Coord],
1185 ensemble_coord: str,
1186 filename: str,
1187 series_coordinate: str,
1188 title: str,
1189 vmin: float,
1190 vmax: float,
1191 **kwargs,
1192):
1193 """Plot and save a 1D line series in vertical.
1195 Parameters
1196 ----------
1197 cubes: CubeList
1198 1 dimensional Cube or CubeList of the data to plot on x-axis.
1199 coord: list[Coord]
1200 Coordinates to plot on the y-axis, one per cube.
1201 ensemble_coord: str
1202 Ensemble coordinate in the cube.
1203 filename: str
1204 Filename of the plot to write.
1205 series_coordinate: str
1206 Coordinate to use as vertical axis.
1207 title: str
1208 Plot title.
1209 vmin: float
1210 Minimum value for the x-axis.
1211 vmax: float
1212 Maximum value for the x-axis.
1213 """
1214 # plot the vertical pressure axis using log scale
1215 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1217 title_suffix = ""
1218 if cubes and cubes[0].attributes.get("cset_ensemble_mean") == "true": 1218 ↛ 1219line 1218 didn't jump to line 1219 because the condition on line 1218 was never true
1219 title_suffix = " (ensemble mean)"
1221 model_colors_map = get_model_colors_map(cubes)
1223 # Check match-up across sequence coords gives consistent sizes
1224 validate_cubes_coords(cubes, coords)
1226 for cube, coord in zip(cubes, coords, strict=True):
1227 label = None
1228 color = "black"
1229 if model_colors_map: 1229 ↛ 1230line 1229 didn't jump to line 1230 because the condition on line 1229 was never true
1230 label = cube.attributes.get("model_name")
1231 color = model_colors_map.get(label)
1233 if not cube.coords(ensemble_coord): 1233 ↛ 1236line 1233 didn't jump to line 1236 because the condition on line 1233 was never true
1234 # No ensemble coordinate (e.g. after ensemble-mean RMSE): plot
1235 # the cube as a single deterministic line.
1236 iplt.plot(cube, coord, color=color, marker="o", ls="-", lw=3, label=label)
1237 else:
1238 for cube_slice in cube.slices_over(ensemble_coord):
1239 # If ensemble data given plot control member with (control)
1240 # unless single forecast.
1241 if cube_slice.coord(ensemble_coord).points == [0]:
1242 iplt.plot(
1243 cube_slice,
1244 coord,
1245 color=color,
1246 marker="o",
1247 ls="-",
1248 lw=3,
1249 label=f"{label} (control)"
1250 if len(cube.coord(ensemble_coord).points) > 1
1251 else label,
1252 )
1253 # If ensemble data given plot perturbed members with (perturbed).
1254 else:
1255 iplt.plot(
1256 cube_slice,
1257 coord,
1258 color=color,
1259 ls="-",
1260 lw=1.5,
1261 alpha=0.75,
1262 label=f"{label} (member)",
1263 )
1265 # Get the current axis
1266 ax = plt.gca()
1268 # Special handling for pressure level data.
1269 if series_coordinate == "pressure": 1269 ↛ 1291line 1269 didn't jump to line 1291 because the condition on line 1269 was always true
1270 # Invert y-axis and set to log scale.
1271 ax.invert_yaxis()
1272 ax.set_yscale("log")
1274 # Define y-ticks and labels for pressure log axis.
1275 y_tick_labels = [
1276 "1000",
1277 "850",
1278 "700",
1279 "500",
1280 "300",
1281 "200",
1282 "100",
1283 ]
1284 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1286 # Set y-axis limits and ticks.
1287 ax.set_ylim(1100, 100)
1289 # Test if series_coordinate is model level data. The UM data uses
1290 # model_level_number and lfric uses full_levels as coordinate.
1291 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1292 # Define y-ticks and labels for vertical axis.
1293 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1294 y_tick_labels = [str(int(i)) for i in y_ticks]
1295 ax.set_ylim(min(y_ticks), max(y_ticks))
1297 ax.set_yticks(y_ticks)
1298 ax.set_yticklabels(y_tick_labels)
1300 # Set x-axis limits.
1301 ax.set_xlim(vmin, vmax)
1302 # Mark y=0 if present in plot.
1303 if vmin < 0.0 and vmax > 0.0: 1303 ↛ 1304line 1303 didn't jump to line 1304 because the condition on line 1303 was never true
1304 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1306 # Add some labels and tweak the style.
1307 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1308 ax.set_xlabel(
1309 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1310 )
1311 ax.set_title(f"{title}{title_suffix}", fontsize=16)
1312 ax.ticklabel_format(axis="x")
1313 ax.tick_params(axis="y")
1314 ax.tick_params(axis="both", labelsize=12)
1316 # Add gridlines
1317 ax.grid(linestyle="--", color="grey", linewidth=1)
1318 # Ientify unique labels for legend
1319 handles = list(
1320 {
1321 label: handle
1322 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1323 }.values()
1324 )
1325 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1327 # Save plot.
1328 _save_close_figure(fig, "vertical line", filename)
1331def _plot_and_save_scatter_plot(
1332 cube_x: iris.cube.Cube | iris.cube.CubeList,
1333 cube_y: iris.cube.Cube | iris.cube.CubeList,
1334 filename: str,
1335 title: str,
1336 one_to_one: bool,
1337 model_names: list[str] | None = None,
1338 **kwargs,
1339):
1340 """Plot and save a 2D scatter plot.
1342 Parameters
1343 ----------
1344 cube_x: Cube | CubeList
1345 1 dimensional Cube or CubeList of the data to plot on x-axis.
1346 cube_y: Cube | CubeList
1347 1 dimensional Cube or CubeList of the data to plot on y-axis.
1348 filename: str
1349 Filename of the plot to write.
1350 title: str
1351 Plot title.
1352 one_to_one: bool
1353 Whether a 1:1 line is plotted.
1354 """
1355 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1356 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1357 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1358 # over the pairs simultaneously.
1360 # Ensure cube_x and cube_y are iterable
1361 cube_x_iterable = iter_maybe(cube_x)
1362 cube_y_iterable = iter_maybe(cube_y)
1364 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1365 iplt.scatter(cube_x_iter, cube_y_iter)
1366 if one_to_one is True:
1367 plt.plot(
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 [
1373 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1374 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1375 ],
1376 "k",
1377 linestyle="--",
1378 )
1379 ax = plt.gca()
1381 # Add some labels and tweak the style.
1382 if model_names is None:
1383 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1384 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1385 else:
1386 # Add the model names, these should be order of base (x) and other (y).
1387 ax.set_xlabel(
1388 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1389 )
1390 ax.set_ylabel(
1391 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1392 )
1393 ax.set_title(title, fontsize=16)
1394 ax.ticklabel_format(axis="y", useOffset=False)
1395 ax.tick_params(axis="x", labelrotation=15)
1396 ax.tick_params(axis="both", labelsize=12)
1397 ax.autoscale()
1399 # Save plot.
1400 _save_close_figure(fig, "scatter", filename)
1403def _plot_and_save_vector_plot(
1404 cube_u: iris.cube.Cube,
1405 cube_v: iris.cube.Cube,
1406 filename: str,
1407 title: str,
1408 method: Literal["contourf", "pcolormesh"],
1409 **kwargs,
1410):
1411 """Plot and save a 2D vector plot.
1413 Parameters
1414 ----------
1415 cube_u: Cube
1416 2 dimensional Cube of u component of the data.
1417 cube_v: Cube
1418 2 dimensional Cube of v component of the data.
1419 filename: str
1420 Filename of the plot to write.
1421 title: str
1422 Plot title.
1423 """
1424 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1425 # Create a cube containing the magnitude of the vector field.
1426 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1427 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1428 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1429 cube_vec_mag.rename(
1430 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1431 )
1433 # Specify the color bar
1434 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1436 # Setup plot map projection, extent and coastlines and borderlines.
1437 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1439 if method == "contourf":
1440 # Filled contour plot of the field.
1441 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1442 elif method == "pcolormesh":
1443 try:
1444 vmin = min(levels)
1445 vmax = max(levels)
1446 except TypeError:
1447 vmin, vmax = None, None
1448 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1449 # if levels are defined.
1450 if norm is not None:
1451 vmin = None
1452 vmax = None
1453 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1454 else:
1455 raise ValueError(f"Unknown plotting method: {method}")
1457 # Check to see if transect, and if so, adjust y axis.
1458 if is_transect(cube_vec_mag):
1459 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1460 axes.invert_yaxis()
1461 axes.set_yscale("log")
1462 axes.set_ylim(1100, 100)
1463 # If both model_level_number and level_height exists, iplt can construct
1464 # plot as a function of height above orography (NOT sea level).
1465 elif {"model_level_number", "level_height"}.issubset(
1466 {coord.name() for coord in cube_vec_mag.coords()}
1467 ):
1468 axes.set_yscale("log")
1470 axes.set_title(
1471 f"{title}\n"
1472 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1473 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1474 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1475 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1476 fontsize=16,
1477 )
1479 else:
1480 # Add title.
1481 axes.set_title(title, fontsize=16)
1483 # Add watermark with min/max/mean. Currently not user togglable.
1484 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1485 axes.annotate(
1486 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}",
1487 xy=(0.05, -0.05),
1488 xycoords="axes fraction",
1489 xytext=(-5, 5),
1490 textcoords="offset points",
1491 ha="right",
1492 va="bottom",
1493 size=11,
1494 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1495 )
1497 # Add colour bar.
1498 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1499 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1500 # add ticks and tick_labels for every levels if less than 20 levels exist
1501 if levels is not None and len(levels) < 20:
1502 cbar.set_ticks(levels)
1503 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1505 # 30 barbs along the longest axis of the plot, or a barb per point for data
1506 # with less than 30 points.
1507 step = max(max(cube_u.shape) // 30, 1)
1508 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1510 # Save plot.
1511 _save_close_figure(fig, "vector", filename)
1514def _plot_and_save_histogram_series(
1515 cubes: iris.cube.Cube | iris.cube.CubeList,
1516 filename: str,
1517 title: str,
1518 vmin: float,
1519 vmax: float,
1520 **kwargs,
1521):
1522 """Plot and save a histogram series.
1524 Parameters
1525 ----------
1526 cubes: Cube or CubeList
1527 2 dimensional Cube or CubeList of the data to plot as histogram.
1528 filename: str
1529 Filename of the plot to write.
1530 title: str
1531 Plot title.
1532 vmin: float
1533 minimum for colorbar
1534 vmax: float
1535 maximum for colorbar
1536 """
1537 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1538 ax = plt.gca()
1540 model_colors_map = get_model_colors_map(cubes)
1542 # Set default that histograms will produce probability density function
1543 # at each bin (integral over range sums to 1).
1544 density = True
1546 for cube in iter_maybe(cubes):
1547 # Easier to check title (where var name originates)
1548 # than seeing if long names exist etc.
1549 # Exception case, where distribution better fits log scales/bins.
1550 if (
1551 ("surface_microphysical" in title)
1552 or ("rain accumulation" in title)
1553 or ("Rainfall rate Composite" in title)
1554 or ("Nimrod_5min" in title)
1555 ):
1556 if "amount" in title:
1557 # Compute histogram following Klingaman et al. (2017): ASoP
1558 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1559 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1560 density = False
1561 else:
1562 bins = 10.0 ** (
1563 np.arange(-10, 27, 1) / 10.0
1564 ) # Suggestion from RMED toolbox.
1565 bins = np.insert(bins, 0, 0)
1566 ax.set_yscale("log")
1567 vmin = bins[1]
1568 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1569 ax.set_xscale("log")
1570 elif "lightning" in title:
1571 bins = [0, 1, 2, 3, 4, 5]
1572 else:
1573 bins = np.linspace(vmin, vmax, 51)
1574 logger.debug(
1575 "Plotting histogram with %s bins %s - %s.",
1576 np.size(bins),
1577 np.min(bins),
1578 np.max(bins),
1579 )
1581 # Reshape cube data into a single array to allow for a single histogram.
1582 # Otherwise we plot xdim histograms stacked.
1583 cube_data_1d = (cube.data).flatten()
1585 label = None
1586 color = "black"
1587 if model_colors_map:
1588 label = cube.attributes.get("model_name")
1589 color = model_colors_map[label]
1590 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1592 # Compute area under curve.
1593 if (
1594 ("surface_microphysical" in title and "amount" in title)
1595 or ("rain_accumulation" in title)
1596 or ("Rainfall rate Composite" in title)
1597 or ("Nimrod_5min" in title)
1598 ):
1599 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1600 x = x * bin_mean / x.sum()
1601 x = x[1:]
1602 y = y[1:]
1604 ax.plot(
1605 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1606 )
1608 # Add some labels and tweak the style.
1609 ax.set_title(title, fontsize=16)
1610 ax.set_xlabel(
1611 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1612 )
1613 ax.set_ylabel("Normalised probability density", fontsize=14)
1614 if (
1615 ("surface_microphysical" in title and "amount" in title)
1616 or ("rain accumulation" in title)
1617 or ("Nimrod_5min" in title)
1618 ):
1619 ax.set_ylabel(
1620 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1621 )
1622 ax.set_xlim(vmin, vmax)
1623 ax.tick_params(axis="both", labelsize=12)
1625 # Overlay grid-lines onto histogram plot.
1626 ax.grid(linestyle="--", color="grey", linewidth=1)
1627 if model_colors_map:
1628 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1630 # Save plot.
1631 _save_close_figure(fig, "histogram", filename)
1634def _plot_and_save_postage_stamp_histogram_series(
1635 cube: iris.cube.Cube,
1636 filename: str,
1637 title: str,
1638 stamp_coordinate: str,
1639 vmin: float,
1640 vmax: float,
1641 **kwargs,
1642):
1643 """Plot and save postage (ensemble members) stamps for a histogram series.
1645 Parameters
1646 ----------
1647 cube: Cube
1648 2 dimensional Cube of the data to plot as histogram.
1649 filename: str
1650 Filename of the plot to write.
1651 title: str
1652 Plot title.
1653 stamp_coordinate: str
1654 Coordinate that becomes different plots.
1655 vmin: float
1656 minimum for pdf x-axis
1657 vmax: float
1658 maximum for pdf x-axis
1659 """
1660 # Use the smallest square grid that will fit the members.
1661 nmember = len(cube.coord(stamp_coordinate).points)
1662 grid_rows = int(math.sqrt(nmember))
1663 grid_size = math.ceil(nmember / grid_rows)
1665 fig = plt.figure(
1666 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1667 )
1668 # Make a subplot for each member.
1669 for member, subplot in zip(
1670 cube.slices_over(stamp_coordinate),
1671 range(1, grid_size * grid_rows + 1),
1672 strict=False,
1673 ):
1674 # Implicit interface is much easier here, due to needing to have the
1675 # cartopy GeoAxes generated.
1676 plt.subplot(grid_rows, grid_size, subplot)
1677 # Reshape cube data into a single array to allow for a single histogram.
1678 # Otherwise we plot xdim histograms stacked.
1679 member_data_1d = (member.data).flatten()
1680 plt.hist(member_data_1d, density=True, stacked=True)
1681 axes = plt.gca()
1682 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1683 axes.set_title(f"{mtitle}")
1684 axes.set_xlim(vmin, vmax)
1686 # Overall figure title.
1687 fig.suptitle(title, fontsize=16)
1689 # Save plot.
1690 _save_close_figure(fig, "histogram postage stamp", filename)
1693def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1694 cube: iris.cube.Cube,
1695 filename: str,
1696 title: str,
1697 stamp_coordinate: str,
1698 vmin: float,
1699 vmax: float,
1700 **kwargs,
1701):
1702 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1703 ax.set_title(title, fontsize=16)
1704 ax.set_xlim(vmin, vmax)
1705 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1706 ax.set_ylabel("normalised probability density", fontsize=14)
1707 # Loop over all slices along the stamp_coordinate
1708 for member in cube.slices_over(stamp_coordinate):
1709 # Flatten the member data to 1D
1710 member_data_1d = member.data.flatten()
1711 # Plot the histogram using plt.hist
1712 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1713 plt.hist(
1714 member_data_1d,
1715 density=True,
1716 stacked=True,
1717 label=f"{mtitle}",
1718 )
1720 # Add a legend
1721 ax.legend(fontsize=16)
1723 # Save plot.
1724 _save_close_figure(fig, "histogram postage stamp", filename)
1727def _plot_and_save_scatter_series(
1728 cubes: iris.cube.Cube | iris.cube.CubeList,
1729 filename: str,
1730 title: str,
1731 vmin: float,
1732 vmax: float,
1733 hexbin: bool,
1734 **kwargs,
1735):
1736 """Plot and save a scatter plot series.
1738 Parameters
1739 ----------
1740 cubes: Cube or CubeList
1741 2 dimensional Cube or CubeList of the data to plot as scatter.
1742 filename: str
1743 Filename of the plot to write.
1744 title: str
1745 Plot title.
1746 vmin: float
1747 minimum for colorbar
1748 vmax: float
1749 maximum for colorbar
1750 hexbin: bool
1751 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1752 Else scatter of all points, with potential to overplot many comparisons on same plot.
1753 """
1754 if hexbin:
1755 # Check cubes using same functionality as the difference operator.
1756 if len(cubes) != 2:
1757 raise ValueError(
1758 "Cubes should contain exactly 2 cubes for hexbin plotting."
1759 )
1760 title = title.replace("scatter", "hexbin")
1761 filename = filename.replace("scatter", "hexbin")
1763 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1764 ax = plt.gca()
1766 model_colors_map = get_model_colors_map(cubes)
1768 percentiles = np.arange(0, 100, 5)
1769 percentiles[0] = 1
1770 percentiles[-1] = 99
1771 quantiles = iris.cube.CubeList()
1773 # Loop through all output cubes for both data points and overplotting quantiles.
1774 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1775 for plottype in ["points", "quantiles"]:
1776 nplot = 0
1777 for cube in iter_maybe(cubes):
1778 label = None
1779 color = "black"
1780 if model_colors_map: 1780 ↛ 1785line 1780 didn't jump to line 1785 because the condition on line 1780 was always true
1781 label = cube.attributes.get("model_name")
1782 color = model_colors_map[label]
1784 # Plot all data points
1785 if plottype == "points":
1786 if nplot > 0:
1787 if hexbin:
1788 hb = plt.hexbin(
1789 cubes[0].data.flatten(),
1790 cube.data.flatten(),
1791 alpha=0.3,
1792 gridsize=100,
1793 mincnt=1,
1794 )
1795 else:
1796 plt.scatter(
1797 cubes[0].data.flatten(),
1798 cube.data.flatten(),
1799 color=color,
1800 marker="+",
1801 label=None,
1802 alpha=0.3,
1803 )
1805 elif plottype == "quantiles": 1805 ↛ 1824line 1805 didn't jump to line 1824 because the condition on line 1805 was always true
1806 # Construct Q-Q plot
1807 quantiles.append(
1808 cube.collapsed(
1809 cube.coords(dim_coords=True),
1810 iris.analysis.PERCENTILE,
1811 percent=percentiles,
1812 )
1813 )
1814 if nplot > 0:
1815 iplt.scatter(
1816 quantiles[0],
1817 quantiles[-1],
1818 color=color,
1819 marker="o",
1820 label=label,
1821 edgecolors="black",
1822 )
1824 nplot = nplot + 1
1826 # Add some labels and tweak the style.
1827 ax.set_title(title, fontsize=16)
1828 ax.set_xlabel(
1829 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1830 )
1831 ax.set_ylabel(
1832 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1833 )
1834 ax.tick_params(axis="both", labelsize=12)
1835 ax.autoscale()
1837 # Set 1:1 line and equal axes if scatter plot of common cube names
1838 nameA = iter_maybe(cubes)[0].name()
1839 nameB = iter_maybe(cubes)[1].name()
1840 if any(part in nameB.split("_") for part in nameA.split("_")): 1840 ↛ 1851line 1840 didn't jump to line 1851 because the condition on line 1840 was always true
1841 lims = [
1842 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1843 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1844 ]
1845 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1846 ax.set_aspect("equal")
1847 ax.set_xlim(lims)
1848 ax.set_ylim(lims)
1850 # Overlay grid-lines onto scatter plot.
1851 ax.grid(linestyle="--", color="grey", linewidth=1)
1852 if model_colors_map: 1852 ↛ 1856line 1852 didn't jump to line 1856 because the condition on line 1852 was always true
1853 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1855 # Add colorbar if hexbin output
1856 if hexbin:
1857 cb = plt.colorbar(
1858 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1859 )
1860 cb.set_label("Number of data points", size=12)
1862 # Save plot.
1863 _save_close_figure(fig, "scatter", filename)
1866def _spatial_plot(
1867 method: Literal["contourf", "pcolormesh", "scatter"],
1868 cube: iris.cube.Cube,
1869 filename: str | None,
1870 sequence_coordinate: str,
1871 stamp_coordinate: str,
1872 overlay_cube: iris.cube.Cube | None = None,
1873 contour_cube: iris.cube.Cube | None = None,
1874 point_cube: iris.cube.Cube | None = None,
1875 **kwargs,
1876):
1877 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1879 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1880 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1881 is present then postage stamp plots will be produced.
1883 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1884 be overplotted on the same figure.
1886 Parameters
1887 ----------
1888 method: "contourf" | "pcolormesh" | "scatter"
1889 The plotting method to use.
1890 Select choice of "contourf" or "pcolormesh" for gridded data.
1891 Use "scatter" for point-based data.
1892 cube: Cube
1893 Iris cube of the data to plot. It should have two spatial dimensions,
1894 such as lat and lon, and may also have a another two dimension to be
1895 plotted sequentially and/or as postage stamp plots.
1896 filename: str | None
1897 Name of the plot to write, used as a prefix for plot sequences. If None
1898 uses the recipe name.
1899 sequence_coordinate: str
1900 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1901 This coordinate must exist in the cube.
1902 stamp_coordinate: str
1903 Coordinate about which to plot postage stamp plots. Defaults to
1904 ``"realization"``.
1905 overlay_cube: Cube | None, optional
1906 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1907 contour_cube: Cube | None, optional
1908 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1909 point_cube: Cube | None, optional
1910 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
1912 Raises
1913 ------
1914 ValueError
1915 If the cube doesn't have the right dimensions.
1916 TypeError
1917 If the cube isn't a single cube.
1918 """
1919 # Ensure we've got a single cube.
1920 cube = check_single_cube(cube)
1922 # Set title based on recipe metadata or use cube name
1923 recipe_title = get_recipe_metadata().get("title", cube.name())
1925 # Check if there is a valid stamp coordinate in cube dimensions.
1926 if stamp_coordinate == "realization": 1926 ↛ 1931line 1926 didn't jump to line 1931 because the condition on line 1926 was always true
1927 stamp_coordinate = check_stamp_coordinate(cube)
1929 # Make postage stamp plots if stamp_coordinate exists and has more than a
1930 # single point.
1931 plotting_func = _plot_and_save_spatial_plot
1932 try:
1933 if cube.coord(stamp_coordinate).shape[0] > 1:
1934 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1935 except iris.exceptions.CoordinateNotFoundError:
1936 pass
1938 # Produce a geographical scatter plot if the data have a
1939 # dimension called observation or model_obs_error
1940 if any(
1941 crd.var_name == "station"
1942 or crd.var_name == "Station_Name"
1943 or crd.var_name == "model_obs_error"
1944 for crd in cube.coords()
1945 ):
1946 plotting_func = _plot_and_save_spatial_plot
1947 method = "scatter"
1949 # Must have a sequence coordinate.
1950 try:
1951 cube.coord(sequence_coordinate)
1952 except iris.exceptions.CoordinateNotFoundError as err:
1953 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1955 # Create a plot for each value of the sequence coordinate.
1956 plot_index = []
1957 nplot = np.size(cube.coord(sequence_coordinate).points)
1959 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1960 # Set plot titles and filename
1961 seq_coord = cube_slice.coord(sequence_coordinate)
1962 plot_title, plot_filename = _set_title_and_filename(
1963 seq_coord, nplot, recipe_title, filename
1964 )
1966 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1967 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1968 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1969 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1971 # Do the actual plotting.
1972 plotting_func(
1973 cube_slice,
1974 filename=plot_filename,
1975 stamp_coordinate=stamp_coordinate,
1976 title=plot_title,
1977 method=method,
1978 overlay_cube=overlay_slice,
1979 contour_cube=contour_slice,
1980 point_cube=point_slice,
1981 **kwargs,
1982 )
1983 plot_index.append(plot_filename)
1985 # Add list of plots to plot metadata.
1986 complete_plot_index = _append_to_plot_index(plot_index)
1988 # Make a page to display the plots.
1989 _make_plot_html_page(complete_plot_index)
1992####################
1993# Public functions #
1994####################
1997def spatial_contour_plot(
1998 cube: iris.cube.Cube,
1999 filename: str | None = None,
2000 sequence_coordinate: str = "time",
2001 stamp_coordinate: str = "realization",
2002 **kwargs,
2003) -> iris.cube.Cube:
2004 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2006 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2007 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2008 is present then postage stamp plots will be produced.
2010 Parameters
2011 ----------
2012 cube: Cube
2013 Iris cube of the data to plot. It should have two spatial dimensions,
2014 such as lat and lon, and may also have a another two dimension to be
2015 plotted sequentially and/or as postage stamp plots.
2016 filename: str, optional
2017 Name of the plot to write, used as a prefix for plot sequences. Defaults
2018 to the recipe name.
2019 sequence_coordinate: str, optional
2020 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2021 This coordinate must exist in the cube.
2022 stamp_coordinate: str, optional
2023 Coordinate about which to plot postage stamp plots. Defaults to
2024 ``"realization"``.
2026 Returns
2027 -------
2028 Cube
2029 The original cube (so further operations can be applied).
2031 Raises
2032 ------
2033 ValueError
2034 If the cube doesn't have the right dimensions.
2035 TypeError
2036 If the cube isn't a single cube.
2037 """
2038 _spatial_plot(
2039 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2040 )
2041 return cube
2044def spatial_pcolormesh_plot(
2045 cube: iris.cube.Cube,
2046 filename: str | None = None,
2047 sequence_coordinate: str = "time",
2048 stamp_coordinate: str = "realization",
2049 **kwargs,
2050) -> iris.cube.Cube:
2051 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2053 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2054 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2055 is present then postage stamp plots will be produced.
2057 This function is significantly faster than ``spatial_contour_plot``,
2058 especially at high resolutions, and should be preferred unless contiguous
2059 contour areas are important.
2061 Parameters
2062 ----------
2063 cube: Cube
2064 Iris cube of the data to plot. It should have two spatial dimensions,
2065 such as lat and lon, and may also have a another two dimension to be
2066 plotted sequentially and/or as postage stamp plots.
2067 filename: str, optional
2068 Name of the plot to write, used as a prefix for plot sequences. Defaults
2069 to the recipe name.
2070 sequence_coordinate: str, optional
2071 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2072 This coordinate must exist in the cube.
2073 stamp_coordinate: str, optional
2074 Coordinate about which to plot postage stamp plots. Defaults to
2075 ``"realization"``.
2077 Returns
2078 -------
2079 Cube
2080 The original cube (so further operations can be applied).
2082 Raises
2083 ------
2084 ValueError
2085 If the cube doesn't have the right dimensions.
2086 TypeError
2087 If the cube isn't a single cube.
2088 """
2089 _spatial_plot(
2090 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2091 )
2092 return cube
2095def spatial_multi_pcolormesh_plot(
2096 cube: iris.cube.Cube,
2097 overlay_cube: iris.cube.Cube | None = None,
2098 contour_cube: iris.cube.Cube | None = None,
2099 point_cube: iris.cube.Cube | None = None,
2100 filename: str | None = None,
2101 sequence_coordinate: str = "time",
2102 stamp_coordinate: str = "realization",
2103 **kwargs,
2104) -> iris.cube.Cube:
2105 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2107 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2108 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2109 is present then postage stamp plots will be produced.
2111 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2113 If specified, contours of a contour_cube can be overplotted on top of those.
2115 If specified, a spatial scatter map of point_cube can be overplotted.
2117 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2119 This function is significantly faster than ``spatial_contour_plot``,
2120 especially at high resolutions, and should be preferred unless contiguous
2121 contour areas are important.
2123 Parameters
2124 ----------
2125 cube: Cube
2126 Iris cube of the data to plot. It should have two spatial dimensions,
2127 such as lat and lon, and may also have two additional dimensions to be
2128 plotted sequentially and/or as postage stamp plots.
2129 overlay_cube: Cube, optional
2130 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2131 such as lat and lon, and may also have two additional dimensions to be
2132 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.
2133 If not provided, output plot generated without overlay cube.
2134 contour_cube: Cube, optional
2135 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,
2136 such as lat and lon, and may also have two additional dimensions to be
2137 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2138 point_cube: Cube, optional
2139 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
2140 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2141 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2142 postage stamp plots. If not provided, output plot generated without point-based layer.
2143 filename: str, optional
2144 Name of the plot to write, used as a prefix for plot sequences. Defaults
2145 to the recipe name.
2146 sequence_coordinate: str, optional
2147 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2148 This coordinate must exist in the cube.
2149 stamp_coordinate: str, optional
2150 Coordinate about which to plot postage stamp plots. Defaults to
2151 ``"realization"``.
2153 Returns
2154 -------
2155 Cube
2156 The original cube (so further operations can be applied).
2158 Raises
2159 ------
2160 ValueError
2161 If the cube doesn't have the right dimensions.
2162 TypeError
2163 If the cube isn't a single cube.
2164 """
2165 _spatial_plot(
2166 "pcolormesh",
2167 cube,
2168 filename,
2169 sequence_coordinate,
2170 stamp_coordinate,
2171 overlay_cube=overlay_cube,
2172 contour_cube=contour_cube,
2173 point_cube=point_cube,
2174 )
2175 return cube, overlay_cube, contour_cube, point_cube
2178# TODO: Expand function to handle ensemble data.
2179# line_coordinate: str, optional
2180# Coordinate about which to plot multiple lines. Defaults to
2181# ``"realization"``.
2182def plot_line_series(
2183 cube: iris.cube.Cube | iris.cube.CubeList,
2184 filename: str | None = None,
2185 series_coordinate: str = "time",
2186 sequence_coordinate: str = "time",
2187 # add the following for ensembles
2188 stamp_coordinate: str = "realization",
2189 single_plot: bool = False,
2190 **kwargs,
2191) -> iris.cube.Cube | iris.cube.CubeList:
2192 """Plot a line plot for the specified coordinate.
2194 The Cube or CubeList must be 1D.
2196 Parameters
2197 ----------
2198 iris.cube | iris.cube.CubeList
2199 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2200 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2201 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2202 filename: str, optional
2203 Name of the plot to write, used as a prefix for plot sequences. Defaults
2204 to the recipe name.
2205 series_coordinate: str, optional
2206 Coordinate about which to make a series. Defaults to ``"time"``. This
2207 coordinate must exist in the cube.
2209 Returns
2210 -------
2211 iris.cube.Cube | iris.cube.CubeList
2212 The original Cube or CubeList (so further operations can be applied).
2214 Raises
2215 ------
2216 ValueError
2217 If the cubes don't have the right dimensions.
2218 TypeError
2219 If the cube isn't a Cube or CubeList.
2220 """
2221 # Ensure we have a name for the plot file.
2222 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2224 num_models = get_num_models(cube)
2226 validate_cube_shape(cube, num_models)
2228 # Iterate over all cubes and extract coordinate to plot.
2229 cubes = iris.cube.CubeList(iter_maybe(cube))
2230 coords = []
2231 for model_cube in cubes:
2232 try:
2233 coords.append(model_cube.coord(series_coordinate))
2234 except iris.exceptions.CoordinateNotFoundError as err:
2235 raise ValueError(
2236 f"Cube must have a {series_coordinate} coordinate."
2237 ) from err
2238 if model_cube.coords("realization") and model_cube.ndim > 2:
2239 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2241 plot_index = []
2243 # Check if this is a spectral plot by looking for spectral coordinates
2244 is_spectral_plot = series_coordinate in [
2245 "frequency",
2246 "physical_wavenumber",
2247 "wavelength",
2248 ]
2250 if is_spectral_plot:
2251 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2252 # coordinate frequency/wavenumber.
2253 # If several power spectra are plotted with time as sequence_coordinate for the
2254 # time slider option.
2256 # Internal plotting function.
2257 plotting_func = _plot_and_save_line_power_spectrum_series
2259 for model_cube in cubes:
2260 try:
2261 model_cube.coord(sequence_coordinate)
2262 except iris.exceptions.CoordinateNotFoundError as err:
2263 raise ValueError(
2264 f"Cube must have a {sequence_coordinate} coordinate."
2265 ) from err
2267 if num_models == 1: 2267 ↛ 2282line 2267 didn't jump to line 2282 because the condition on line 2267 was always true
2268 # check for ensembles
2269 if ( 2269 ↛ 2273line 2269 didn't jump to line 2273 because the condition on line 2269 was never true
2270 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2271 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2272 ):
2273 if single_plot:
2274 # Plot spectra, mean and ensemble spread on 1 plot
2275 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2276 else:
2277 # Plot postage stamps
2278 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2279 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2280 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2281 else:
2282 all_points = sorted(
2283 set(
2284 itertools.chain.from_iterable(
2285 cb.coord(sequence_coordinate).points for cb in cubes
2286 )
2287 )
2288 )
2289 all_slices = list(
2290 itertools.chain.from_iterable(
2291 cb.slices_over(sequence_coordinate) for cb in cubes
2292 )
2293 )
2294 # Matched slices (matched by seq coord point; it may happen that
2295 # evaluated models do not cover the same seq coord range, hence matching
2296 # necessary)
2297 cube_iterables = [
2298 iris.cube.CubeList(
2299 s
2300 for s in all_slices
2301 if s.coord(sequence_coordinate).points[0] == point
2302 )
2303 for point in all_points
2304 ]
2305 nplot = len(all_points)
2307 # Create a plot for each value of the sequence coordinate. Allowing for
2308 # multiple cubes in a CubeList to be plotted in the same plot for similar
2309 # sequence values. Passing a CubeList into the internal plotting function
2310 # for similar values of the sequence coordinate. cube_slice can be an
2311 # iris.cube.Cube or an iris.cube.CubeList.
2313 for cube_slice in cube_iterables:
2314 # Normalize cube_slice to a list of cubes
2315 if isinstance(cube_slice, iris.cube.CubeList): 2315 ↛ 2316line 2315 didn't jump to line 2316 because the condition on line 2315 was never true
2316 cubes = list(cube_slice)
2317 elif isinstance(cube_slice, iris.cube.Cube): 2317 ↛ 2320line 2317 didn't jump to line 2320 because the condition on line 2317 was always true
2318 cubes = [cube_slice]
2319 else:
2320 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2322 # Use sequence value so multiple sequences can merge.
2323 seq_coord = cube_slice[0].coord(sequence_coordinate)
2324 plot_title, plot_filename = _set_title_and_filename(
2325 seq_coord, nplot, recipe_title, filename
2326 )
2328 # Format the coordinate value in a unit appropriate way.
2329 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2331 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2332 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2332 ↛ 2333line 2332 didn't jump to line 2333 because the condition on line 2332 was never true
2333 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2335 # Do the actual plotting.
2336 plotting_func(
2337 cube_slice,
2338 coords,
2339 stamp_coordinate,
2340 plot_filename,
2341 title,
2342 series_coordinate,
2343 )
2345 plot_index.append(plot_filename)
2346 else:
2347 # Format the title and filename using plotted series coordinate
2348 nplot = 1
2349 seq_coord = coords[0]
2350 plot_title, plot_filename = _set_title_and_filename(
2351 seq_coord, nplot, recipe_title, filename
2352 )
2354 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2355 if (
2356 "station" in [c.name() for c in cubes[0].coords()]
2357 and len(cubes[0].coord("station").points) > 1
2358 ):
2359 for station in cubes[0].coord("station").points:
2360 station_cubes = cubes.extract(iris.Constraint(station=station))
2361 station_name = station_cubes[0].coord("Station_Name").points[0]
2362 station_plotname = plot_filename.replace(
2363 ".png", "_" + station_name + ".png"
2364 )
2365 _plot_and_save_line_series(
2366 station_cubes,
2367 coords,
2368 "realization",
2369 station_plotname,
2370 f"{plot_title} {station_name}",
2371 )
2372 plot_index.append(station_plotname)
2374 else:
2375 # Do the actual plotting for all other series coordinate options.
2376 _plot_and_save_line_series(
2377 cubes, coords, stamp_coordinate, plot_filename, plot_title
2378 )
2380 plot_index.append(plot_filename)
2382 # append plot to list of plots
2383 complete_plot_index = _append_to_plot_index(plot_index)
2385 # Make a page to display the plots.
2386 _make_plot_html_page(complete_plot_index)
2388 return cube
2391def plot_vertical_line_series(
2392 cubes: iris.cube.Cube | iris.cube.CubeList,
2393 filename: str | None = None,
2394 series_coordinate: str = "model_level_number",
2395 sequence_coordinate: str = "time",
2396 # line_coordinate: str = "realization",
2397 **kwargs,
2398) -> iris.cube.Cube | iris.cube.CubeList:
2399 """Plot a line plot against a type of vertical coordinate.
2401 The Cube or CubeList must be 1D.
2403 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2404 then a sequence of plots will be produced.
2406 Parameters
2407 ----------
2408 iris.cube | iris.cube.CubeList
2409 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2410 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2411 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2412 filename: str, optional
2413 Name of the plot to write, used as a prefix for plot sequences. Defaults
2414 to the recipe name.
2415 series_coordinate: str, optional
2416 Coordinate to plot on the y-axis. Can be ``pressure`` or
2417 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2418 for LFRic. Defaults to ``model_level_number``.
2419 This coordinate must exist in the cube.
2420 sequence_coordinate: str, optional
2421 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2422 This coordinate must exist in the cube.
2424 Returns
2425 -------
2426 iris.cube.Cube | iris.cube.CubeList
2427 The original Cube or CubeList (so further operations can be applied).
2428 Plotted data.
2430 Raises
2431 ------
2432 ValueError
2433 If the cubes doesn't have the right dimensions.
2434 TypeError
2435 If the cube isn't a Cube or CubeList.
2436 """
2437 # Ensure we have a name for the plot file.
2438 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2440 cubes = iter_maybe(cubes)
2441 # Initialise empty list to hold all data from all cubes in a CubeList
2442 all_data = []
2444 # Store min/max ranges for x range.
2445 x_levels = []
2447 num_models = get_num_models(cubes)
2449 validate_cube_shape(cubes, num_models)
2451 # Iterate over all cubes in cube or CubeList and plot.
2452 coords = []
2453 for cube in cubes:
2454 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2455 try:
2456 coords.append(cube.coord(series_coordinate))
2457 except iris.exceptions.CoordinateNotFoundError as err:
2458 raise ValueError(
2459 f"Cube must have a {series_coordinate} coordinate."
2460 ) from err
2462 try:
2463 if cube.ndim > 1 or not cube.coords("realization"): 2463 ↛ 2471line 2463 didn't jump to line 2471 because the condition on line 2463 was always true
2464 cube.coord(sequence_coordinate)
2465 except iris.exceptions.CoordinateNotFoundError as err:
2466 raise ValueError(
2467 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2468 ) from err
2470 # Get minimum and maximum from levels information.
2471 _, levels, _ = colorbar_map_levels(cube, axis="x")
2472 if levels is not None: 2472 ↛ 2476line 2472 didn't jump to line 2476 because the condition on line 2472 was always true
2473 x_levels.append(min(levels))
2474 x_levels.append(max(levels))
2475 else:
2476 all_data.append(cube.data)
2478 if len(x_levels) == 0: 2478 ↛ 2480line 2478 didn't jump to line 2480 because the condition on line 2478 was never true
2479 # Combine all data into a single NumPy array
2480 combined_data = np.concatenate(all_data)
2482 # Set the lower and upper limit for the x-axis to ensure all plots have
2483 # same range. This needs to read the whole cube over the range of the
2484 # sequence and if applicable postage stamp coordinate.
2485 vmin = np.floor(combined_data.min())
2486 vmax = np.ceil(combined_data.max())
2487 else:
2488 vmin = min(x_levels)
2489 vmax = max(x_levels)
2491 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2492 # a single profile directly without iterating over a sequence.
2493 sequence_coords = [
2494 cube.coord(sequence_coordinate)
2495 for cube in cubes
2496 if cube.coords(sequence_coordinate)
2497 ]
2498 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2499 np.size(coord.points) > 1 for coord in sequence_coords
2500 )
2501 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2502 np.size(coord.points) == 1 for coord in sequence_coords
2503 )
2505 plot_index = []
2506 if has_sequence_coord: 2506 ↛ 2531line 2506 didn't jump to line 2531 because the condition on line 2506 was always true
2507 # Matching the slices (matching by seq coord point; it may happen that
2508 # evaluated models do not cover the same seq coord range, hence matching
2509 # necessary)
2510 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2511 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2512 for cubes_slice in cube_iterables:
2513 # Format the coordinate value in a unit appropriate way.
2514 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2515 plot_title, plot_filename = _set_title_and_filename(
2516 seq_coord, nplot, recipe_title, filename
2517 )
2519 # Do the actual plotting.
2520 _plot_and_save_vertical_line_series(
2521 cubes_slice,
2522 coords,
2523 "realization",
2524 plot_filename,
2525 series_coordinate,
2526 title=plot_title,
2527 vmin=vmin,
2528 vmax=vmax,
2529 )
2530 plot_index.append(plot_filename)
2531 elif has_scalar_sequence_coord:
2532 # Scalar sequence coordinate (typically aggregated time bounds):
2533 # make one plot and include sequence period in title/filename.
2534 plot_title, plot_filename = _set_title_and_filename(
2535 sequence_coords[0], 1, recipe_title, filename
2536 )
2538 _plot_and_save_vertical_line_series(
2539 cubes,
2540 coords,
2541 "realization",
2542 plot_filename,
2543 series_coordinate,
2544 title=plot_title,
2545 vmin=vmin,
2546 vmax=vmax,
2547 )
2548 plot_index.append(plot_filename)
2549 else:
2550 # 1D case: no sequence coordinate, plot a single profile.
2551 plot_title = recipe_title
2552 if filename:
2553 plot_filename = filename
2554 else:
2555 plot_filename = f"{slugify(plot_title)}.png"
2557 _plot_and_save_vertical_line_series(
2558 cubes,
2559 coords,
2560 "realization",
2561 plot_filename,
2562 series_coordinate,
2563 title=plot_title,
2564 vmin=vmin,
2565 vmax=vmax,
2566 )
2567 plot_index.append(plot_filename)
2569 # Add list of plots to plot metadata.
2570 complete_plot_index = _append_to_plot_index(plot_index)
2572 # Make a page to display the plots.
2573 _make_plot_html_page(complete_plot_index)
2575 return cubes
2578def qq_plot(
2579 cubes: iris.cube.CubeList,
2580 coordinates: list[str],
2581 percentiles: list[float],
2582 model_names: list[str],
2583 filename: str | None = None,
2584 one_to_one: bool = True,
2585 **kwargs,
2586) -> iris.cube.CubeList:
2587 """Plot a Quantile-Quantile plot between two models for common time points.
2589 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2590 collapsed within the operator over all specified coordinates such as
2591 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2592 ensemble members to ensure a 1D cube (array).
2594 Parameters
2595 ----------
2596 cubes: iris.cube.CubeList
2597 Two cubes of the same variable with different models.
2598 coordinate: list[str]
2599 The list of coordinates to collapse over. This list should be
2600 every coordinate within the cube to result in a 1D cube around
2601 the percentile coordinate.
2602 percent: list[float]
2603 A list of percentiles to appear in the plot.
2604 model_names: list[str]
2605 A list of model names to appear on the axis of the plot.
2606 filename: str, optional
2607 Filename of the plot to write.
2608 one_to_one: bool, optional
2609 If True a 1:1 line is plotted; if False it is not. Default is True.
2611 Raises
2612 ------
2613 ValueError
2614 When the cubes are not compatible.
2616 Notes
2617 -----
2618 The quantile-quantile plot is a variant on the scatter plot representing
2619 two datasets by their quantiles (percentiles) for common time points.
2620 This plot does not use a theoretical distribution to compare against, but
2621 compares percentiles of two datasets. This plot does
2622 not use all raw data points, but plots the selected percentiles (quantiles) of
2623 each variable instead for the two datasets, thereby normalising the data for a
2624 direct comparison between the selected percentiles of the two dataset distributions.
2626 Quantile-quantile plots are valuable for comparing against
2627 observations and other models. Identical percentiles between the variables
2628 will lie on the one-to-one line implying the values correspond well to each
2629 other. Where there is a deviation from the one-to-one line a range of
2630 possibilities exist depending on how and where the data is shifted (e.g.,
2631 Wilks 2011 [Wilks2011]_).
2633 For distributions above the one-to-one line the distribution is left-skewed;
2634 below is right-skewed. A distinct break implies a bimodal distribution, and
2635 closer values/values further apart at the tails imply poor representation of
2636 the extremes.
2638 References
2639 ----------
2640 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2641 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2642 """
2643 # Check cubes using same functionality as the difference operator.
2644 if len(cubes) != 2:
2645 raise ValueError("cubes should contain exactly 2 cubes.")
2646 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2647 other: Cube = cubes.extract_cube(
2648 iris.Constraint(
2649 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2650 )
2651 )
2653 # Get spatial coord names.
2654 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2655 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2657 # Ensure cubes to compare are on common differencing grid.
2658 # This is triggered if either
2659 # i) latitude and longitude shapes are not the same. Note grid points
2660 # are not compared directly as these can differ through rounding
2661 # errors.
2662 # ii) or variables are known to often sit on different grid staggering
2663 # in different models (e.g. cell center vs cell edge), as is the case
2664 # for UM and LFRic comparisons.
2665 # In future greater choice of regridding method might be applied depending
2666 # on variable type. Linear regridding can in general be appropriate for smooth
2667 # variables. Care should be taken with interpretation of differences
2668 # given this dependency on regridding.
2669 if (
2670 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2671 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2672 ) or (
2673 base.long_name
2674 in [
2675 "eastward_wind_at_10m",
2676 "northward_wind_at_10m",
2677 "northward_wind_at_cell_centres",
2678 "eastward_wind_at_cell_centres",
2679 "zonal_wind_at_pressure_levels",
2680 "meridional_wind_at_pressure_levels",
2681 "potential_vorticity_at_pressure_levels",
2682 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2683 ]
2684 ):
2685 logger.debug("Linear regridding base cube to other grid to compute differences")
2686 base = regrid_onto_cube(base, other, method="Linear")
2688 # Extract just common time points.
2689 base, other = _extract_common_time_points(base, other)
2691 # Equalise attributes so we can merge.
2692 fully_equalise_attributes([base, other])
2693 logger.debug("Base: %s\nOther: %s", base, other)
2695 # Collapse cubes.
2696 base = collapse(
2697 base,
2698 coordinate=coordinates,
2699 method="PERCENTILE",
2700 additional_percent=percentiles,
2701 )
2702 other = collapse(
2703 other,
2704 coordinate=coordinates,
2705 method="PERCENTILE",
2706 additional_percent=percentiles,
2707 )
2709 # Ensure we have a name for the plot file.
2710 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2711 title = f"{recipe_title}"
2713 if filename is None:
2714 filename = slugify(recipe_title)
2716 # Add file extension.
2717 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2719 # Do the actual plotting on a scatter plot
2720 _plot_and_save_scatter_plot(
2721 base, other, plot_filename, title, one_to_one, model_names
2722 )
2724 # Add list of plots to plot metadata.
2725 plot_index = _append_to_plot_index([plot_filename])
2727 # Make a page to display the plots.
2728 _make_plot_html_page(plot_index)
2730 return iris.cube.CubeList([base, other])
2733def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2734 """
2735 Plot a Hinton style triangle/scorecard plot.
2737 This plot type can be useful for summarising high level information, such as comparing
2738 how 'skillful' two models are when verified against observations for a variety of metrics,
2739 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2740 than customisable by the user as input arguments; many have been designed to automatically
2741 scale the plot depending on the number of x and y components.
2743 Parameters
2744 ----------
2745 change: np.ndarray
2746 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2747 size/direction.
2748 signif: np.ndarray
2749 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2750 xaxis_labels: list
2751 List of labels for the xaxis (must match the second dimension length of signif and change,
2752 along with magnitude if not None).
2753 yaxis_labels: list
2754 List of labels for the yaxis (must match the first dimension length of signif and change,
2755 along with magnitude if not None).
2756 magnitude: np.ndarray | None
2757 Optional 2D array, matching the shape of change, signif, which contains numerical values
2758 the user wishes to display under each respective triangle.
2760 Returns
2761 -------
2762 matplotlib axes object to either display or do further modifications to.
2763 """
2764 # Setup colors of triangles
2765 color_pos = "#7CAE00"
2766 color_neg = "#7B68EE"
2768 # Setup cell/text size ratios
2769 figsize = None
2770 cell_size_in = 0.35
2771 text_row_ratio = 0.25
2773 # Ensure arrays, and change to bool for sig.
2774 change = np.asarray(change)
2775 signif = np.asarray(signif).astype(bool)
2776 if magnitude is not None: 2776 ↛ 2777line 2776 didn't jump to line 2777 because the condition on line 2776 was never true
2777 magnitude = np.asarray(magnitude)
2779 # Get the number of x and y elements
2780 ny, nx = change.shape
2782 # Build non-uniform y coordinates
2783 tri_height = 1.0
2784 txt_height = text_row_ratio
2786 tri_y = []
2787 txt_y = []
2788 y_edges = [0.0]
2790 y = 0.0
2791 for _j in range(ny):
2792 tri_y.append(y + tri_height / 2)
2793 y += tri_height
2794 y_edges.append(y)
2796 if magnitude is not None: 2796 ↛ 2797line 2796 didn't jump to line 2797 because the condition on line 2796 was never true
2797 txt_y.append(y + txt_height / 2)
2798 y += txt_height
2799 y_edges.append(y)
2801 total_height = y
2803 # Dynamic figure size
2804 if figsize is None: 2804 ↛ 2809line 2804 didn't jump to line 2809 because the condition on line 2804 was always true
2805 width = nx * cell_size_in
2806 height = total_height * cell_size_in + 2
2807 figsize = (width, height)
2809 fig, ax = plt.subplots(figsize=figsize)
2811 # Setup axes and grid.
2812 ax.set_aspect("equal", adjustable="box")
2813 ax.set_xlim(-0.5, nx - 0.5)
2814 ax.set_ylim(0, total_height)
2816 ax.set_xticks(np.arange(nx))
2817 ax.set_xticklabels(xaxis_labels, rotation=90)
2819 ax.set_yticks(tri_y)
2820 ax.set_yticklabels(yaxis_labels)
2822 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2823 ax.set_yticks(y_edges, minor=True)
2825 ax.set_axisbelow(True)
2826 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2827 ax.grid(False, which="major")
2828 ax.tick_params(which="minor", length=0)
2830 ax.invert_yaxis()
2832 # Compute marker scaling (fixed overlap)
2833 fig.canvas.draw()
2835 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2836 width_in, height_in = bbox.width, bbox.height
2838 cell_w = (width_in * fig.dpi) / nx
2839 cell_h = (height_in * fig.dpi) / total_height
2840 cell_pixels = min(cell_w, cell_h)
2842 max_marker_size = (0.6 * cell_pixels) ** 2
2844 text_fontsize = cell_pixels * 0.15
2846 # Plot triangles + text
2847 for j in range(ny):
2848 for i in range(nx):
2849 val = change[j, i]
2850 if np.isnan(val): 2850 ↛ 2851line 2850 didn't jump to line 2851 because the condition on line 2850 was never true
2851 continue
2853 if abs(val) < 0.01: 2853 ↛ 2854line 2853 didn't jump to line 2854 because the condition on line 2853 was never true
2854 continue
2856 sig = signif[j, i]
2857 size = max_marker_size * abs(val)
2859 # Triangle style
2860 if val >= 0:
2861 marker = "^"
2862 color = color_pos
2863 else:
2864 marker = "v"
2865 color = color_neg
2867 if sig:
2868 edgecolor = "black"
2869 linewidth = 0.6
2870 else:
2871 edgecolor = "none"
2872 linewidth = 0.0
2874 # Triangle
2875 ax.scatter(
2876 i,
2877 tri_y[j],
2878 s=size,
2879 marker=marker,
2880 c=color,
2881 edgecolors=edgecolor,
2882 linewidths=linewidth,
2883 zorder=3,
2884 clip_on=True, # ensures no rendering bleed
2885 )
2887 # Text row
2888 if magnitude is not None: 2888 ↛ 2889line 2888 didn't jump to line 2889 because the condition on line 2888 was never true
2889 mag_val = magnitude[j, i]
2891 if not np.isnan(mag_val):
2892 ax.text(
2893 i,
2894 txt_y[j],
2895 f"{mag_val:.1f}",
2896 ha="center",
2897 va="center",
2898 fontsize=text_fontsize,
2899 color="black",
2900 zorder=4,
2901 )
2903 plt.tight_layout()
2904 return fig, ax
2907def scatter_plot(
2908 cube_x: iris.cube.Cube | iris.cube.CubeList,
2909 cube_y: iris.cube.Cube | iris.cube.CubeList,
2910 filename: str | None = None,
2911 one_to_one: bool = True,
2912 **kwargs,
2913) -> iris.cube.CubeList:
2914 """Plot a scatter plot between two variables.
2916 Both cubes must be 1D.
2918 Parameters
2919 ----------
2920 cube_x: Cube | CubeList
2921 1 dimensional Cube of the data to plot on y-axis.
2922 cube_y: Cube | CubeList
2923 1 dimensional Cube of the data to plot on x-axis.
2924 filename: str, optional
2925 Filename of the plot to write.
2926 one_to_one: bool, optional
2927 If True a 1:1 line is plotted; if False it is not. Default is True.
2929 Returns
2930 -------
2931 cubes: CubeList
2932 CubeList of the original x and y cubes for further processing.
2934 Raises
2935 ------
2936 ValueError
2937 If the cube doesn't have the right dimensions and cubes not the same
2938 size.
2939 TypeError
2940 If the cube isn't a single cube.
2942 Notes
2943 -----
2944 Scatter plots are used for determining if there is a relationship between
2945 two variables. Positive relations have a slope going from bottom left to top
2946 right; Negative relations have a slope going from top left to bottom right.
2947 """
2948 # Iterate over all cubes in cube or CubeList and plot.
2949 for cube_iter in iter_maybe(cube_x):
2950 # Check cubes are correct shape.
2951 cube_iter = check_single_cube(cube_iter)
2952 if cube_iter.ndim > 1:
2953 raise ValueError("cube_x must be 1D.")
2955 # Iterate over all cubes in cube or CubeList and plot.
2956 for cube_iter in iter_maybe(cube_y):
2957 # Check cubes are correct shape.
2958 cube_iter = check_single_cube(cube_iter)
2959 if cube_iter.ndim > 1:
2960 raise ValueError("cube_y must be 1D.")
2962 # Ensure we have a name for the plot file.
2963 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2964 title = f"{recipe_title}"
2966 if filename is None:
2967 filename = slugify(recipe_title)
2969 # Add file extension.
2970 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2972 # Do the actual plotting.
2973 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2975 # Add list of plots to plot metadata.
2976 plot_index = _append_to_plot_index([plot_filename])
2978 # Make a page to display the plots.
2979 _make_plot_html_page(plot_index)
2981 return iris.cube.CubeList([cube_x, cube_y])
2984def vector_plot(
2985 cube_u: iris.cube.Cube,
2986 cube_v: iris.cube.Cube,
2987 filename: str | None = None,
2988 sequence_coordinate: str = "time",
2989 **kwargs,
2990) -> iris.cube.CubeList:
2991 """Plot a vector plot based on the input u and v components."""
2992 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2994 # Cubes must have a matching sequence coordinate.
2995 try:
2996 # Check that the u and v cubes have the same sequence coordinate.
2997 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2997 ↛ anywhereline 2997 didn't jump anywhere: it always raised an exception.
2998 raise ValueError("Coordinates do not match.")
2999 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3000 raise ValueError(
3001 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3002 ) from err
3004 # Create a plot for each value of the sequence coordinate.
3005 plot_index = []
3006 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3007 for cube_u_slice, cube_v_slice in zip(
3008 cube_u.slices_over(sequence_coordinate),
3009 cube_v.slices_over(sequence_coordinate),
3010 strict=True,
3011 ):
3012 # Format the coordinate value in a unit appropriate way.
3013 seq_coord = cube_u_slice.coord(sequence_coordinate)
3014 plot_title, plot_filename = _set_title_and_filename(
3015 seq_coord, nplot, recipe_title, filename
3016 )
3018 # Do the actual plotting.
3019 _plot_and_save_vector_plot(
3020 cube_u_slice,
3021 cube_v_slice,
3022 filename=plot_filename,
3023 title=plot_title,
3024 method="pcolormesh",
3025 )
3026 plot_index.append(plot_filename)
3028 # Add list of plots to plot metadata.
3029 complete_plot_index = _append_to_plot_index(plot_index)
3031 # Make a page to display the plots.
3032 _make_plot_html_page(complete_plot_index)
3034 return iris.cube.CubeList([cube_u, cube_v])
3037def plot_histogram_series(
3038 cubes: iris.cube.Cube | iris.cube.CubeList,
3039 filename: str | None = None,
3040 sequence_coordinate: str = "time",
3041 stamp_coordinate: str = "realization",
3042 single_plot: bool = False,
3043 **kwargs,
3044) -> iris.cube.Cube | iris.cube.CubeList:
3045 """Plot a histogram plot for each vertical level provided.
3047 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3048 is present then a sequence of plots will be produced using the time slider
3049 functionality to scroll through histograms against time. If a
3050 stamp_coordinate is present then postage stamp plots will be produced. If
3051 stamp_coordinate and single_plot is True, all postage stamp plots will be
3052 plotted in a single plot instead of separate postage stamp plots.
3054 Parameters
3055 ----------
3056 cubes: Cube | iris.cube.CubeList
3057 Iris cube or CubeList of the data to plot. It should have a single dimension other
3058 than the stamp coordinate.
3059 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3060 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3061 filename: str, optional
3062 Name of the plot to write, used as a prefix for plot sequences. Defaults
3063 to the recipe name.
3064 sequence_coordinate: str, optional
3065 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3066 This coordinate must exist in the cube and will be used for the time
3067 slider.
3068 stamp_coordinate: str, optional
3069 Coordinate about which to plot postage stamp plots. Defaults to
3070 ``"realization"``.
3071 single_plot: bool, optional
3072 If True, all postage stamp plots will be plotted in a single plot. If
3073 False, each postage stamp plot will be plotted separately. Is only valid
3074 if stamp_coordinate exists and has more than a single point.
3076 Returns
3077 -------
3078 iris.cube.Cube | iris.cube.CubeList
3079 The original Cube or CubeList (so further operations can be applied).
3080 Plotted data.
3082 Raises
3083 ------
3084 ValueError
3085 If the cube doesn't have the right dimensions.
3086 TypeError
3087 If the cube isn't a Cube or CubeList.
3088 """
3089 recipe_title = get_recipe_metadata().get("title", "Histogram")
3091 cubes = iter_maybe(cubes)
3093 # Internal plotting function.
3094 plotting_func = _plot_and_save_histogram_series
3096 num_models = get_num_models(cubes)
3098 validate_cube_shape(cubes, num_models)
3100 # If several histograms are plotted, check sequence_coordinate
3101 check_sequence_coordinate(cubes, sequence_coordinate)
3103 # Get axis minimum and maximum from levels information.
3104 # If no levels set, derive minima and maxima from data in CubeList.
3105 vmin, vmax = _set_axis_range(cubes)
3107 # Make postage stamp plots if stamp_coordinate exists and has more than a
3108 # single point. If single_plot is True:
3109 # -- all postage stamp plots will be plotted in a single plot instead of
3110 # separate postage stamp plots.
3111 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3112 # produced per single model only
3113 if num_models == 1:
3114 if ( 3114 ↛ 3118line 3114 didn't jump to line 3118 because the condition on line 3114 was never true
3115 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3116 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3117 ):
3118 if single_plot:
3119 plotting_func = (
3120 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3121 )
3122 else:
3123 plotting_func = _plot_and_save_postage_stamp_histogram_series
3124 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3125 else:
3126 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3128 plot_index = []
3129 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3130 # Create a plot for each value of the sequence coordinate. Allowing for
3131 # multiple cubes in a CubeList to be plotted in the same plot for similar
3132 # sequence values. Passing a CubeList into the internal plotting function
3133 # for similar values of the sequence coordinate. cube_slice can be an
3134 # iris.cube.Cube or an iris.cube.CubeList.
3135 for cube_slice in cube_iterables:
3136 single_cube = cube_slice
3137 if isinstance(cube_slice, iris.cube.CubeList):
3138 single_cube = cube_slice[0]
3140 # Ensure valid stamp coordinate in cube dimensions
3141 if stamp_coordinate == "realization": 3141 ↛ 3144line 3141 didn't jump to line 3144 because the condition on line 3141 was always true
3142 stamp_coordinate = check_stamp_coordinate(single_cube)
3143 # Set plot titles and filename, based on sequence coordinate
3144 seq_coord = single_cube.coord(sequence_coordinate)
3145 # Use time coordinate in title and filename if single histogram output.
3146 if sequence_coordinate == "realization" and nplot == 1: 3146 ↛ 3147line 3146 didn't jump to line 3147 because the condition on line 3146 was never true
3147 seq_coord = single_cube.coord("time")
3148 # Use station name in title and filename if model vs obs comparison
3149 if sequence_coordinate == "station": 3149 ↛ 3150line 3149 didn't jump to line 3150 because the condition on line 3149 was never true
3150 seq_coord = single_cube.coord("Station_Name")
3152 plot_title, plot_filename = _set_title_and_filename(
3153 seq_coord, nplot, recipe_title, filename
3154 )
3156 # Do the actual plotting.
3157 plotting_func(
3158 cube_slice,
3159 filename=plot_filename,
3160 stamp_coordinate=stamp_coordinate,
3161 title=plot_title,
3162 vmin=vmin,
3163 vmax=vmax,
3164 )
3165 plot_index.append(plot_filename)
3167 # Add list of plots to plot metadata.
3168 complete_plot_index = _append_to_plot_index(plot_index)
3170 # Make a page to display the plots.
3171 _make_plot_html_page(complete_plot_index)
3173 return cubes
3176def plot_scatter_series(
3177 cubes: iris.cube.Cube | iris.cube.CubeList,
3178 filename: str | None = None,
3179 sequence_coordinate: str = "time",
3180 stamp_coordinate: str = "realization",
3181 hexbin: bool = False,
3182 **kwargs,
3183) -> iris.cube.Cube | iris.cube.CubeList:
3184 """Plot a scatter plot for each sequence coordinate provided.
3186 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3187 is present then a sequence of plots will be produced using the time slider
3188 functionality to scroll through scatter against time. If a
3189 stamp_coordinate is present then postage stamp plots will be produced. If
3190 stamp_coordinate and single_plot is True, all postage stamp plots will be
3191 plotted in a single plot instead of separate postage stamp plots.
3193 Parameters
3194 ----------
3195 cubes: Cube | iris.cube.CubeList
3196 Iris cube or CubeList of the data to plot. It should have a single dimension other
3197 than the stamp coordinate.
3198 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3199 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3200 filename: str, optional
3201 Name of the plot to write, used as a prefix for plot sequences. Defaults
3202 to the recipe name.
3203 sequence_coordinate: str, optional
3204 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3205 This coordinate must exist in the cube and will be used for the time
3206 slider.
3207 stamp_coordinate: str, optional
3208 Coordinate about which to plot postage stamp plots. Defaults to
3209 ``"realization"``.
3210 hexbin: bool, optional
3211 If True, generate hexbin comparison plot.
3212 If False, generate point-by-point scatter plot.
3214 Returns
3215 -------
3216 iris.cube.Cube | iris.cube.CubeList
3217 The original Cube or CubeList (so further operations can be applied).
3218 Plotted data.
3220 Raises
3221 ------
3222 ValueError
3223 If the cube doesn't have the right dimensions.
3224 TypeError
3225 If the cube isn't a Cube or CubeList.
3226 """
3227 recipe_title = get_recipe_metadata().get("title", "Scatter")
3229 cubes = iter_maybe(cubes)
3231 # Internal plotting function.
3232 plotting_func = _plot_and_save_scatter_series
3234 num_models = get_num_models(cubes)
3236 validate_cube_shape(cubes, num_models)
3238 check_sequence_coordinate(cubes, sequence_coordinate)
3240 vmin, vmax = _set_axis_range(cubes)
3242 # Require >1 models to compare on scatter plot
3243 if num_models > 1:
3244 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3245 else:
3246 raise ValueError(
3247 "Scatter plot series requires multiple number of models in input data."
3248 )
3250 plot_index = []
3251 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3252 # Create a plot for each value of the sequence coordinate. Allowing for
3253 # multiple cubes in a CubeList to be plotted in the same plot for similar
3254 # sequence values. Passing a CubeList into the internal plotting function
3255 # for similar values of the sequence coordinate. cube_slice can be an
3256 # iris.cube.Cube or an iris.cube.CubeList.
3257 for cube_slice in cube_iterables:
3258 single_cube = cube_slice
3259 if isinstance(cube_slice, iris.cube.CubeList): 3259 ↛ 3263line 3259 didn't jump to line 3263 because the condition on line 3259 was always true
3260 single_cube = cube_slice[0]
3262 # Ensure valid stamp coordinate in cube dimensions
3263 if stamp_coordinate == "realization": 3263 ↛ 3266line 3263 didn't jump to line 3266 because the condition on line 3263 was always true
3264 stamp_coordinate = check_stamp_coordinate(single_cube)
3265 # Set plot titles and filename, based on sequence coordinate
3266 seq_coord = single_cube.coord(sequence_coordinate)
3267 # Use time coordinate in title and filename if single histogram output.
3268 if sequence_coordinate == "realization" and nplot == 1:
3269 seq_coord = single_cube.coord("time")
3270 # Use station name in title and filename if model vs obs comparison
3271 if sequence_coordinate == "station":
3272 seq_coord = single_cube.coord("Station_Name")
3274 plot_title, plot_filename = _set_title_and_filename(
3275 seq_coord, nplot, recipe_title, filename
3276 )
3278 # Do the actual plotting.
3279 plotting_func(
3280 cube_slice,
3281 filename=plot_filename,
3282 stamp_coordinate=stamp_coordinate,
3283 title=plot_title,
3284 vmin=vmin,
3285 vmax=vmax,
3286 hexbin=hexbin,
3287 )
3288 plot_index.append(plot_filename)
3290 # Add list of plots to plot metadata.
3291 complete_plot_index = _append_to_plot_index(plot_index)
3293 # Make a page to display the plots.
3294 _make_plot_html_page(complete_plot_index)
3296 return cubes
3299def _plot_and_save_postage_stamp_power_spectrum_series(
3300 cubes: iris.cube.Cube,
3301 coords: list[iris.coords.Coord],
3302 stamp_coordinate: str,
3303 filename: str,
3304 title: str,
3305 series_coordinate: str | None = None,
3306 **kwargs,
3307):
3308 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3310 Parameters
3311 ----------
3312 cubes: Cube or CubeList
3313 Cube or Cubelist of the power spectrum data.
3314 coords: list[Coord]
3315 Coordinates to plot on the x-axis, one per cube.
3316 stamp_coordinate: str
3317 Coordinate that becomes different plots.
3318 filename: str
3319 Filename of the plot to write.
3320 title: str
3321 Plot title.
3322 series_coordinate: str, optional
3323 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3325 """
3326 # Use the smallest square grid that will fit the members.
3327 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3329 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3330 model_colors_map = get_model_colors_map(cubes)
3331 # ax = plt.gca()
3332 # Make a subplot for each member.
3333 for member, subplot in zip(
3334 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3335 ):
3336 ax = plt.subplot(grid_size, grid_size, subplot)
3338 # Store min/max ranges.
3339 y_levels = []
3341 line_marker = None
3342 line_width = 1
3344 for cube in iter_maybe(member):
3345 xcoord = _select_series_coord(cube, series_coordinate)
3346 xname = xcoord.points
3348 yfield = cube.data # power spectrum
3349 label = None
3350 color = "black"
3351 if model_colors_map: 3351 ↛ 3352line 3351 didn't jump to line 3352 because the condition on line 3351 was never true
3352 label = cube.attributes.get("model_name")
3353 color = model_colors_map.get(label)
3355 if member.coord(stamp_coordinate).points == [0]:
3356 ax.plot(
3357 xname,
3358 yfield,
3359 color=color,
3360 marker=line_marker,
3361 ls="-",
3362 lw=line_width,
3363 label=f"{label} (control)"
3364 if len(cube.coord(stamp_coordinate).points) > 1
3365 else label,
3366 )
3367 # Label with member if part of an ensemble and not the control.
3368 else:
3369 ax.plot(
3370 xname,
3371 yfield,
3372 color=color,
3373 ls="-",
3374 lw=1.5,
3375 alpha=0.75,
3376 label=f"{label} (member)",
3377 )
3379 # Calculate the global min/max if multiple cubes are given.
3380 _, levels, _ = colorbar_map_levels(cube, axis="y")
3381 if levels is not None: 3381 ↛ 3382line 3381 didn't jump to line 3382 because the condition on line 3381 was never true
3382 y_levels.append(min(levels))
3383 y_levels.append(max(levels))
3385 # Add some labels and tweak the style.
3386 title = f"{title}"
3387 ax.set_title(title, fontsize=16)
3389 # Set appropriate x-axis label based on coordinate
3390 if series_coordinate == "wavelength" or ( 3390 ↛ 3393line 3390 didn't jump to line 3393 because the condition on line 3390 was never true
3391 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3392 ):
3393 ax.set_xlabel("Wavelength (km)", fontsize=14)
3394 elif series_coordinate == "physical_wavenumber" or ( 3394 ↛ 3399line 3394 didn't jump to line 3399 because the condition on line 3394 was always true
3395 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3396 ):
3397 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3398 else: # frequency or check units
3399 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3400 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3401 else:
3402 ax.set_xlabel("Wavenumber", fontsize=14)
3404 ax.set_ylabel("Power Spectral Density", fontsize=14)
3405 ax.tick_params(axis="both", labelsize=12)
3407 # Set log-log scale
3408 ax.set_xscale("log")
3409 ax.set_yscale("log")
3411 # Add gridlines
3412 ax.grid(linestyle="--", color="grey", linewidth=1)
3413 # Ientify unique labels for legend
3414 handles = list(
3415 {
3416 label: handle
3417 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3418 }.values()
3419 )
3420 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3422 ax = plt.gca()
3423 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3425 # Save plot.
3426 _save_close_figure(fig, "histogram postage stamp", filename)
3429def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3430 cubes: iris.cube.Cube,
3431 coords: list[iris.coords.Coord],
3432 stamp_coordinate: str,
3433 filename: str,
3434 title: str,
3435 series_coordinate: str | None = None,
3436 **kwargs,
3437):
3438 """Plot and save power spectra for ensemble members in single plot.
3440 Parameters
3441 ----------
3442 cubes: Cube or CubeList
3443 Cube or Cubelist of the power spectrum data.
3444 coords: list[Coord]
3445 Coordinates to plot on the x-axis, one per cube.
3446 stamp_coordinate: str
3447 Coordinate that becomes different plots.
3448 filename: str
3449 Filename of the plot to write.
3450 title: str
3451 Plot title.
3452 series_coordinate: str, optional
3453 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3455 """
3456 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3457 model_colors_map = get_model_colors_map(cubes)
3459 line_marker = None
3460 line_width = 1
3462 # Compute ensemble statistics to show spread
3463 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3464 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3465 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3467 xcoord_global = mean_cube.coord(series_coordinate)
3468 x_global = xcoord_global.points
3470 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3471 xcoord = _select_series_coord(member, series_coordinate)
3472 xname = xcoord.points
3474 yfield = member.data # power spectrum
3475 color = "black"
3476 if model_colors_map: 3476 ↛ 3480line 3476 didn't jump to line 3480 because the condition on line 3476 was always true
3477 label = member.attributes.get("model_name") if i == 0 else None
3478 color = model_colors_map.get(label)
3480 if member.coord(stamp_coordinate).points == [0]:
3481 ax.plot(
3482 xname,
3483 yfield,
3484 color=color,
3485 marker=line_marker,
3486 ls="-",
3487 lw=line_width,
3488 label=f"{label} (control)"
3489 if len(member.coord(stamp_coordinate).points) > 1
3490 else label,
3491 )
3492 # Label with member number if part of an ensemble and not the control.
3493 else:
3494 ax.plot(
3495 xname,
3496 yfield,
3497 color=color,
3498 ls="-",
3499 lw=1.5,
3500 alpha=0.75,
3501 label=label,
3502 )
3504 # Set appropriate x-axis label based on coordinate
3505 if series_coordinate == "wavelength" or ( 3505 ↛ 3508line 3505 didn't jump to line 3508 because the condition on line 3505 was never true
3506 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3507 ):
3508 ax.set_xlabel("Wavelength (km)", fontsize=14)
3509 elif series_coordinate == "physical_wavenumber" or ( 3509 ↛ 3514line 3509 didn't jump to line 3514 because the condition on line 3509 was always true
3510 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3511 ):
3512 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3513 else: # frequency or check units
3514 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3515 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3516 else:
3517 ax.set_xlabel("Wavenumber", fontsize=14)
3519 # Add ensemble spread shading
3520 ax.fill_between(
3521 x_global,
3522 min_cube.data,
3523 max_cube.data,
3524 color="grey",
3525 alpha=0.3,
3526 label="Ensemble spread",
3527 )
3529 # Add ensemble mean line
3530 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3532 ax.set_ylabel("Power Spectral Density", fontsize=14)
3533 ax.tick_params(axis="both", labelsize=12)
3535 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3536 # Set log-log scale
3537 ax.set_xscale("log")
3538 ax.set_yscale("log")
3540 # Add gridlines
3541 ax.grid(linestyle="--", color="grey", linewidth=1)
3542 # Identify unique labels for legend
3543 handles = list(
3544 {
3545 label: handle
3546 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3547 }.values()
3548 )
3549 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3551 # Figure title.
3552 ax.set_title(title, fontsize=16)
3554 # Save plot.
3555 _save_close_figure(fig, "power spectra postage stamp", filename)