Coverage for src/CSET/operators/plot.py: 83%
1120 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-20 11:10 +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.nanmin(cube.data.filled(np.nan)):.3g} Max: {np.nanmax(cube.data.filled(np.nan)):.3g} Mean: {np.nanmean(cube.data.filled(np.nan)):.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.nanmin(cube_vec_mag.data.filled(np.nan)):.3g} Max: {np.nanmax(cube_vec_mag.data.filled(np.nan)):.3g} Mean: {np.nanmean(cube_vec_mag.data.filled(np.nan)):.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 try:
1619 ax.set_xlim(vmin, vmax)
1620 except ValueError:
1621 pass
1622 ax.tick_params(axis="both", labelsize=12)
1624 # Overlay grid-lines onto histogram plot.
1625 ax.grid(linestyle="--", color="grey", linewidth=1)
1626 if model_colors_map:
1627 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1629 # Save plot.
1630 _save_close_figure(fig, "histogram", filename)
1633def _plot_and_save_postage_stamp_histogram_series(
1634 cube: iris.cube.Cube,
1635 filename: str,
1636 title: str,
1637 stamp_coordinate: str,
1638 vmin: float,
1639 vmax: float,
1640 **kwargs,
1641):
1642 """Plot and save postage (ensemble members) stamps for a histogram series.
1644 Parameters
1645 ----------
1646 cube: Cube
1647 2 dimensional Cube of the data to plot as histogram.
1648 filename: str
1649 Filename of the plot to write.
1650 title: str
1651 Plot title.
1652 stamp_coordinate: str
1653 Coordinate that becomes different plots.
1654 vmin: float
1655 minimum for pdf x-axis
1656 vmax: float
1657 maximum for pdf x-axis
1658 """
1659 # Use the smallest square grid that will fit the members.
1660 nmember = len(cube.coord(stamp_coordinate).points)
1661 grid_rows = int(math.sqrt(nmember))
1662 grid_size = math.ceil(nmember / grid_rows)
1664 fig = plt.figure(
1665 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1666 )
1667 # Make a subplot for each member.
1668 for member, subplot in zip(
1669 cube.slices_over(stamp_coordinate),
1670 range(1, grid_size * grid_rows + 1),
1671 strict=False,
1672 ):
1673 # Implicit interface is much easier here, due to needing to have the
1674 # cartopy GeoAxes generated.
1675 plt.subplot(grid_rows, grid_size, subplot)
1676 # Reshape cube data into a single array to allow for a single histogram.
1677 # Otherwise we plot xdim histograms stacked.
1678 member_data_1d = (member.data).flatten()
1679 plt.hist(member_data_1d, density=True, stacked=True)
1680 axes = plt.gca()
1681 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1682 axes.set_title(f"{mtitle}")
1683 axes.set_xlim(vmin, vmax)
1685 # Overall figure title.
1686 fig.suptitle(title, fontsize=16)
1688 # Save plot.
1689 _save_close_figure(fig, "histogram postage stamp", filename)
1692def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1693 cube: iris.cube.Cube,
1694 filename: str,
1695 title: str,
1696 stamp_coordinate: str,
1697 vmin: float,
1698 vmax: float,
1699 **kwargs,
1700):
1701 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1702 ax.set_title(title, fontsize=16)
1703 ax.set_xlim(vmin, vmax)
1704 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1705 ax.set_ylabel("normalised probability density", fontsize=14)
1706 # Loop over all slices along the stamp_coordinate
1707 for member in cube.slices_over(stamp_coordinate):
1708 # Flatten the member data to 1D
1709 member_data_1d = member.data.flatten()
1710 # Plot the histogram using plt.hist
1711 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1712 plt.hist(
1713 member_data_1d,
1714 density=True,
1715 stacked=True,
1716 label=f"{mtitle}",
1717 )
1719 # Add a legend
1720 ax.legend(fontsize=16)
1722 # Save plot.
1723 _save_close_figure(fig, "histogram postage stamp", filename)
1726def _plot_and_save_scatter_series(
1727 cubes: iris.cube.Cube | iris.cube.CubeList,
1728 filename: str,
1729 title: str,
1730 vmin: float,
1731 vmax: float,
1732 hexbin: bool,
1733 **kwargs,
1734):
1735 """Plot and save a scatter plot series.
1737 Parameters
1738 ----------
1739 cubes: Cube or CubeList
1740 2 dimensional Cube or CubeList of the data to plot as scatter.
1741 filename: str
1742 Filename of the plot to write.
1743 title: str
1744 Plot title.
1745 vmin: float
1746 minimum for colorbar
1747 vmax: float
1748 maximum for colorbar
1749 hexbin: bool
1750 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1751 Else scatter of all points, with potential to overplot many comparisons on same plot.
1752 """
1753 if hexbin:
1754 # Check cubes using same functionality as the difference operator.
1755 if len(cubes) != 2:
1756 raise ValueError(
1757 "Cubes should contain exactly 2 cubes for hexbin plotting."
1758 )
1759 title = title.replace("scatter", "hexbin")
1760 filename = filename.replace("scatter", "hexbin")
1762 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1763 ax = plt.gca()
1765 model_colors_map = get_model_colors_map(cubes)
1767 percentiles = np.arange(0, 100, 5)
1768 percentiles[0] = 1
1769 percentiles[-1] = 99
1770 quantiles = iris.cube.CubeList()
1772 # Loop through all output cubes for both data points and overplotting quantiles.
1773 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1774 for plottype in ["points", "quantiles"]:
1775 nplot = 0
1776 for cube in iter_maybe(cubes):
1777 label = None
1778 color = "black"
1779 if model_colors_map: 1779 ↛ 1784line 1779 didn't jump to line 1784 because the condition on line 1779 was always true
1780 label = cube.attributes.get("model_name")
1781 color = model_colors_map[label]
1783 # Plot all data points
1784 if plottype == "points":
1785 if nplot > 0:
1786 if hexbin:
1787 hb = plt.hexbin(
1788 cubes[0].data.flatten(),
1789 cube.data.flatten(),
1790 alpha=0.3,
1791 gridsize=100,
1792 mincnt=1,
1793 )
1794 else:
1795 plt.scatter(
1796 cubes[0].data.flatten(),
1797 cube.data.flatten(),
1798 color=color,
1799 marker="+",
1800 label=None,
1801 alpha=0.3,
1802 )
1804 elif plottype == "quantiles": 1804 ↛ 1823line 1804 didn't jump to line 1823 because the condition on line 1804 was always true
1805 # Construct Q-Q plot
1806 quantiles.append(
1807 cube.collapsed(
1808 cube.coords(dim_coords=True),
1809 iris.analysis.PERCENTILE,
1810 percent=percentiles,
1811 )
1812 )
1813 if nplot > 0:
1814 iplt.scatter(
1815 quantiles[0],
1816 quantiles[-1],
1817 color=color,
1818 marker="o",
1819 label=label,
1820 edgecolors="black",
1821 )
1823 nplot = nplot + 1
1825 # Add some labels and tweak the style.
1826 ax.set_title(title, fontsize=16)
1827 ax.set_xlabel(
1828 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1829 )
1830 ax.set_ylabel(
1831 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1832 )
1833 ax.tick_params(axis="both", labelsize=12)
1834 ax.autoscale()
1836 # Set 1:1 line and equal axes if scatter plot of common cube names
1837 nameA = iter_maybe(cubes)[0].name()
1838 nameB = iter_maybe(cubes)[1].name()
1839 if any(part in nameB.split("_") for part in nameA.split("_")): 1839 ↛ 1850line 1839 didn't jump to line 1850 because the condition on line 1839 was always true
1840 lims = [
1841 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1842 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1843 ]
1844 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1845 ax.set_aspect("equal")
1846 ax.set_xlim(lims)
1847 ax.set_ylim(lims)
1849 # Overlay grid-lines onto scatter plot.
1850 ax.grid(linestyle="--", color="grey", linewidth=1)
1851 if model_colors_map: 1851 ↛ 1855line 1851 didn't jump to line 1855 because the condition on line 1851 was always true
1852 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1854 # Add colorbar if hexbin output
1855 if hexbin:
1856 cb = plt.colorbar(
1857 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1858 )
1859 cb.set_label("Number of data points", size=12)
1861 # Save plot.
1862 _save_close_figure(fig, "scatter", filename)
1865def _spatial_plot(
1866 method: Literal["contourf", "pcolormesh", "scatter"],
1867 cube: iris.cube.Cube,
1868 filename: str | None,
1869 sequence_coordinate: str,
1870 stamp_coordinate: str,
1871 overlay_cube: iris.cube.Cube | None = None,
1872 contour_cube: iris.cube.Cube | None = None,
1873 point_cube: iris.cube.Cube | None = None,
1874 **kwargs,
1875):
1876 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1878 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1879 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1880 is present then postage stamp plots will be produced.
1882 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1883 be overplotted on the same figure.
1885 Parameters
1886 ----------
1887 method: "contourf" | "pcolormesh" | "scatter"
1888 The plotting method to use.
1889 Select choice of "contourf" or "pcolormesh" for gridded data.
1890 Use "scatter" for point-based data.
1891 cube: Cube
1892 Iris cube of the data to plot. It should have two spatial dimensions,
1893 such as lat and lon, and may also have a another two dimension to be
1894 plotted sequentially and/or as postage stamp plots.
1895 filename: str | None
1896 Name of the plot to write, used as a prefix for plot sequences. If None
1897 uses the recipe name.
1898 sequence_coordinate: str
1899 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1900 This coordinate must exist in the cube.
1901 stamp_coordinate: str
1902 Coordinate about which to plot postage stamp plots. Defaults to
1903 ``"realization"``.
1904 overlay_cube: Cube | None, optional
1905 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1906 contour_cube: Cube | None, optional
1907 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1908 point_cube: Cube | None, optional
1909 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
1911 Raises
1912 ------
1913 ValueError
1914 If the cube doesn't have the right dimensions.
1915 TypeError
1916 If the cube isn't a single cube.
1917 """
1918 # Ensure we've got a single cube.
1919 cube = check_single_cube(cube)
1921 # Set title based on recipe metadata or use cube name
1922 recipe_title = get_recipe_metadata().get("title", cube.name())
1924 # Check if there is a valid stamp coordinate in cube dimensions.
1925 if stamp_coordinate == "realization": 1925 ↛ 1930line 1925 didn't jump to line 1930 because the condition on line 1925 was always true
1926 stamp_coordinate = check_stamp_coordinate(cube)
1928 # Make postage stamp plots if stamp_coordinate exists and has more than a
1929 # single point.
1930 plotting_func = _plot_and_save_spatial_plot
1931 try:
1932 if cube.coord(stamp_coordinate).shape[0] > 1:
1933 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1934 except iris.exceptions.CoordinateNotFoundError:
1935 pass
1937 # Produce a geographical scatter plot if the data have a
1938 # dimension called observation or model_obs_error
1939 if any(
1940 crd.var_name == "station"
1941 or crd.var_name == "Station_Name"
1942 or crd.var_name == "model_obs_error"
1943 for crd in cube.coords()
1944 ):
1945 plotting_func = _plot_and_save_spatial_plot
1946 method = "scatter"
1948 # Must have a sequence coordinate.
1949 try:
1950 cube.coord(sequence_coordinate)
1951 except iris.exceptions.CoordinateNotFoundError as err:
1952 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1954 # Create a plot for each value of the sequence coordinate.
1955 plot_index = []
1956 nplot = np.size(cube.coord(sequence_coordinate).points)
1958 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1959 # Set plot titles and filename
1960 seq_coord = cube_slice.coord(sequence_coordinate)
1962 if "model_name" in cube.attributes: 1962 ↛ 1963line 1962 didn't jump to line 1963 because the condition on line 1962 was never true
1963 model_name = cube.attributes["model_name"]
1964 else:
1965 model_name = None
1967 plot_title, plot_filename = _set_title_and_filename(
1968 seq_coord, nplot, recipe_title, filename, model_name=model_name
1969 )
1971 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1972 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1973 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1974 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1976 # Do the actual plotting.
1977 plotting_func(
1978 cube_slice,
1979 filename=plot_filename,
1980 stamp_coordinate=stamp_coordinate,
1981 title=plot_title,
1982 method=method,
1983 overlay_cube=overlay_slice,
1984 contour_cube=contour_slice,
1985 point_cube=point_slice,
1986 **kwargs,
1987 )
1988 plot_index.append(plot_filename)
1990 # Add list of plots to plot metadata.
1991 complete_plot_index = _append_to_plot_index(plot_index)
1993 # Make a page to display the plots.
1994 _make_plot_html_page(complete_plot_index)
1997####################
1998# Public functions #
1999####################
2002def spatial_contour_plot(
2003 cube: iris.cube.Cube,
2004 filename: str | None = None,
2005 sequence_coordinate: str = "time",
2006 stamp_coordinate: str = "realization",
2007 **kwargs,
2008) -> iris.cube.Cube:
2009 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2011 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2012 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2013 is present then postage stamp plots will be produced.
2015 Parameters
2016 ----------
2017 cube: Cube
2018 Iris cube of the data to plot. It should have two spatial dimensions,
2019 such as lat and lon, and may also have a another two dimension to be
2020 plotted sequentially and/or as postage stamp plots.
2021 filename: str, optional
2022 Name of the plot to write, used as a prefix for plot sequences. Defaults
2023 to the recipe name.
2024 sequence_coordinate: str, optional
2025 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2026 This coordinate must exist in the cube.
2027 stamp_coordinate: str, optional
2028 Coordinate about which to plot postage stamp plots. Defaults to
2029 ``"realization"``.
2031 Returns
2032 -------
2033 Cube
2034 The original cube (so further operations can be applied).
2036 Raises
2037 ------
2038 ValueError
2039 If the cube doesn't have the right dimensions.
2040 TypeError
2041 If the cube isn't a single cube.
2042 """
2043 _spatial_plot(
2044 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2045 )
2046 return cube
2049def spatial_pcolormesh_plot(
2050 cubes: iris.cube.Cube | iris.cube.CubeList,
2051 filename: str | None = None,
2052 sequence_coordinate: str = "time",
2053 stamp_coordinate: str = "realization",
2054 **kwargs,
2055) -> iris.cube.Cube:
2056 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2058 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2059 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2060 is present then postage stamp plots will be produced.
2062 This function is significantly faster than ``spatial_contour_plot``,
2063 especially at high resolutions, and should be preferred unless contiguous
2064 contour areas are important.
2066 Parameters
2067 ----------
2068 cube: Cubes
2069 Iris cube or cubelist of the data to plot. Each cube should have two spatial dimensions,
2070 such as lat and lon, and may also have a another two dimension to be
2071 plotted sequentially and/or as postage stamp plots.
2072 filename: str, optional
2073 Name of the plot to write, used as a prefix for plot sequences. Defaults
2074 to the recipe name.
2075 sequence_coordinate: str, optional
2076 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2077 This coordinate must exist in the cube.
2078 stamp_coordinate: str, optional
2079 Coordinate about which to plot postage stamp plots. Defaults to
2080 ``"realization"``.
2082 Returns
2083 -------
2084 Cubes
2085 The original cube/cubelist (so further operations can be applied).
2087 Raises
2088 ------
2089 ValueError
2090 If the cube doesn't have the right dimensions.
2091 """
2092 if isinstance(cubes, iris.cube.CubeList): 2092 ↛ 2093line 2092 didn't jump to line 2093 because the condition on line 2092 was never true
2093 for model_cube in cubes:
2094 _spatial_plot(
2095 "pcolormesh",
2096 model_cube,
2097 filename,
2098 sequence_coordinate,
2099 stamp_coordinate,
2100 **kwargs,
2101 )
2102 elif isinstance(cubes, iris.cube.Cube): 2102 ↛ 2111line 2102 didn't jump to line 2111 because the condition on line 2102 was always true
2103 _spatial_plot(
2104 "pcolormesh",
2105 cubes,
2106 filename,
2107 sequence_coordinate,
2108 stamp_coordinate,
2109 **kwargs,
2110 )
2111 return cubes
2114def spatial_multi_pcolormesh_plot(
2115 cube: iris.cube.Cube,
2116 overlay_cube: iris.cube.Cube | None = None,
2117 contour_cube: iris.cube.Cube | None = None,
2118 point_cube: iris.cube.Cube | None = None,
2119 filename: str | None = None,
2120 sequence_coordinate: str = "time",
2121 stamp_coordinate: str = "realization",
2122 **kwargs,
2123) -> iris.cube.Cube:
2124 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2126 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2127 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2128 is present then postage stamp plots will be produced.
2130 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2132 If specified, contours of a contour_cube can be overplotted on top of those.
2134 If specified, a spatial scatter map of point_cube can be overplotted.
2136 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2138 This function is significantly faster than ``spatial_contour_plot``,
2139 especially at high resolutions, and should be preferred unless contiguous
2140 contour areas are important.
2142 Parameters
2143 ----------
2144 cube: Cube
2145 Iris cube of the data to plot. It should have two spatial dimensions,
2146 such as lat and lon, and may also have two additional dimensions to be
2147 plotted sequentially and/or as postage stamp plots.
2148 overlay_cube: Cube, optional
2149 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2150 such as lat and lon, and may also have two additional dimensions to be
2151 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.
2152 If not provided, output plot generated without overlay cube.
2153 contour_cube: Cube, optional
2154 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,
2155 such as lat and lon, and may also have two additional dimensions to be
2156 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2157 point_cube: Cube, optional
2158 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
2159 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2160 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2161 postage stamp plots. If not provided, output plot generated without point-based layer.
2162 filename: str, optional
2163 Name of the plot to write, used as a prefix for plot sequences. Defaults
2164 to the recipe name.
2165 sequence_coordinate: str, optional
2166 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2167 This coordinate must exist in the cube.
2168 stamp_coordinate: str, optional
2169 Coordinate about which to plot postage stamp plots. Defaults to
2170 ``"realization"``.
2172 Returns
2173 -------
2174 Cube
2175 The original cube (so further operations can be applied).
2177 Raises
2178 ------
2179 ValueError
2180 If the cube doesn't have the right dimensions.
2181 TypeError
2182 If the cube isn't a single cube.
2183 """
2184 _spatial_plot(
2185 "pcolormesh",
2186 cube,
2187 filename,
2188 sequence_coordinate,
2189 stamp_coordinate,
2190 overlay_cube=overlay_cube,
2191 contour_cube=contour_cube,
2192 point_cube=point_cube,
2193 )
2194 return cube, overlay_cube, contour_cube, point_cube
2197# TODO: Expand function to handle ensemble data.
2198# line_coordinate: str, optional
2199# Coordinate about which to plot multiple lines. Defaults to
2200# ``"realization"``.
2201def plot_line_series(
2202 cube: iris.cube.Cube | iris.cube.CubeList,
2203 filename: str | None = None,
2204 series_coordinate: str = "time",
2205 sequence_coordinate: str = "time",
2206 # add the following for ensembles
2207 stamp_coordinate: str = "realization",
2208 single_plot: bool = False,
2209 **kwargs,
2210) -> iris.cube.Cube | iris.cube.CubeList:
2211 """Plot a line plot for the specified coordinate.
2213 The Cube or CubeList must be 1D.
2215 Parameters
2216 ----------
2217 iris.cube | iris.cube.CubeList
2218 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2219 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2220 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2221 filename: str, optional
2222 Name of the plot to write, used as a prefix for plot sequences. Defaults
2223 to the recipe name.
2224 series_coordinate: str, optional
2225 Coordinate about which to make a series. Defaults to ``"time"``. This
2226 coordinate must exist in the cube.
2228 Returns
2229 -------
2230 iris.cube.Cube | iris.cube.CubeList
2231 The original Cube or CubeList (so further operations can be applied).
2233 Raises
2234 ------
2235 ValueError
2236 If the cubes don't have the right dimensions.
2237 TypeError
2238 If the cube isn't a Cube or CubeList.
2239 """
2240 # Ensure we have a name for the plot file.
2241 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2243 num_models = get_num_models(cube)
2245 validate_cube_shape(cube, num_models)
2247 # Iterate over all cubes and extract coordinate to plot.
2248 cubes = iris.cube.CubeList(iter_maybe(cube))
2249 coords = []
2250 for model_cube in cubes:
2251 try:
2252 coords.append(model_cube.coord(series_coordinate))
2253 except iris.exceptions.CoordinateNotFoundError as err:
2254 raise ValueError(
2255 f"Cube must have a {series_coordinate} coordinate."
2256 ) from err
2257 # Count dimensions excluding realization
2258 ndim = model_cube.ndim
2260 if model_cube.coords("realization"): 2260 ↛ 2267line 2260 didn't jump to line 2267 because the condition on line 2260 was always true
2261 realization_dims = model_cube.coord_dims("realization")
2263 # Only subtract if realization is a dimension coordinate
2264 if realization_dims:
2265 ndim -= len(realization_dims)
2267 if ndim > 2:
2268 raise ValueError(
2269 "Cube must be 1D or 2D (excluding any realization dimension)."
2270 )
2272 plot_index = []
2274 # Check if this is a spectral plot by looking for spectral coordinates
2275 is_spectral_plot = series_coordinate in [
2276 "frequency",
2277 "physical_wavenumber",
2278 "wavelength",
2279 ]
2281 if is_spectral_plot:
2282 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2283 # coordinate frequency/wavenumber.
2284 # If several power spectra are plotted with time as sequence_coordinate for the
2285 # time slider option.
2287 # Internal plotting function.
2288 plotting_func = _plot_and_save_line_power_spectrum_series
2290 for model_cube in cubes:
2291 try:
2292 model_cube.coord(sequence_coordinate)
2293 except iris.exceptions.CoordinateNotFoundError as err:
2294 raise ValueError(
2295 f"Cube must have a {sequence_coordinate} coordinate."
2296 ) from err
2298 if num_models == 1: 2298 ↛ 2313line 2298 didn't jump to line 2313 because the condition on line 2298 was always true
2299 # check for ensembles
2300 if ( 2300 ↛ 2304line 2300 didn't jump to line 2304 because the condition on line 2300 was never true
2301 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2302 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2303 ):
2304 if single_plot:
2305 # Plot spectra, mean and ensemble spread on 1 plot
2306 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2307 else:
2308 # Plot postage stamps
2309 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2310 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2311 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2312 else:
2313 all_points = sorted(
2314 set(
2315 itertools.chain.from_iterable(
2316 cb.coord(sequence_coordinate).points for cb in cubes
2317 )
2318 )
2319 )
2320 all_slices = list(
2321 itertools.chain.from_iterable(
2322 cb.slices_over(sequence_coordinate) for cb in cubes
2323 )
2324 )
2325 # Matched slices (matched by seq coord point; it may happen that
2326 # evaluated models do not cover the same seq coord range, hence matching
2327 # necessary)
2328 cube_iterables = [
2329 iris.cube.CubeList(
2330 s
2331 for s in all_slices
2332 if s.coord(sequence_coordinate).points[0] == point
2333 )
2334 for point in all_points
2335 ]
2336 nplot = len(all_points)
2338 # Create a plot for each value of the sequence coordinate. Allowing for
2339 # multiple cubes in a CubeList to be plotted in the same plot for similar
2340 # sequence values. Passing a CubeList into the internal plotting function
2341 # for similar values of the sequence coordinate. cube_slice can be an
2342 # iris.cube.Cube or an iris.cube.CubeList.
2344 for cube_slice in cube_iterables:
2345 # Normalize cube_slice to a list of cubes
2346 if isinstance(cube_slice, iris.cube.CubeList): 2346 ↛ 2347line 2346 didn't jump to line 2347 because the condition on line 2346 was never true
2347 cubes = list(cube_slice)
2348 elif isinstance(cube_slice, iris.cube.Cube): 2348 ↛ 2351line 2348 didn't jump to line 2351 because the condition on line 2348 was always true
2349 cubes = [cube_slice]
2350 else:
2351 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2353 # Use sequence value so multiple sequences can merge.
2354 seq_coord = cube_slice[0].coord(sequence_coordinate)
2355 plot_title, plot_filename = _set_title_and_filename(
2356 seq_coord, nplot, recipe_title, filename
2357 )
2359 # Format the coordinate value in a unit appropriate way.
2360 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2362 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2363 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2363 ↛ 2364line 2363 didn't jump to line 2364 because the condition on line 2363 was never true
2364 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2366 # Do the actual plotting.
2367 plotting_func(
2368 cube_slice,
2369 coords,
2370 stamp_coordinate,
2371 plot_filename,
2372 title,
2373 series_coordinate,
2374 )
2376 plot_index.append(plot_filename)
2377 else:
2378 # Format the title and filename using plotted series coordinate
2379 nplot = 1
2380 seq_coord = coords[0]
2381 plot_title, plot_filename = _set_title_and_filename(
2382 seq_coord, nplot, recipe_title, filename
2383 )
2385 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2386 if (
2387 "station" in [c.name() for c in cubes[0].coords()]
2388 and len(cubes[0].coord("station").points) > 1
2389 ):
2390 for station in cubes[0].coord("station").points:
2391 station_cubes = cubes.extract(iris.Constraint(station=station))
2392 station_name = station_cubes[0].coord("Station_Name").points[0]
2393 station_plotname = plot_filename.replace(
2394 ".png", "_" + station_name + ".png"
2395 )
2396 _plot_and_save_line_series(
2397 station_cubes,
2398 coords,
2399 "realization",
2400 station_plotname,
2401 f"{plot_title} {station_name}",
2402 )
2403 plot_index.append(station_plotname)
2405 else:
2406 # Do the actual plotting for all other series coordinate options.
2407 _plot_and_save_line_series(
2408 cubes, coords, stamp_coordinate, plot_filename, plot_title
2409 )
2411 plot_index.append(plot_filename)
2413 # append plot to list of plots
2414 complete_plot_index = _append_to_plot_index(plot_index)
2416 # Make a page to display the plots.
2417 _make_plot_html_page(complete_plot_index)
2419 return cube
2422def plot_vertical_line_series(
2423 cubes: iris.cube.Cube | iris.cube.CubeList,
2424 filename: str | None = None,
2425 series_coordinate: str = "model_level_number",
2426 sequence_coordinate: str = "time",
2427 # line_coordinate: str = "realization",
2428 **kwargs,
2429) -> iris.cube.Cube | iris.cube.CubeList:
2430 """Plot a line plot against a type of vertical coordinate.
2432 The Cube or CubeList must be 1D.
2434 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2435 then a sequence of plots will be produced.
2437 Parameters
2438 ----------
2439 iris.cube | iris.cube.CubeList
2440 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2441 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2442 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2443 filename: str, optional
2444 Name of the plot to write, used as a prefix for plot sequences. Defaults
2445 to the recipe name.
2446 series_coordinate: str, optional
2447 Coordinate to plot on the y-axis. Can be ``pressure`` or
2448 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2449 for LFRic. Defaults to ``model_level_number``.
2450 This coordinate must exist in the cube.
2451 sequence_coordinate: str, optional
2452 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2453 This coordinate must exist in the cube.
2455 Returns
2456 -------
2457 iris.cube.Cube | iris.cube.CubeList
2458 The original Cube or CubeList (so further operations can be applied).
2459 Plotted data.
2461 Raises
2462 ------
2463 ValueError
2464 If the cubes doesn't have the right dimensions.
2465 TypeError
2466 If the cube isn't a Cube or CubeList.
2467 """
2468 # Ensure we have a name for the plot file.
2469 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2471 cubes = iter_maybe(cubes)
2472 # Initialise empty list to hold all data from all cubes in a CubeList
2473 all_data = []
2475 # Store min/max ranges for x range.
2476 x_levels = []
2478 num_models = get_num_models(cubes)
2480 validate_cube_shape(cubes, num_models)
2482 # Iterate over all cubes in cube or CubeList and plot.
2483 coords = []
2484 for cube in cubes:
2485 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2486 try:
2487 coords.append(cube.coord(series_coordinate))
2488 except iris.exceptions.CoordinateNotFoundError as err:
2489 raise ValueError(
2490 f"Cube must have a {series_coordinate} coordinate."
2491 ) from err
2493 try:
2494 if cube.ndim > 1 or not cube.coords("realization"): 2494 ↛ 2502line 2494 didn't jump to line 2502 because the condition on line 2494 was always true
2495 cube.coord(sequence_coordinate)
2496 except iris.exceptions.CoordinateNotFoundError as err:
2497 raise ValueError(
2498 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2499 ) from err
2501 # Get minimum and maximum from levels information.
2502 _, levels, _ = colorbar_map_levels(cube, axis="x")
2503 if levels is not None: 2503 ↛ 2507line 2503 didn't jump to line 2507 because the condition on line 2503 was always true
2504 x_levels.append(min(levels))
2505 x_levels.append(max(levels))
2506 else:
2507 all_data.append(cube.data)
2509 if len(x_levels) == 0: 2509 ↛ 2511line 2509 didn't jump to line 2511 because the condition on line 2509 was never true
2510 # Combine all data into a single NumPy array
2511 combined_data = np.concatenate(all_data)
2513 # Set the lower and upper limit for the x-axis to ensure all plots have
2514 # same range. This needs to read the whole cube over the range of the
2515 # sequence and if applicable postage stamp coordinate.
2516 vmin = np.floor(combined_data.min())
2517 vmax = np.ceil(combined_data.max())
2518 else:
2519 vmin = min(x_levels)
2520 vmax = max(x_levels)
2522 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2523 # a single profile directly without iterating over a sequence.
2524 sequence_coords = [
2525 cube.coord(sequence_coordinate)
2526 for cube in cubes
2527 if cube.coords(sequence_coordinate)
2528 ]
2529 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2530 np.size(coord.points) > 1 for coord in sequence_coords
2531 )
2532 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2533 np.size(coord.points) == 1 for coord in sequence_coords
2534 )
2536 plot_index = []
2537 if has_sequence_coord: 2537 ↛ 2562line 2537 didn't jump to line 2562 because the condition on line 2537 was always true
2538 # Matching the slices (matching by seq coord point; it may happen that
2539 # evaluated models do not cover the same seq coord range, hence matching
2540 # necessary)
2541 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2542 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2543 for cubes_slice in cube_iterables:
2544 # Format the coordinate value in a unit appropriate way.
2545 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2546 plot_title, plot_filename = _set_title_and_filename(
2547 seq_coord, nplot, recipe_title, filename
2548 )
2550 # Do the actual plotting.
2551 _plot_and_save_vertical_line_series(
2552 cubes_slice,
2553 coords,
2554 "realization",
2555 plot_filename,
2556 series_coordinate,
2557 title=plot_title,
2558 vmin=vmin,
2559 vmax=vmax,
2560 )
2561 plot_index.append(plot_filename)
2562 elif has_scalar_sequence_coord:
2563 # Scalar sequence coordinate (typically aggregated time bounds):
2564 # make one plot and include sequence period in title/filename.
2565 plot_title, plot_filename = _set_title_and_filename(
2566 sequence_coords[0], 1, recipe_title, filename
2567 )
2569 _plot_and_save_vertical_line_series(
2570 cubes,
2571 coords,
2572 "realization",
2573 plot_filename,
2574 series_coordinate,
2575 title=plot_title,
2576 vmin=vmin,
2577 vmax=vmax,
2578 )
2579 plot_index.append(plot_filename)
2580 else:
2581 # 1D case: no sequence coordinate, plot a single profile.
2582 plot_title = recipe_title
2583 if filename:
2584 plot_filename = filename
2585 else:
2586 plot_filename = f"{slugify(plot_title)}.png"
2588 _plot_and_save_vertical_line_series(
2589 cubes,
2590 coords,
2591 "realization",
2592 plot_filename,
2593 series_coordinate,
2594 title=plot_title,
2595 vmin=vmin,
2596 vmax=vmax,
2597 )
2598 plot_index.append(plot_filename)
2600 # Add list of plots to plot metadata.
2601 complete_plot_index = _append_to_plot_index(plot_index)
2603 # Make a page to display the plots.
2604 _make_plot_html_page(complete_plot_index)
2606 return cubes
2609def qq_plot(
2610 cubes: iris.cube.CubeList,
2611 coordinates: list[str],
2612 percentiles: list[float],
2613 model_names: list[str],
2614 filename: str | None = None,
2615 one_to_one: bool = True,
2616 **kwargs,
2617) -> iris.cube.CubeList:
2618 """Plot a Quantile-Quantile plot between two models for common time points.
2620 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2621 collapsed within the operator over all specified coordinates such as
2622 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2623 ensemble members to ensure a 1D cube (array).
2625 Parameters
2626 ----------
2627 cubes: iris.cube.CubeList
2628 Two cubes of the same variable with different models.
2629 coordinate: list[str]
2630 The list of coordinates to collapse over. This list should be
2631 every coordinate within the cube to result in a 1D cube around
2632 the percentile coordinate.
2633 percent: list[float]
2634 A list of percentiles to appear in the plot.
2635 model_names: list[str]
2636 A list of model names to appear on the axis of the plot.
2637 filename: str, optional
2638 Filename of the plot to write.
2639 one_to_one: bool, optional
2640 If True a 1:1 line is plotted; if False it is not. Default is True.
2642 Raises
2643 ------
2644 ValueError
2645 When the cubes are not compatible.
2647 Notes
2648 -----
2649 The quantile-quantile plot is a variant on the scatter plot representing
2650 two datasets by their quantiles (percentiles) for common time points.
2651 This plot does not use a theoretical distribution to compare against, but
2652 compares percentiles of two datasets. This plot does
2653 not use all raw data points, but plots the selected percentiles (quantiles) of
2654 each variable instead for the two datasets, thereby normalising the data for a
2655 direct comparison between the selected percentiles of the two dataset distributions.
2657 Quantile-quantile plots are valuable for comparing against
2658 observations and other models. Identical percentiles between the variables
2659 will lie on the one-to-one line implying the values correspond well to each
2660 other. Where there is a deviation from the one-to-one line a range of
2661 possibilities exist depending on how and where the data is shifted (e.g.,
2662 Wilks 2011 [Wilks2011]_).
2664 For distributions above the one-to-one line the distribution is left-skewed;
2665 below is right-skewed. A distinct break implies a bimodal distribution, and
2666 closer values/values further apart at the tails imply poor representation of
2667 the extremes.
2669 """
2670 # Check cubes using same functionality as the difference operator.
2671 if len(cubes) != 2:
2672 raise ValueError("cubes should contain exactly 2 cubes.")
2673 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2674 other: Cube = cubes.extract_cube(
2675 iris.Constraint(
2676 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2677 )
2678 )
2680 # Get spatial coord names.
2681 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2682 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2684 # Ensure cubes to compare are on common differencing grid.
2685 # This is triggered if either
2686 # i) latitude and longitude shapes are not the same. Note grid points
2687 # are not compared directly as these can differ through rounding
2688 # errors.
2689 # ii) or variables are known to often sit on different grid staggering
2690 # in different models (e.g. cell center vs cell edge), as is the case
2691 # for UM and LFRic comparisons.
2692 # In future greater choice of regridding method might be applied depending
2693 # on variable type. Linear regridding can in general be appropriate for smooth
2694 # variables. Care should be taken with interpretation of differences
2695 # given this dependency on regridding.
2696 if (
2697 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2698 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2699 ) or (
2700 base.long_name
2701 in [
2702 "eastward_wind_at_10m",
2703 "northward_wind_at_10m",
2704 "northward_wind_at_cell_centres",
2705 "eastward_wind_at_cell_centres",
2706 "zonal_wind_at_pressure_levels",
2707 "meridional_wind_at_pressure_levels",
2708 "potential_vorticity_at_pressure_levels",
2709 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2710 ]
2711 ):
2712 logger.debug("Linear regridding base cube to other grid to compute differences")
2713 base = regrid_onto_cube(base, other, method="Linear")
2715 # Extract just common time points.
2716 base, other = _extract_common_time_points(base, other)
2718 # Equalise attributes so we can merge.
2719 fully_equalise_attributes([base, other])
2720 logger.debug("Base: %s\nOther: %s", base, other)
2722 # Collapse cubes.
2723 base = collapse(
2724 base,
2725 coordinate=coordinates,
2726 method="PERCENTILE",
2727 additional_percent=percentiles,
2728 )
2729 other = collapse(
2730 other,
2731 coordinate=coordinates,
2732 method="PERCENTILE",
2733 additional_percent=percentiles,
2734 )
2736 # Ensure we have a name for the plot file.
2737 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2738 title = f"{recipe_title}"
2740 if filename is None:
2741 filename = slugify(recipe_title)
2743 # Add file extension.
2744 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2746 # Do the actual plotting on a scatter plot
2747 _plot_and_save_scatter_plot(
2748 base, other, plot_filename, title, one_to_one, model_names
2749 )
2751 # Add list of plots to plot metadata.
2752 plot_index = _append_to_plot_index([plot_filename])
2754 # Make a page to display the plots.
2755 _make_plot_html_page(plot_index)
2757 return iris.cube.CubeList([base, other])
2760def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2761 """
2762 Plot a Hinton style triangle/scorecard plot.
2764 This plot type can be useful for summarising high level information, such as comparing
2765 how 'skillful' two models are when verified against observations for a variety of metrics,
2766 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2767 than customisable by the user as input arguments; many have been designed to automatically
2768 scale the plot depending on the number of x and y components.
2770 Parameters
2771 ----------
2772 change: np.ndarray
2773 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2774 size/direction.
2775 signif: np.ndarray
2776 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2777 xaxis_labels: list
2778 List of labels for the xaxis (must match the second dimension length of signif and change,
2779 along with magnitude if not None).
2780 yaxis_labels: list
2781 List of labels for the yaxis (must match the first dimension length of signif and change,
2782 along with magnitude if not None).
2783 magnitude: np.ndarray | None
2784 Optional 2D array, matching the shape of change, signif, which contains numerical values
2785 the user wishes to display under each respective triangle.
2787 Returns
2788 -------
2789 matplotlib axes object to either display or do further modifications to.
2790 """
2791 # Setup colors of triangles
2792 color_pos = "#7CAE00"
2793 color_neg = "#7B68EE"
2795 # Setup cell/text size ratios
2796 figsize = None
2797 cell_size_in = 0.35
2798 text_row_ratio = 0.25
2800 # Ensure arrays, and change to bool for sig.
2801 change = np.asarray(change)
2802 signif = np.asarray(signif).astype(bool)
2803 if magnitude is not None: 2803 ↛ 2804line 2803 didn't jump to line 2804 because the condition on line 2803 was never true
2804 magnitude = np.asarray(magnitude)
2806 # Get the number of x and y elements
2807 ny, nx = change.shape
2809 # Build non-uniform y coordinates
2810 tri_height = 1.0
2811 txt_height = text_row_ratio
2813 tri_y = []
2814 txt_y = []
2815 y_edges = [0.0]
2817 y = 0.0
2818 for _j in range(ny):
2819 tri_y.append(y + tri_height / 2)
2820 y += tri_height
2821 y_edges.append(y)
2823 if magnitude is not None: 2823 ↛ 2824line 2823 didn't jump to line 2824 because the condition on line 2823 was never true
2824 txt_y.append(y + txt_height / 2)
2825 y += txt_height
2826 y_edges.append(y)
2828 total_height = y
2830 # Dynamic figure size
2831 if figsize is None: 2831 ↛ 2836line 2831 didn't jump to line 2836 because the condition on line 2831 was always true
2832 width = nx * cell_size_in
2833 height = total_height * cell_size_in + 2
2834 figsize = (width, height)
2836 fig, ax = plt.subplots(figsize=figsize)
2838 # Setup axes and grid.
2839 ax.set_aspect("equal", adjustable="box")
2840 ax.set_xlim(-0.5, nx - 0.5)
2841 ax.set_ylim(0, total_height)
2843 ax.set_xticks(np.arange(nx))
2844 ax.set_xticklabels(xaxis_labels, rotation=90)
2846 ax.set_yticks(tri_y)
2847 ax.set_yticklabels(yaxis_labels)
2849 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2850 ax.set_yticks(y_edges, minor=True)
2852 ax.set_axisbelow(True)
2853 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2854 ax.grid(False, which="major")
2855 ax.tick_params(which="minor", length=0)
2857 ax.invert_yaxis()
2859 # Compute marker scaling (fixed overlap)
2860 fig.canvas.draw()
2862 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2863 width_in, height_in = bbox.width, bbox.height
2865 cell_w = (width_in * fig.dpi) / nx
2866 cell_h = (height_in * fig.dpi) / total_height
2867 cell_pixels = min(cell_w, cell_h)
2869 max_marker_size = (0.6 * cell_pixels) ** 2
2871 text_fontsize = cell_pixels * 0.15
2873 # Plot triangles + text
2874 for j in range(ny):
2875 for i in range(nx):
2876 val = change[j, i]
2877 if np.isnan(val): 2877 ↛ 2878line 2877 didn't jump to line 2878 because the condition on line 2877 was never true
2878 continue
2880 if abs(val) < 0.01: 2880 ↛ 2881line 2880 didn't jump to line 2881 because the condition on line 2880 was never true
2881 continue
2883 sig = signif[j, i]
2884 size = max_marker_size * abs(val)
2886 # Triangle style
2887 if val >= 0:
2888 marker = "^"
2889 color = color_pos
2890 else:
2891 marker = "v"
2892 color = color_neg
2894 if sig:
2895 edgecolor = "black"
2896 linewidth = 0.6
2897 else:
2898 edgecolor = "none"
2899 linewidth = 0.0
2901 # Triangle
2902 ax.scatter(
2903 i,
2904 tri_y[j],
2905 s=size,
2906 marker=marker,
2907 c=color,
2908 edgecolors=edgecolor,
2909 linewidths=linewidth,
2910 zorder=3,
2911 clip_on=True, # ensures no rendering bleed
2912 )
2914 # Text row
2915 if magnitude is not None: 2915 ↛ 2916line 2915 didn't jump to line 2916 because the condition on line 2915 was never true
2916 mag_val = magnitude[j, i]
2918 if not np.isnan(mag_val):
2919 ax.text(
2920 i,
2921 txt_y[j],
2922 f"{mag_val:.1f}",
2923 ha="center",
2924 va="center",
2925 fontsize=text_fontsize,
2926 color="black",
2927 zorder=4,
2928 )
2930 plt.tight_layout()
2931 return fig, ax
2934def scatter_plot(
2935 cube_x: iris.cube.Cube | iris.cube.CubeList,
2936 cube_y: iris.cube.Cube | iris.cube.CubeList,
2937 filename: str | None = None,
2938 one_to_one: bool = True,
2939 **kwargs,
2940) -> iris.cube.CubeList:
2941 """Plot a scatter plot between two variables.
2943 Both cubes must be 1D.
2945 Parameters
2946 ----------
2947 cube_x: Cube | CubeList
2948 1 dimensional Cube of the data to plot on y-axis.
2949 cube_y: Cube | CubeList
2950 1 dimensional Cube of the data to plot on x-axis.
2951 filename: str, optional
2952 Filename of the plot to write.
2953 one_to_one: bool, optional
2954 If True a 1:1 line is plotted; if False it is not. Default is True.
2956 Returns
2957 -------
2958 cubes: CubeList
2959 CubeList of the original x and y cubes for further processing.
2961 Raises
2962 ------
2963 ValueError
2964 If the cube doesn't have the right dimensions and cubes not the same
2965 size.
2966 TypeError
2967 If the cube isn't a single cube.
2969 Notes
2970 -----
2971 Scatter plots are used for determining if there is a relationship between
2972 two variables. Positive relations have a slope going from bottom left to top
2973 right; Negative relations have a slope going from top left to bottom right.
2974 """
2975 # Iterate over all cubes in cube or CubeList and plot.
2976 for cube_iter in iter_maybe(cube_x):
2977 # Check cubes are correct shape.
2978 cube_iter = check_single_cube(cube_iter)
2979 if cube_iter.ndim > 1:
2980 raise ValueError("cube_x must be 1D.")
2982 # Iterate over all cubes in cube or CubeList and plot.
2983 for cube_iter in iter_maybe(cube_y):
2984 # Check cubes are correct shape.
2985 cube_iter = check_single_cube(cube_iter)
2986 if cube_iter.ndim > 1:
2987 raise ValueError("cube_y must be 1D.")
2989 # Ensure we have a name for the plot file.
2990 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2991 title = f"{recipe_title}"
2993 if filename is None:
2994 filename = slugify(recipe_title)
2996 # Add file extension.
2997 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2999 # Do the actual plotting.
3000 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3002 # Add list of plots to plot metadata.
3003 plot_index = _append_to_plot_index([plot_filename])
3005 # Make a page to display the plots.
3006 _make_plot_html_page(plot_index)
3008 return iris.cube.CubeList([cube_x, cube_y])
3011def vector_plot(
3012 cube_u: iris.cube.Cube,
3013 cube_v: iris.cube.Cube,
3014 filename: str | None = None,
3015 sequence_coordinate: str = "time",
3016 **kwargs,
3017) -> iris.cube.CubeList:
3018 """Plot a vector plot based on the input u and v components."""
3019 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
3021 # Cubes must have a matching sequence coordinate.
3022 try:
3023 # Check that the u and v cubes have the same sequence coordinate.
3024 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3024 ↛ anywhereline 3024 didn't jump anywhere: it always raised an exception.
3025 raise ValueError("Coordinates do not match.")
3026 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3027 raise ValueError(
3028 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3029 ) from err
3031 # Create a plot for each value of the sequence coordinate.
3032 plot_index = []
3033 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3034 for cube_u_slice, cube_v_slice in zip(
3035 cube_u.slices_over(sequence_coordinate),
3036 cube_v.slices_over(sequence_coordinate),
3037 strict=True,
3038 ):
3039 # Format the coordinate value in a unit appropriate way.
3040 seq_coord = cube_u_slice.coord(sequence_coordinate)
3041 plot_title, plot_filename = _set_title_and_filename(
3042 seq_coord, nplot, recipe_title, filename
3043 )
3045 # Do the actual plotting.
3046 _plot_and_save_vector_plot(
3047 cube_u_slice,
3048 cube_v_slice,
3049 filename=plot_filename,
3050 title=plot_title,
3051 method="pcolormesh",
3052 )
3053 plot_index.append(plot_filename)
3055 # Add list of plots to plot metadata.
3056 complete_plot_index = _append_to_plot_index(plot_index)
3058 # Make a page to display the plots.
3059 _make_plot_html_page(complete_plot_index)
3061 return iris.cube.CubeList([cube_u, cube_v])
3064def plot_histogram_series(
3065 cubes: iris.cube.Cube | iris.cube.CubeList,
3066 filename: str | None = None,
3067 sequence_coordinate: str = "time",
3068 stamp_coordinate: str = "realization",
3069 single_plot: bool = False,
3070 **kwargs,
3071) -> iris.cube.Cube | iris.cube.CubeList:
3072 """Plot a histogram plot for each vertical level provided.
3074 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3075 is present then a sequence of plots will be produced using the time slider
3076 functionality to scroll through histograms against time. If a
3077 stamp_coordinate is present then postage stamp plots will be produced. If
3078 stamp_coordinate and single_plot is True, all postage stamp plots will be
3079 plotted in a single plot instead of separate postage stamp plots.
3081 Parameters
3082 ----------
3083 cubes: Cube | iris.cube.CubeList
3084 Iris cube or CubeList of the data to plot. It should have a single dimension other
3085 than the stamp coordinate.
3086 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3087 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3088 filename: str, optional
3089 Name of the plot to write, used as a prefix for plot sequences. Defaults
3090 to the recipe name.
3091 sequence_coordinate: str, optional
3092 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3093 This coordinate must exist in the cube and will be used for the time
3094 slider.
3095 stamp_coordinate: str, optional
3096 Coordinate about which to plot postage stamp plots. Defaults to
3097 ``"realization"``.
3098 single_plot: bool, optional
3099 If True, all postage stamp plots will be plotted in a single plot. If
3100 False, each postage stamp plot will be plotted separately. Is only valid
3101 if stamp_coordinate exists and has more than a single point.
3103 Returns
3104 -------
3105 iris.cube.Cube | iris.cube.CubeList
3106 The original Cube or CubeList (so further operations can be applied).
3107 Plotted data.
3109 Raises
3110 ------
3111 ValueError
3112 If the cube doesn't have the right dimensions.
3113 TypeError
3114 If the cube isn't a Cube or CubeList.
3115 """
3116 recipe_title = get_recipe_metadata().get("title", "Histogram")
3118 cubes = iter_maybe(cubes)
3120 # Internal plotting function.
3121 plotting_func = _plot_and_save_histogram_series
3123 num_models = get_num_models(cubes)
3125 validate_cube_shape(cubes, num_models)
3127 # If several histograms are plotted, check sequence_coordinate
3128 check_sequence_coordinate(cubes, sequence_coordinate)
3130 # Get axis minimum and maximum from levels information.
3131 # If no levels set, derive minima and maxima from data in CubeList.
3132 vmin, vmax = _set_axis_range(cubes)
3134 # Make postage stamp plots if stamp_coordinate exists and has more than a
3135 # single point. If single_plot is True:
3136 # -- all postage stamp plots will be plotted in a single plot instead of
3137 # separate postage stamp plots.
3138 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3139 # produced per single model only
3140 if num_models == 1:
3141 if ( 3141 ↛ 3145line 3141 didn't jump to line 3145 because the condition on line 3141 was never true
3142 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3143 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3144 ):
3145 if single_plot:
3146 plotting_func = (
3147 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3148 )
3149 else:
3150 plotting_func = _plot_and_save_postage_stamp_histogram_series
3151 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3152 else:
3153 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3155 plot_index = []
3156 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3157 # Create a plot for each value of the sequence coordinate. Allowing for
3158 # multiple cubes in a CubeList to be plotted in the same plot for similar
3159 # sequence values. Passing a CubeList into the internal plotting function
3160 # for similar values of the sequence coordinate. cube_slice can be an
3161 # iris.cube.Cube or an iris.cube.CubeList.
3162 for cube_slice in cube_iterables:
3163 single_cube = cube_slice
3164 if isinstance(cube_slice, iris.cube.CubeList):
3165 single_cube = cube_slice[0]
3167 # Ensure valid stamp coordinate in cube dimensions
3168 if stamp_coordinate == "realization": 3168 ↛ 3171line 3168 didn't jump to line 3171 because the condition on line 3168 was always true
3169 stamp_coordinate = check_stamp_coordinate(single_cube)
3170 # Set plot titles and filename, based on sequence coordinate
3171 seq_coord = single_cube.coord(sequence_coordinate)
3172 # Use time coordinate in title and filename if single histogram output.
3173 if sequence_coordinate == "realization" and nplot == 1: 3173 ↛ 3174line 3173 didn't jump to line 3174 because the condition on line 3173 was never true
3174 seq_coord = single_cube.coord("time")
3175 # Use station name in title and filename if model vs obs comparison
3176 if sequence_coordinate == "station": 3176 ↛ 3177line 3176 didn't jump to line 3177 because the condition on line 3176 was never true
3177 seq_coord = single_cube.coord("Station_Name")
3179 plot_title, plot_filename = _set_title_and_filename(
3180 seq_coord, nplot, recipe_title, filename
3181 )
3183 # Do the actual plotting.
3184 plotting_func(
3185 cube_slice,
3186 filename=plot_filename,
3187 stamp_coordinate=stamp_coordinate,
3188 title=plot_title,
3189 vmin=vmin,
3190 vmax=vmax,
3191 )
3192 plot_index.append(plot_filename)
3194 # Add list of plots to plot metadata.
3195 complete_plot_index = _append_to_plot_index(plot_index)
3197 # Make a page to display the plots.
3198 _make_plot_html_page(complete_plot_index)
3200 return cubes
3203def plot_scatter_series(
3204 cubes: iris.cube.Cube | iris.cube.CubeList,
3205 filename: str | None = None,
3206 sequence_coordinate: str = "time",
3207 stamp_coordinate: str = "realization",
3208 hexbin: bool = False,
3209 **kwargs,
3210) -> iris.cube.Cube | iris.cube.CubeList:
3211 """Plot a scatter plot for each sequence coordinate provided.
3213 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3214 is present then a sequence of plots will be produced using the time slider
3215 functionality to scroll through scatter against time. If a
3216 stamp_coordinate is present then postage stamp plots will be produced. If
3217 stamp_coordinate and single_plot is True, all postage stamp plots will be
3218 plotted in a single plot instead of separate postage stamp plots.
3220 Parameters
3221 ----------
3222 cubes: Cube | iris.cube.CubeList
3223 Iris cube or CubeList of the data to plot. It should have a single dimension other
3224 than the stamp coordinate.
3225 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3226 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3227 filename: str, optional
3228 Name of the plot to write, used as a prefix for plot sequences. Defaults
3229 to the recipe name.
3230 sequence_coordinate: str, optional
3231 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3232 This coordinate must exist in the cube and will be used for the time
3233 slider.
3234 stamp_coordinate: str, optional
3235 Coordinate about which to plot postage stamp plots. Defaults to
3236 ``"realization"``.
3237 hexbin: bool, optional
3238 If True, generate hexbin comparison plot.
3239 If False, generate point-by-point scatter plot.
3241 Returns
3242 -------
3243 iris.cube.Cube | iris.cube.CubeList
3244 The original Cube or CubeList (so further operations can be applied).
3245 Plotted data.
3247 Raises
3248 ------
3249 ValueError
3250 If the cube doesn't have the right dimensions.
3251 TypeError
3252 If the cube isn't a Cube or CubeList.
3253 """
3254 recipe_title = get_recipe_metadata().get("title", "Scatter")
3256 cubes = iter_maybe(cubes)
3258 # Internal plotting function.
3259 plotting_func = _plot_and_save_scatter_series
3261 num_models = get_num_models(cubes)
3263 validate_cube_shape(cubes, num_models)
3265 check_sequence_coordinate(cubes, sequence_coordinate)
3267 vmin, vmax = _set_axis_range(cubes)
3269 # Require >1 models to compare on scatter plot
3270 if num_models > 1:
3271 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3272 else:
3273 raise ValueError(
3274 "Scatter plot series requires multiple number of models in input data."
3275 )
3277 plot_index = []
3278 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3279 # Create a plot for each value of the sequence coordinate. Allowing for
3280 # multiple cubes in a CubeList to be plotted in the same plot for similar
3281 # sequence values. Passing a CubeList into the internal plotting function
3282 # for similar values of the sequence coordinate. cube_slice can be an
3283 # iris.cube.Cube or an iris.cube.CubeList.
3284 for cube_slice in cube_iterables:
3285 single_cube = cube_slice
3286 if isinstance(cube_slice, iris.cube.CubeList): 3286 ↛ 3290line 3286 didn't jump to line 3290 because the condition on line 3286 was always true
3287 single_cube = cube_slice[0]
3289 # Ensure valid stamp coordinate in cube dimensions
3290 if stamp_coordinate == "realization": 3290 ↛ 3293line 3290 didn't jump to line 3293 because the condition on line 3290 was always true
3291 stamp_coordinate = check_stamp_coordinate(single_cube)
3292 # Set plot titles and filename, based on sequence coordinate
3293 seq_coord = single_cube.coord(sequence_coordinate)
3294 # Use time coordinate in title and filename if single histogram output.
3295 if sequence_coordinate == "realization" and nplot == 1:
3296 seq_coord = single_cube.coord("time")
3297 # Use station name in title and filename if model vs obs comparison
3298 if sequence_coordinate == "station":
3299 seq_coord = single_cube.coord("Station_Name")
3301 plot_title, plot_filename = _set_title_and_filename(
3302 seq_coord, nplot, recipe_title, filename
3303 )
3305 # Do the actual plotting.
3306 plotting_func(
3307 cube_slice,
3308 filename=plot_filename,
3309 stamp_coordinate=stamp_coordinate,
3310 title=plot_title,
3311 vmin=vmin,
3312 vmax=vmax,
3313 hexbin=hexbin,
3314 )
3315 plot_index.append(plot_filename)
3317 # Add list of plots to plot metadata.
3318 complete_plot_index = _append_to_plot_index(plot_index)
3320 # Make a page to display the plots.
3321 _make_plot_html_page(complete_plot_index)
3323 return cubes
3326def _plot_and_save_postage_stamp_power_spectrum_series(
3327 cubes: iris.cube.Cube,
3328 coords: list[iris.coords.Coord],
3329 stamp_coordinate: str,
3330 filename: str,
3331 title: str,
3332 series_coordinate: str | None = None,
3333 **kwargs,
3334):
3335 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3337 Parameters
3338 ----------
3339 cubes: Cube or CubeList
3340 Cube or Cubelist of the power spectrum data.
3341 coords: list[Coord]
3342 Coordinates to plot on the x-axis, one per cube.
3343 stamp_coordinate: str
3344 Coordinate that becomes different plots.
3345 filename: str
3346 Filename of the plot to write.
3347 title: str
3348 Plot title.
3349 series_coordinate: str, optional
3350 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3352 """
3353 # Use the smallest square grid that will fit the members.
3354 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3356 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3357 model_colors_map = get_model_colors_map(cubes)
3358 # ax = plt.gca()
3359 # Make a subplot for each member.
3360 for member, subplot in zip(
3361 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3362 ):
3363 ax = plt.subplot(grid_size, grid_size, subplot)
3365 # Store min/max ranges.
3366 y_levels = []
3368 line_marker = None
3369 line_width = 1
3371 for cube in iter_maybe(member):
3372 xcoord = _select_series_coord(cube, series_coordinate)
3373 xname = xcoord.points
3375 yfield = cube.data # power spectrum
3376 label = None
3377 color = "black"
3378 if model_colors_map: 3378 ↛ 3379line 3378 didn't jump to line 3379 because the condition on line 3378 was never true
3379 label = cube.attributes.get("model_name")
3380 color = model_colors_map.get(label)
3382 if member.coord(stamp_coordinate).points == [0]:
3383 ax.plot(
3384 xname,
3385 yfield,
3386 color=color,
3387 marker=line_marker,
3388 ls="-",
3389 lw=line_width,
3390 label=f"{label} (control)"
3391 if len(cube.coord(stamp_coordinate).points) > 1
3392 else label,
3393 )
3394 # Label with member if part of an ensemble and not the control.
3395 else:
3396 ax.plot(
3397 xname,
3398 yfield,
3399 color=color,
3400 ls="-",
3401 lw=1.5,
3402 alpha=0.75,
3403 label=f"{label} (member)",
3404 )
3406 # Calculate the global min/max if multiple cubes are given.
3407 _, levels, _ = colorbar_map_levels(cube, axis="y")
3408 if levels is not None: 3408 ↛ 3409line 3408 didn't jump to line 3409 because the condition on line 3408 was never true
3409 y_levels.append(min(levels))
3410 y_levels.append(max(levels))
3412 # Add some labels and tweak the style.
3413 title = f"{title}"
3414 ax.set_title(title, fontsize=16)
3416 # Set appropriate x-axis label based on coordinate
3417 if series_coordinate == "wavelength" or ( 3417 ↛ 3420line 3417 didn't jump to line 3420 because the condition on line 3417 was never true
3418 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3419 ):
3420 ax.set_xlabel("Wavelength (km)", fontsize=14)
3421 elif series_coordinate == "physical_wavenumber" or ( 3421 ↛ 3426line 3421 didn't jump to line 3426 because the condition on line 3421 was always true
3422 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3423 ):
3424 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3425 else: # frequency or check units
3426 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3427 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3428 else:
3429 ax.set_xlabel("Wavenumber", fontsize=14)
3431 ax.set_ylabel("Power Spectral Density", fontsize=14)
3432 ax.tick_params(axis="both", labelsize=12)
3434 # Set log-log scale
3435 ax.set_xscale("log")
3436 ax.set_yscale("log")
3438 # Add gridlines
3439 ax.grid(linestyle="--", color="grey", linewidth=1)
3440 # Ientify unique labels for legend
3441 handles = list(
3442 {
3443 label: handle
3444 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3445 }.values()
3446 )
3447 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3449 ax = plt.gca()
3450 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3452 # Save plot.
3453 _save_close_figure(fig, "histogram postage stamp", filename)
3456def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3457 cubes: iris.cube.Cube,
3458 coords: list[iris.coords.Coord],
3459 stamp_coordinate: str,
3460 filename: str,
3461 title: str,
3462 series_coordinate: str | None = None,
3463 **kwargs,
3464):
3465 """Plot and save power spectra for ensemble members in single plot.
3467 Parameters
3468 ----------
3469 cubes: Cube or CubeList
3470 Cube or Cubelist of the power spectrum data.
3471 coords: list[Coord]
3472 Coordinates to plot on the x-axis, one per cube.
3473 stamp_coordinate: str
3474 Coordinate that becomes different plots.
3475 filename: str
3476 Filename of the plot to write.
3477 title: str
3478 Plot title.
3479 series_coordinate: str, optional
3480 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3482 """
3483 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3484 model_colors_map = get_model_colors_map(cubes)
3486 line_marker = None
3487 line_width = 1
3489 # Compute ensemble statistics to show spread
3490 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3491 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3492 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3494 xcoord_global = mean_cube.coord(series_coordinate)
3495 x_global = xcoord_global.points
3497 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3498 xcoord = _select_series_coord(member, series_coordinate)
3499 xname = xcoord.points
3501 yfield = member.data # power spectrum
3502 color = "black"
3503 if model_colors_map: 3503 ↛ 3507line 3503 didn't jump to line 3507 because the condition on line 3503 was always true
3504 label = member.attributes.get("model_name") if i == 0 else None
3505 color = model_colors_map.get(label)
3507 if member.coord(stamp_coordinate).points == [0]:
3508 ax.plot(
3509 xname,
3510 yfield,
3511 color=color,
3512 marker=line_marker,
3513 ls="-",
3514 lw=line_width,
3515 label=f"{label} (control)"
3516 if len(member.coord(stamp_coordinate).points) > 1
3517 else label,
3518 )
3519 # Label with member number if part of an ensemble and not the control.
3520 else:
3521 ax.plot(
3522 xname,
3523 yfield,
3524 color=color,
3525 ls="-",
3526 lw=1.5,
3527 alpha=0.75,
3528 label=label,
3529 )
3531 # Set appropriate x-axis label based on coordinate
3532 if series_coordinate == "wavelength" or ( 3532 ↛ 3535line 3532 didn't jump to line 3535 because the condition on line 3532 was never true
3533 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3534 ):
3535 ax.set_xlabel("Wavelength (km)", fontsize=14)
3536 elif series_coordinate == "physical_wavenumber" or ( 3536 ↛ 3541line 3536 didn't jump to line 3541 because the condition on line 3536 was always true
3537 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3538 ):
3539 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3540 else: # frequency or check units
3541 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3542 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3543 else:
3544 ax.set_xlabel("Wavenumber", fontsize=14)
3546 # Add ensemble spread shading
3547 ax.fill_between(
3548 x_global,
3549 min_cube.data,
3550 max_cube.data,
3551 color="grey",
3552 alpha=0.3,
3553 label="Ensemble spread",
3554 )
3556 # Add ensemble mean line
3557 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3559 ax.set_ylabel("Power Spectral Density", fontsize=14)
3560 ax.tick_params(axis="both", labelsize=12)
3562 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3563 # Set log-log scale
3564 ax.set_xscale("log")
3565 ax.set_yscale("log")
3567 # Add gridlines
3568 ax.grid(linestyle="--", color="grey", linewidth=1)
3569 # Identify unique labels for legend
3570 handles = list(
3571 {
3572 label: handle
3573 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3574 }.values()
3575 )
3576 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3578 # Figure title.
3579 ax.set_title(title, fontsize=16)
3581 # Save plot.
3582 _save_close_figure(fig, "power spectra postage stamp", filename)