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