Coverage for src/CSET/operators/plot.py: 84%
1171 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-24 15:03 +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 calc_array_stats,
55 check_sequence_coordinate,
56 check_single_cube,
57 check_stamp_coordinate,
58 fully_equalise_attributes,
59 get_cube_yxcoordname,
60 get_num_models,
61 is_transect,
62 slice_over_maybe,
63 validate_cube_shape,
64 validate_cubes_coords,
65)
66from CSET.operators.collapse import collapse
67from CSET.operators.misc import _extract_common_time_points
68from CSET.operators.regrid import regrid_onto_cube
70logger = logging.getLogger(__name__)
72# Use a non-interactive plotting backend.
73mpl.use("agg")
76############################
77# Private helper functions #
78############################
81def in_sphinx_gallery():
82 """Test if running plot code in sphinx-gallery context."""
83 return "sphinx_gallery" in sys.modules
86def _append_to_plot_index(plot_index: list) -> list:
87 """Add plots into the plot index, returning the complete plot index."""
88 with open("meta.json", "r+t", encoding="UTF-8") as fp:
89 fcntl.flock(fp, fcntl.LOCK_EX)
90 fp.seek(0)
91 meta = json.load(fp)
92 complete_plot_index = meta.get("plots", [])
93 complete_plot_index = complete_plot_index + plot_index
94 meta["plots"] = complete_plot_index
95 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
96 os.getenv("DO_CASE_AGGREGATION")
97 ):
98 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
99 fp.seek(0)
100 fp.truncate()
101 json.dump(meta, fp, indent=2)
102 return complete_plot_index
105def _make_plot_html_page(plots: list):
106 """Create a HTML page to display a plot image."""
107 # Debug check that plots actually contains some strings.
108 assert isinstance(plots[0], str)
110 # Load HTML template file.
111 operator_files = importlib.resources.files()
112 template_file = operator_files.joinpath("_plot_page_template.html")
114 # Get some metadata.
115 meta = get_recipe_metadata()
116 title = meta.get("title", "Untitled")
117 description = MarkdownIt().render(meta.get("description", "*No description.*"))
119 # Prepare template variables.
120 variables = {
121 "title": title,
122 "description": description,
123 "initial_plot": plots[0],
124 "plots": plots,
125 "title_slug": slugify(title),
126 }
128 # Render template.
129 html = render_file(template_file, **variables)
131 # Save completed HTML.
132 with open("index.html", "wt", encoding="UTF-8") as fp:
133 fp.write(html)
136def _save_close_figure(figure, plot_type: str, filename: str):
137 """Save generated plot figure file and close figure.
139 If running documentation gallery generation, avoid saving to file.
141 Parameters
142 ----------
143 figure:
144 Matplotlib Figure object holding all plot elements.
145 plot_type: str
146 String identifier for plot type for logging information.
147 filename: str
148 Filename for saved figure.
149 """
150 if not in_sphinx_gallery():
151 figure.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
152 logger.info("Saved %s plot to %s", plot_type, filename)
153 plt.close(figure)
156def _setup_spatial_map(
157 cube: iris.cube.Cube,
158 figure,
159 cmap,
160 grid_size: tuple[int, int] | None = None,
161 subplot: int | None = None,
162):
163 """Define map projections, extent and add coastlines and borderlines for spatial plots.
165 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
166 is specified, and map extent defined based on the input data.
168 Parameters
169 ----------
170 cube: Cube
171 2 dimensional (lat and lon) Cube of the data to plot.
172 figure:
173 Matplotlib Figure object holding all plot elements.
174 cmap:
175 Matplotlib colormap.
176 grid_size: (int, int), optional
177 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
178 subplot: int, optional
179 Subplot index if multiple spatial subplots in figure.
181 Returns
182 -------
183 axes:
184 Matplotlib GeoAxes definition.
185 """
186 # Identify min/max plot bounds.
187 try:
188 lat_axis, lon_axis = get_cube_yxcoordname(cube)
189 xmin = np.nanmin(cube.coord(lon_axis).points)
190 xmax = np.nanmax(cube.coord(lon_axis).points)
191 ymin = np.nanmin(cube.coord(lat_axis).points)
192 ymax = np.nanmax(cube.coord(lat_axis).points)
194 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
195 if np.abs(xmax - xmin) > 180.0:
196 xmin = xmin - 180.0
197 xmax = xmax - 180.0
198 logger.debug("Adjusting plot bounds to fit global extent.")
200 # Consider map projection orientation.
201 # Adapting orientation enables plotting across international dateline.
202 # Users can adapt the default central_longitude if alternative projections views.
203 if xmax > 180.0 or xmin < -180.0:
204 central_longitude = 180.0
205 else:
206 central_longitude = 0.0
208 # Define spatial map projection.
209 coord_system = cube.coord(lat_axis).coord_system
210 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
211 # Define rotated pole map projection for rotated pole inputs.
212 projection = ccrs.RotatedPole(
213 pole_longitude=coord_system.grid_north_pole_longitude,
214 pole_latitude=coord_system.grid_north_pole_latitude,
215 central_rotated_longitude=central_longitude,
216 )
217 crs = projection
218 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 218 ↛ 220line 218 didn't jump to line 220 because the condition on line 218 was never true
219 # Define Transverse Mercator projection for TM inputs.
220 projection = ccrs.TransverseMercator(
221 central_longitude=coord_system.longitude_of_central_meridian,
222 central_latitude=coord_system.latitude_of_projection_origin,
223 false_easting=coord_system.false_easting,
224 false_northing=coord_system.false_northing,
225 scale_factor=coord_system.scale_factor_at_central_meridian,
226 )
227 crs = projection
228 else:
229 # Assume polar projection for regional grids encompassing N. Pole
230 if ymin > 20.0 and ymax > 80.0:
231 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
232 elif ymin < -80.0 and ymax < -20.0:
233 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
234 # Define regular map projection for non-rotated pole inputs.
235 # Alternatives might include e.g. for global model outputs:
236 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
237 # projection = ccrs.NearsidePerspective(
238 # central_longitude=180.0,
239 # central_latitude=0,
240 # satellite_height=35785831,
241 # )
242 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
243 else:
244 projection = ccrs.PlateCarree(central_longitude=central_longitude)
245 crs = ccrs.PlateCarree()
247 # Define axes for plot (or subplot) with required map projection.
248 if subplot is not None:
249 axes = figure.add_subplot(
250 grid_size[0], grid_size[1], subplot, projection=projection
251 )
252 else:
253 axes = figure.add_subplot(projection=projection)
255 # Add coastlines and borderlines if cube contains x and y map coordinates.
256 # Avoid adding lines for specific fixed ancillary spatial plots
257 if any(name in cube.name() for name in ("land_", "orography", "altitude")):
258 pass
259 else:
260 if cmap.name in ["viridis", "Greys"]:
261 coastcol = "magenta"
262 else:
263 coastcol = "black"
264 logger.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
265 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
266 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
268 # Add gridlines.
269 gl = axes.gridlines(
270 alpha=0.3,
271 draw_labels=True,
272 dms=False,
273 x_inline=False,
274 y_inline=False,
275 )
276 gl.top_labels = False
277 gl.right_labels = False
278 if subplot:
279 gl.bottom_labels = False
280 gl.left_labels = False
281 if subplot % grid_size[1] == 1:
282 gl.left_labels = True
283 if subplot > ((grid_size[0] - 1) * grid_size[1]): 283 ↛ 288line 283 didn't jump to line 288 because the condition on line 283 was always true
284 gl.bottom_labels = True
286 # If is lat/lon spatial map, fix extent to keep plot tight.
287 # Specifying crs within set_extent helps ensure only data region is shown.
288 if isinstance(
289 coord_system, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
290 ):
291 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
293 except ValueError:
294 # Skip if not both x and y map coordinates.
295 axes = figure.gca()
297 return axes
300def _get_plot_resolution() -> int:
301 """Get resolution of rasterised plots in pixels per inch."""
302 return get_recipe_metadata().get("plot_resolution", 100)
305def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
306 """Return title and filename based on start and end points or bounds."""
307 if use_bounds and seq_coord.has_bounds():
308 vals = seq_coord.bounds.flatten()
309 else:
310 vals = seq_coord.points
311 start = seq_coord.units.title(vals[0])
312 end = seq_coord.units.title(vals[-1])
314 if start == end:
315 sequence_title = f"\n [{start}]"
316 sequence_fname = f"_{filename_slugify(start)}"
317 else:
318 sequence_title = f"\n [{start} to {end}]"
319 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
321 # Do not include time if coord set to zero.
322 if (
323 seq_coord.units == "hours since 0001-01-01 00:00:00"
324 and vals[0] == 0
325 and vals[-1] == 0
326 ):
327 sequence_title = ""
328 sequence_fname = ""
330 return sequence_title, sequence_fname
333def _set_title_and_filename(
334 seq_coord: iris.coords.Coord,
335 nplot: int,
336 recipe_title: str,
337 filename: str,
338 model_name: str | None = None,
339):
340 """Set plot title and filename based on cube coordinate.
342 Parameters
343 ----------
344 sequence_coordinate: iris.coords.Coord
345 Coordinate about which to make a plot sequence.
346 nplot: int
347 Number of output plots to generate - controls title/naming.
348 recipe_title: str
349 Default plot title, potentially to update.
350 filename: str
351 Input plot filename, potentially to update.
353 Returns
354 -------
355 plot_title: str
356 Output formatted plot title string, based on plotted data.
357 plot_filename: str
358 Output formatted plot filename string.
359 """
360 ndim = seq_coord.ndim
361 npoints = np.size(seq_coord.points)
362 sequence_title = ""
363 sequence_fname = ""
365 # Case 1: Multiple dimension sequence input - list number of aggregated cases
366 # (e.g. aggregation histogram plots)
367 if ndim > 1:
368 ncase = np.shape(seq_coord)[0]
369 sequence_title = f"\n [{ncase} cases]"
370 sequence_fname = f"_{ncase}cases"
372 # Case 2: Single dimension input
373 else:
374 # Single sequence point
375 if npoints == 1:
376 if nplot > 1:
377 # Default labels for sequence inputs
378 sequence_value = seq_coord.units.title(seq_coord.points[0])
379 sequence_value = sequence_value.replace(" unknown", "")
380 sequence_title = f"\n [{sequence_value}]"
381 sequence_fname = f"_{filename_slugify(sequence_value)}"
382 else:
383 # Aggregated attribute available where input collapsed over aggregation
384 try:
385 ncase = seq_coord.attributes["number_reference_times"]
386 sequence_title = f"\n [{ncase} cases]"
387 sequence_fname = f"_{ncase}cases"
388 except KeyError:
389 sequence_title, sequence_fname = _get_start_end_strings(
390 seq_coord, use_bounds=seq_coord.has_bounds()
391 )
392 # Multiple sequence (e.g. time) points
393 else:
394 sequence_title, sequence_fname = _get_start_end_strings(
395 seq_coord, use_bounds=False
396 )
398 # Set plot title and filename
399 plot_title = f"{recipe_title}{sequence_title}"
401 # Set plot filename, defaulting to user input if provided.
402 if filename is None:
403 filename = slugify(recipe_title)
404 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
405 else:
406 if nplot > 1:
407 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
408 else:
409 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
411 if model_name: 411 ↛ 412line 411 didn't jump to line 412 because the condition on line 411 was never true
412 plot_filename = f"{model_name}_{plot_filename}"
413 plot_title = f"{model_name}_{plot_title}"
415 return plot_title, plot_filename
418def _select_series_coord(cube, series_coordinate):
419 """Determine the grid coordinates to use to calculate grid spacing."""
420 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
421 if series_coordinate in spacing_coordinates: 421 ↛ 427line 421 didn't jump to line 427 because the condition on line 421 was always true
422 # Try the requested coordinate first then the fallbacks in order.
423 fallbacks = [series_coordinate] + [
424 c for c in spacing_coordinates if c != series_coordinate
425 ]
426 else:
427 fallbacks = {series_coordinate}
429 # Try each possible coordinate.
430 for coord in fallbacks:
431 try:
432 return cube.coord(coord)
433 except iris.exceptions.CoordinateNotFoundError:
434 logger.debug("Coordinate %s not found.", coord)
436 # If we get here, none of the fallback options were found.
437 raise iris.exceptions.CoordinateNotFoundError(
438 f"No valid coordinate found for '{series_coordinate}' "
439 f"or fallback options {fallbacks}"
440 )
443def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
444 """Control postage stamp plot output titles based on stamp coordinate."""
445 if stamp_coord.name() == "realization":
446 mtitle = "Member"
447 else:
448 mtitle = stamp_coord.name().capitalize()
450 if stamp_coord.name() == "time":
451 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
452 else:
453 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
455 return mtitle
458def _set_axis_range(cubes):
459 """Get minimum and maximum from levels information."""
460 levels = None
461 for cube in cubes: 461 ↛ 477line 461 didn't jump to line 477 because the loop on line 461 didn't complete
462 # First check if user-specified "auto" range variable.
463 # This maintains the value of levels as None, so proceed.
464 _, levels, _ = colorbar_map_levels(cube, axis="y")
465 if levels is None:
466 break
467 # If levels is changed, recheck to use the vmin,vmax or
468 # levels-based ranges for histogram plots.
469 _, levels, _ = colorbar_map_levels(cube)
470 logger.debug("levels: %s", levels)
471 if levels is not None: 471 ↛ 461line 471 didn't jump to line 461 because the condition on line 471 was always true
472 vmin = min(levels)
473 vmax = max(levels)
474 logger.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
475 break
477 if levels is None:
478 vmin = min(cb.data.min() for cb in cubes)
479 vmax = max(cb.data.max() for cb in cubes)
481 return vmin, vmax
484def _find_matched_slices(cubes, sequence_coordinate):
485 """Identify matched cubes in CubeList by sequence_coordinate values.
487 Ensures common points are compared for multiple cube inputs.
488 """
489 all_points = sorted(
490 set(
491 itertools.chain.from_iterable(
492 cb.coord(sequence_coordinate).points for cb in cubes
493 )
494 )
495 )
496 all_slices = list(
497 itertools.chain.from_iterable(
498 cb.slices_over(sequence_coordinate) for cb in cubes
499 )
500 )
501 # Matched slices (matched by seq coord point; it may happen that
502 # evaluated models do not cover the same seq coord range, hence matching
503 # necessary)
504 cube_iterables = [
505 iris.cube.CubeList(
506 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
507 )
508 for point in all_points
509 ]
511 return cube_iterables
514def _plot_and_save_spatial_plot(
515 cube: iris.cube.Cube,
516 filename: str,
517 title: str,
518 method: Literal["contourf", "pcolormesh", "scatter"],
519 overlay_cube: iris.cube.Cube | None = None,
520 contour_cube: iris.cube.Cube | None = None,
521 point_cube: iris.cube.Cube | None = None,
522 **kwargs,
523):
524 """Plot and save a spatial plot.
526 Parameters
527 ----------
528 cube: Cube
529 2 dimensional (lat and lon) Cube of the data to plot.
530 filename: str
531 Filename of the plot to write.
532 title: str
533 Plot title.
534 method: "contourf" | "pcolormesh" | "scatter"
535 The plotting method to use
536 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
537 overlay_cube: Cube, optional
538 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
539 contour_cube: Cube, optional
540 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
541 point_cube: Cube, optional
542 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
543 """
544 # Setup plot details, size, resolution, etc.
545 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
547 # Specify the color bar
548 cmap, levels, norm = colorbar_map_levels(cube)
550 # If overplotting, set required colorbars
551 if overlay_cube:
552 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
553 if contour_cube:
554 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
556 # Setup plot map projection, extent and coastlines and borderlines.
557 axes = _setup_spatial_map(cube, fig, cmap)
559 # Set colorscale bounds
560 try:
561 vmin = min(levels)
562 vmax = max(levels)
563 except TypeError:
564 vmin, vmax = None, None
565 # Ensure to use norm and not vmin/vmax if levels are defined.
566 if norm is not None:
567 vmin = None
568 vmax = None
569 logger.debug("Plotting using defined levels.")
571 # Plot the field.
572 if method == "contourf":
573 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
574 elif method == "pcolormesh":
575 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
576 elif method == "scatter":
577 # Scatter plot of the field. The marker size is chosen to give
578 # symbols that decrease in size as the number of data points
579 # increases, although the fraction of the figure covered by
580 # symbols increases roughly as N^(1/2), disregarding overlaps,
581 # and has been selected for the default figure size of (10, 10).
582 # Should this be changed, the marker size should be adjusted in
583 # proportion to the area of the figure.
584 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
585 lat_axis, lon_axis = get_cube_yxcoordname(cube)
586 plot = iplt.scatter(
587 cube.coord(lon_axis),
588 cube.coord(lat_axis),
589 c=cube.data[:],
590 s=mrk_size,
591 cmap=cmap,
592 edgecolors="k",
593 norm=norm,
594 vmin=vmin,
595 vmax=vmax,
596 )
597 else:
598 raise ValueError(f"Unknown plotting method: {method}")
600 # Overplot overlay field, if required
601 if overlay_cube:
602 try:
603 over_vmin = min(over_levels)
604 over_vmax = max(over_levels)
605 except TypeError:
606 over_vmin, over_vmax = None, None
607 if over_norm is not None: 607 ↛ 608line 607 didn't jump to line 608 because the condition on line 607 was never true
608 over_vmin = None
609 over_vmax = None
610 overlay = iplt.pcolormesh(
611 overlay_cube,
612 cmap=over_cmap,
613 norm=over_norm,
614 alpha=0.8,
615 vmin=over_vmin,
616 vmax=over_vmax,
617 )
618 # Overplot contour field, if required, with contour labelling.
619 if contour_cube:
620 contour = iplt.contour(
621 contour_cube,
622 colors="darkgray",
623 levels=cntr_levels,
624 norm=cntr_norm,
625 alpha=0.5,
626 linestyles="--",
627 linewidths=1,
628 )
629 plt.clabel(contour)
630 # Overplot valid elements of point-based field, if required.
631 # Check for non-masked points only to avoid plotting missing data.
632 if point_cube:
633 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
634 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
635 lon_coord = point_cube.coord(lon_axis)
636 lat_coord = point_cube.coord(lat_axis)
637 valid = ~point_cube.data.mask
638 valid_lon = iris.coords.AuxCoord(
639 lon_coord.points[valid],
640 standard_name=lon_coord.standard_name,
641 units=lon_coord.units,
642 coord_system=lon_coord.coord_system,
643 )
644 valid_lat = iris.coords.AuxCoord(
645 lat_coord.points[valid],
646 standard_name=lat_coord.standard_name,
647 units=lat_coord.units,
648 coord_system=lat_coord.coord_system,
649 )
650 iplt.scatter(
651 valid_lon,
652 valid_lat,
653 c=point_cube.data[valid],
654 s=mrk_size,
655 cmap=cmap,
656 edgecolors="k",
657 norm=norm,
658 vmin=vmin,
659 vmax=vmax,
660 )
662 # Check to see if transect, and if so, adjust y axis.
663 if is_transect(cube):
664 if "pressure" in [coord.name() for coord in cube.coords()]:
665 axes.invert_yaxis()
666 axes.set_yscale("log")
667 axes.set_ylim(1100, 100)
668 # If both model_level_number and level_height exists, iplt can construct
669 # plot as a function of height above orography (NOT sea level).
670 elif {"model_level_number", "level_height"}.issubset( 670 ↛ 675line 670 didn't jump to line 675 because the condition on line 670 was always true
671 {coord.name() for coord in cube.coords()}
672 ):
673 axes.set_yscale("log")
675 axes.set_title(
676 f"{title}\n"
677 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
678 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
679 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
680 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
681 fontsize=16,
682 )
684 # Inset code
685 axins = inset_axes(
686 axes,
687 width="20%",
688 height="20%",
689 loc="upper right",
690 axes_class=GeoAxes,
691 axes_kwargs={"map_projection": ccrs.PlateCarree()},
692 )
694 # Slightly transparent to reduce plot blocking.
695 axins.patch.set_alpha(0.4)
697 axins.coastlines(resolution="50m")
698 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
700 SLat, SLon, ELat, ELon = (
701 float(coord) for coord in cube.attributes["transect_coords"].split("_")
702 )
704 # Draw line between them
705 axins.plot(
706 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
707 )
709 # Plot points (note: lon, lat order for Cartopy)
710 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
711 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
713 lon_min, lon_max = sorted([SLon, ELon])
714 lat_min, lat_max = sorted([SLat, ELat])
716 # Midpoints
717 lon_mid = (lon_min + lon_max) / 2
718 lat_mid = (lat_min + lat_max) / 2
720 # Maximum half-range
721 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
722 if half_range == 0: # points identical → provide small default 722 ↛ 726line 722 didn't jump to line 726 because the condition on line 722 was always true
723 half_range = 1
725 # Set square extent
726 axins.set_extent(
727 [
728 lon_mid - half_range,
729 lon_mid + half_range,
730 lat_mid - half_range,
731 lat_mid + half_range,
732 ],
733 crs=ccrs.PlateCarree(),
734 )
736 # Ensure square aspect
737 axins.set_aspect("equal")
739 else:
740 # Add title.
741 axes.set_title(title, fontsize=16)
743 # Adjust padding if spatial plot or transect
744 if is_transect(cube):
745 yinfopad = -0.1
746 ycbarpad = 0.1
747 else:
748 yinfopad = 0.01
749 ycbarpad = 0.042
751 # Add watermark with min/max/mean. Currently not user togglable.
752 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
753 cube_min, cube_max, cube_mean = calc_array_stats(cube.data)
754 axes.annotate(
755 f"Min: {np.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(cube.data):.3g}",
756 xy=(0.025, yinfopad),
757 xycoords="axes fraction",
758 xytext=(-5, 5),
759 textcoords="offset points",
760 ha="left",
761 va="bottom",
762 size=11,
763 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
764 )
766 # Add secondary colour bar for overlay_cube field if required.
767 if overlay_cube:
768 cbarB = fig.colorbar(
769 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
770 )
771 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
772 # add ticks and tick_labels for every levels if less than 20 levels exist
773 if over_levels is not None and len(over_levels) < 20: 773 ↛ 774line 773 didn't jump to line 774 because the condition on line 773 was never true
774 cbarB.set_ticks(over_levels)
775 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
776 if any(
777 var in overlay_cube.name()
778 for var in ("rainfall", "snowfall", "visibility")
779 ):
780 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
781 logger.debug("Set secondary colorbar ticks and labels.")
783 # Add main colour bar.
784 cbar = fig.colorbar(
785 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
786 )
788 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
789 # add ticks and tick_labels for every levels if less than 20 levels exist
790 if levels is not None and len(levels) < 20:
791 cbar.set_ticks(levels)
792 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
793 if any(var in cube.name() for var in ("rainfall", "snowfall", "visibility")): 793 ↛ 796line 793 didn't jump to line 796 because the condition on line 793 was always true
794 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
795 # Tick labels for rainfall rates from Nimrod radar data.
796 if "rainfall rate composite" in cube.name(): 796 ↛ 797line 796 didn't jump to line 797 because the condition on line 796 was never true
797 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
798 # Tick labels for rain accumulations from Nimrod radar data.
799 if "rain accumulation" in cube.name(): 799 ↛ 800line 799 didn't jump to line 800 because the condition on line 799 was never true
800 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
801 if "wts accumulation" in cube.name(): 801 ↛ 802line 801 didn't jump to line 802 because the condition on line 801 was never true
802 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
803 cbar.minorticks_off()
804 cbar.set_ticks(tick_levels)
805 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
806 cbar.set_label(label=f"{cube.name()}", size=14)
807 # Tick labels for model rainfall data.
808 if "surface_microphysical" in cube.name(): 808 ↛ 811line 808 didn't jump to line 811 because the condition on line 808 was always true
809 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
810 # Tick labels for Nimrod weights data.
811 logger.debug("Set colorbar ticks and labels.")
813 # Save plot.
814 _save_close_figure(fig, "spatial", filename)
817def _plot_and_save_postage_stamp_spatial_plot(
818 cube: iris.cube.Cube,
819 filename: str,
820 stamp_coordinate: str,
821 title: str,
822 method: Literal["contourf", "pcolormesh"],
823 overlay_cube: iris.cube.Cube | None = None,
824 contour_cube: iris.cube.Cube | None = None,
825 **kwargs,
826):
827 """Plot postage stamp spatial plots from an ensemble.
829 Parameters
830 ----------
831 cube: Cube
832 Iris cube of data to be plotted. It must have the stamp coordinate.
833 filename: str
834 Filename of the plot to write.
835 stamp_coordinate: str
836 Coordinate that becomes different plots.
837 method: "contourf" | "pcolormesh"
838 The plotting method to use.
839 overlay_cube: Cube, optional
840 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
841 contour_cube: Cube, optional
842 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
844 Raises
845 ------
846 ValueError
847 If the cube doesn't have the right dimensions.
848 """
849 # Use the smallest square grid that will fit the members.
850 nmember = len(cube.coord(stamp_coordinate).points)
851 grid_rows = int(math.sqrt(nmember))
852 grid_size = math.ceil(nmember / grid_rows)
854 fig = plt.figure(
855 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
856 )
858 # Specify the color bar
859 cmap, levels, norm = colorbar_map_levels(cube)
860 # If overplotting, set required colorbars
861 if overlay_cube: 861 ↛ 862line 861 didn't jump to line 862 because the condition on line 861 was never true
862 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
863 if contour_cube: 863 ↛ 864line 863 didn't jump to line 864 because the condition on line 863 was never true
864 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
866 # Make a subplot for each member.
867 for member, subplot in zip(
868 cube.slices_over(stamp_coordinate),
869 range(1, grid_size * grid_rows + 1),
870 strict=False,
871 ):
872 # Setup subplot map projection, extent and coastlines and borderlines.
873 axes = _setup_spatial_map(
874 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
875 )
876 if method == "contourf":
877 # Filled contour plot of the field.
878 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
879 elif method == "pcolormesh":
880 if levels is not None:
881 vmin = min(levels)
882 vmax = max(levels)
883 else:
884 raise TypeError("Unknown vmin and vmax range.")
885 vmin, vmax = None, None
886 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
887 # if levels are defined.
888 if norm is not None: 888 ↛ 889line 888 didn't jump to line 889 because the condition on line 888 was never true
889 vmin = None
890 vmax = None
891 # pcolormesh plot of the field.
892 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
893 else:
894 raise ValueError(f"Unknown plotting method: {method}")
896 # Overplot overlay field, if required
897 if overlay_cube: 897 ↛ 898line 897 didn't jump to line 898 because the condition on line 897 was never true
898 try:
899 over_vmin = min(over_levels)
900 over_vmax = max(over_levels)
901 except TypeError:
902 over_vmin, over_vmax = None, None
903 if over_norm is not None:
904 over_vmin = None
905 over_vmax = None
906 iplt.pcolormesh(
907 overlay_cube[member.coord(stamp_coordinate).points[0]],
908 cmap=over_cmap,
909 norm=over_norm,
910 alpha=0.6,
911 vmin=over_vmin,
912 vmax=over_vmax,
913 )
914 # Overplot contour field, if required
915 if contour_cube: 915 ↛ 916line 915 didn't jump to line 916 because the condition on line 915 was never true
916 iplt.contour(
917 contour_cube[member.coord(stamp_coordinate).points[0]],
918 colors="darkgray",
919 levels=cntr_levels,
920 norm=cntr_norm,
921 alpha=0.6,
922 linestyles="--",
923 linewidths=1,
924 )
925 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
926 axes.set_title(f"{mtitle}")
928 # Put the shared colorbar in its own axes.
929 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
930 colorbar = fig.colorbar(
931 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
932 )
933 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
935 # Overall figure title.
936 fig.suptitle(title, fontsize=16)
938 # Save plot.
939 _save_close_figure(fig, "contour postate stamp", filename)
942def _plot_and_save_line_series(
943 cubes: iris.cube.CubeList,
944 coords: list[iris.coords.Coord],
945 ensemble_coord: str,
946 filename: str,
947 title: str,
948 **kwargs,
949):
950 """Plot and save a 1D line series.
952 Parameters
953 ----------
954 cubes: Cube or CubeList
955 Cube or CubeList containing the cubes to plot on the y-axis.
956 coords: list[Coord]
957 Coordinates to plot on the x-axis, one per cube.
958 ensemble_coord: str
959 Ensemble coordinate in the cube.
960 filename: str
961 Filename of the plot to write.
962 title: str
963 Plot title.
964 """
965 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
967 model_colors_map = get_model_colors_map(cubes)
969 # Store min/max ranges.
970 y_levels = []
972 # Check match-up across sequence coords gives consistent sizes
973 validate_cubes_coords(cubes, coords)
975 for cube, coord in zip(cubes, coords, strict=True):
976 label = None
977 color = "black"
978 if model_colors_map:
979 label = cube.attributes.get("model_name")
980 color = model_colors_map.get(label)
981 if not cube.coords(ensemble_coord): 981 ↛ 983line 981 didn't jump to line 983 because the condition on line 981 was never true
982 # No ensemble coordinate — plot the cube directly as a single line.
983 iplt.plot(coord, cube, color=color, marker="o", ls="-", lw=3, label=label)
984 else:
985 for cube_slice in cube.slices_over(ensemble_coord):
986 # Label with (control) if part of an ensemble or not otherwise.
987 if cube_slice.coord(ensemble_coord).points == [0]:
988 iplt.plot(
989 coord,
990 cube_slice,
991 color=color,
992 marker="o",
993 ls="-",
994 lw=3,
995 label=f"{label} (control)"
996 if len(cube.coord(ensemble_coord).points) > 1
997 else label,
998 )
999 # Label with (perturbed) if part of an ensemble and not the control.
1000 else:
1001 iplt.plot(
1002 coord,
1003 cube_slice,
1004 color=color,
1005 ls="-",
1006 lw=1.5,
1007 alpha=0.75,
1008 label=f"{label} (member)",
1009 )
1011 # Calculate the global min/max if multiple cubes are given.
1012 _, levels, _ = colorbar_map_levels(cube, axis="y")
1013 if levels is not None: 1013 ↛ 1014line 1013 didn't jump to line 1014 because the condition on line 1013 was never true
1014 y_levels.append(min(levels))
1015 y_levels.append(max(levels))
1017 # Get the current axes.
1018 ax = plt.gca()
1020 # Add some labels and tweak the style.
1021 # check if cubes[0] works for single cube if not CubeList
1022 if coords[0].name() == "time":
1023 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1024 else:
1025 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1026 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1027 ax.set_title(title, fontsize=16)
1029 ax.ticklabel_format(axis="y", useOffset=False)
1030 ax.tick_params(axis="x", labelrotation=15)
1031 ax.tick_params(axis="both", labelsize=12)
1033 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1034 if y_levels: 1034 ↛ 1035line 1034 didn't jump to line 1035 because the condition on line 1034 was never true
1035 ax.set_ylim(min(y_levels), max(y_levels))
1036 logger.debug("Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels))
1037 else:
1038 ax.autoscale()
1040 # Add gridlines
1041 ax.grid(linestyle="--", color="grey", linewidth=1)
1042 # Add zero line
1043 ymin, ymax = ax.get_ylim()
1044 if ymin < 0.0 and ymax > 0.0:
1045 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1046 # Identify unique labels for legend
1047 handles = list(
1048 {
1049 label: handle
1050 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1051 }.values()
1052 )
1053 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1055 # Save plot.
1056 _save_close_figure(fig, "line", filename)
1059def _plot_and_save_line_power_spectrum_series(
1060 cubes: iris.cube.Cube | iris.cube.CubeList,
1061 coords: list[iris.coords.Coord],
1062 ensemble_coord: str,
1063 filename: str,
1064 title: str,
1065 series_coordinate: str,
1066 **kwargs,
1067):
1068 """Plot and save a 1D line series.
1070 Parameters
1071 ----------
1072 cubes: Cube or CubeList
1073 Cube or CubeList containing the cubes to plot on the y-axis.
1074 coords: list[Coord]
1075 Coordinates to plot on the x-axis, one per cube.
1076 ensemble_coord: str
1077 Ensemble coordinate in the cube.
1078 filename: str
1079 Filename of the plot to write.
1080 title: str
1081 Plot title.
1082 series_coordinate: str
1083 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1084 """
1085 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1086 model_colors_map = get_model_colors_map(cubes)
1087 ax = plt.gca()
1089 # Store min/max ranges.
1090 y_levels = []
1092 line_marker = None
1093 line_width = 1
1095 for cube in iter_maybe(cubes):
1096 # next 2 lines replace chunk of code.
1097 xcoord = _select_series_coord(cube, series_coordinate)
1098 xname = xcoord.points
1100 yfield = cube.data # power spectrum
1102 # If data from power spectra is all np.nans (like T+0h rainfall field which
1103 # might be full of zeros), then set yfield to zeros so it doesn't crash the
1104 # plotting.
1105 if np.all(np.isnan(yfield)):
1106 yfield = np.zeros_like(yfield)
1108 label = None
1109 color = "black"
1110 if model_colors_map: 1110 ↛ 1113line 1110 didn't jump to line 1113 because the condition on line 1110 was always true
1111 label = cube.attributes.get("model_name")
1112 color = model_colors_map.get(label)
1113 for cube_slice in cube.slices_over(ensemble_coord):
1114 # Label with (control) if part of an ensemble or not otherwise.
1115 if cube_slice.coord(ensemble_coord).points == [0]: 1115 ↛ 1129line 1115 didn't jump to line 1129 because the condition on line 1115 was always true
1116 ax.plot(
1117 xname,
1118 yfield,
1119 color=color,
1120 marker=line_marker,
1121 ls="-",
1122 lw=line_width,
1123 label=f"{label} (control)"
1124 if len(cube.coord(ensemble_coord).points) > 1
1125 else label,
1126 )
1127 # Label with (perturbed) if part of an ensemble and not the control.
1128 else:
1129 ax.plot(
1130 xname,
1131 yfield,
1132 color=color,
1133 ls="-",
1134 lw=1.5,
1135 alpha=0.75,
1136 label=f"{label} (member)",
1137 )
1139 # Calculate the global min/max if multiple cubes are given.
1140 _, levels, _ = colorbar_map_levels(cube, axis="y")
1141 if levels is not None: 1141 ↛ 1142line 1141 didn't jump to line 1142 because the condition on line 1141 was never true
1142 y_levels.append(min(levels))
1143 y_levels.append(max(levels))
1145 # Add some labels and tweak the style.
1147 title = f"{title}"
1148 ax.set_title(title, fontsize=16)
1150 # Set appropriate x-axis label based on coordinate
1151 if series_coordinate == "wavelength" or ( 1151 ↛ 1154line 1151 didn't jump to line 1154 because the condition on line 1151 was never true
1152 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1153 ):
1154 ax.set_xlabel("Wavelength (km)", fontsize=14)
1155 elif series_coordinate == "physical_wavenumber" or ( 1155 ↛ 1158line 1155 didn't jump to line 1158 because the condition on line 1155 was never true
1156 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1157 ):
1158 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1159 else: # frequency or check units
1160 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1160 ↛ 1161line 1160 didn't jump to line 1161 because the condition on line 1160 was never true
1161 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1162 else:
1163 ax.set_xlabel("Wavenumber", fontsize=14)
1165 ax.set_ylabel("Power Spectral Density", fontsize=14)
1166 ax.tick_params(axis="both", labelsize=12)
1168 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1170 # Set log-log scale
1171 ax.set_xscale("log")
1172 ax.set_yscale("log")
1174 # Add gridlines
1175 ax.grid(linestyle="--", color="grey", linewidth=1)
1176 # Ientify unique labels for legend
1177 handles = list(
1178 {
1179 label: handle
1180 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1181 }.values()
1182 )
1183 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1185 # Save plot.
1186 _save_close_figure(fig, "line power spectrum", filename)
1189def _plot_and_save_vertical_line_series(
1190 cubes: iris.cube.CubeList,
1191 coords: list[iris.coords.Coord],
1192 ensemble_coord: str,
1193 filename: str,
1194 series_coordinate: str,
1195 title: str,
1196 vmin: float,
1197 vmax: float,
1198 **kwargs,
1199):
1200 """Plot and save a 1D line series in vertical.
1202 Parameters
1203 ----------
1204 cubes: CubeList
1205 1 dimensional Cube or CubeList of the data to plot on x-axis.
1206 coord: list[Coord]
1207 Coordinates to plot on the y-axis, one per cube.
1208 ensemble_coord: str
1209 Ensemble coordinate in the cube.
1210 filename: str
1211 Filename of the plot to write.
1212 series_coordinate: str
1213 Coordinate to use as vertical axis.
1214 title: str
1215 Plot title.
1216 vmin: float
1217 Minimum value for the x-axis.
1218 vmax: float
1219 Maximum value for the x-axis.
1220 """
1221 # plot the vertical pressure axis using log scale
1222 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1224 model_colors_map = get_model_colors_map(cubes)
1226 # Check match-up across sequence coords gives consistent sizes
1227 validate_cubes_coords(cubes, coords)
1229 for cube, coord in zip(cubes, coords, strict=True):
1230 label = None
1231 color = "black"
1232 if model_colors_map: 1232 ↛ 1233line 1232 didn't jump to line 1233 because the condition on line 1232 was never true
1233 label = cube.attributes.get("model_name")
1234 color = model_colors_map.get(label)
1236 for cube_slice in cube.slices_over(ensemble_coord):
1237 # If ensemble data given plot control member with (control)
1238 # unless single forecast.
1239 if cube_slice.coord(ensemble_coord).points == [0]:
1240 iplt.plot(
1241 cube_slice,
1242 coord,
1243 color=color,
1244 marker="o",
1245 ls="-",
1246 lw=3,
1247 label=f"{label} (control)"
1248 if len(cube.coord(ensemble_coord).points) > 1
1249 else label,
1250 )
1251 # If ensemble data given plot perturbed members with (perturbed).
1252 else:
1253 iplt.plot(
1254 cube_slice,
1255 coord,
1256 color=color,
1257 ls="-",
1258 lw=1.5,
1259 alpha=0.75,
1260 label=f"{label} (member)",
1261 )
1263 # Get the current axis
1264 ax = plt.gca()
1266 # Special handling for pressure level data.
1267 if series_coordinate == "pressure": 1267 ↛ 1289line 1267 didn't jump to line 1289 because the condition on line 1267 was always true
1268 # Invert y-axis and set to log scale.
1269 ax.invert_yaxis()
1270 ax.set_yscale("log")
1272 # Define y-ticks and labels for pressure log axis.
1273 y_tick_labels = [
1274 "1000",
1275 "850",
1276 "700",
1277 "500",
1278 "300",
1279 "200",
1280 "100",
1281 ]
1282 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1284 # Set y-axis limits and ticks.
1285 ax.set_ylim(1100, 100)
1287 # Test if series_coordinate is model level data. The UM data uses
1288 # model_level_number and lfric uses full_levels as coordinate.
1289 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1290 # Define y-ticks and labels for vertical axis.
1291 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1292 y_tick_labels = [str(int(i)) for i in y_ticks]
1293 ax.set_ylim(min(y_ticks), max(y_ticks))
1295 ax.set_yticks(y_ticks)
1296 ax.set_yticklabels(y_tick_labels)
1298 # Set x-axis limits.
1299 ax.set_xlim(vmin, vmax)
1300 # Mark y=0 if present in plot.
1301 if vmin < 0.0 and vmax > 0.0: 1301 ↛ 1302line 1301 didn't jump to line 1302 because the condition on line 1301 was never true
1302 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1304 # Add some labels and tweak the style.
1305 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1306 ax.set_xlabel(
1307 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1308 )
1309 ax.set_title(title, fontsize=16)
1310 ax.ticklabel_format(axis="x")
1311 ax.tick_params(axis="y")
1312 ax.tick_params(axis="both", labelsize=12)
1314 # Add gridlines
1315 ax.grid(linestyle="--", color="grey", linewidth=1)
1316 # Ientify unique labels for legend
1317 handles = list(
1318 {
1319 label: handle
1320 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1321 }.values()
1322 )
1323 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1325 # Save plot.
1326 _save_close_figure(fig, "vertical line", filename)
1329def _plot_and_save_scatter_plot(
1330 cube_x: iris.cube.Cube | iris.cube.CubeList,
1331 cube_y: iris.cube.Cube | iris.cube.CubeList,
1332 filename: str,
1333 title: str,
1334 one_to_one: bool,
1335 model_names: list[str] | None = None,
1336 **kwargs,
1337):
1338 """Plot and save a 2D scatter plot.
1340 Parameters
1341 ----------
1342 cube_x: Cube | CubeList
1343 1 dimensional Cube or CubeList of the data to plot on x-axis.
1344 cube_y: Cube | CubeList
1345 1 dimensional Cube or CubeList of the data to plot on y-axis.
1346 filename: str
1347 Filename of the plot to write.
1348 title: str
1349 Plot title.
1350 one_to_one: bool
1351 Whether a 1:1 line is plotted.
1352 """
1353 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1354 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1355 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1356 # over the pairs simultaneously.
1358 # Ensure cube_x and cube_y are iterable
1359 cube_x_iterable = iter_maybe(cube_x)
1360 cube_y_iterable = iter_maybe(cube_y)
1362 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1363 iplt.scatter(cube_x_iter, cube_y_iter)
1364 if one_to_one is True:
1365 plt.plot(
1366 [
1367 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1368 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1369 ],
1370 [
1371 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1372 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1373 ],
1374 "k",
1375 linestyle="--",
1376 )
1377 ax = plt.gca()
1379 # Add some labels and tweak the style.
1380 if model_names is None:
1381 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1382 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1383 else:
1384 # Add the model names, these should be order of base (x) and other (y).
1385 ax.set_xlabel(
1386 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1387 )
1388 ax.set_ylabel(
1389 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1390 )
1391 ax.set_title(title, fontsize=16)
1392 ax.ticklabel_format(axis="y", useOffset=False)
1393 ax.tick_params(axis="x", labelrotation=15)
1394 ax.tick_params(axis="both", labelsize=12)
1395 ax.autoscale()
1397 # Save plot.
1398 _save_close_figure(fig, "scatter", filename)
1401def _plot_and_save_vector_plot(
1402 cube_u: iris.cube.Cube,
1403 cube_v: iris.cube.Cube,
1404 filename: str,
1405 title: str,
1406 method: Literal["contourf", "pcolormesh"],
1407 **kwargs,
1408):
1409 """Plot and save a 2D vector plot.
1411 Parameters
1412 ----------
1413 cube_u: Cube
1414 2 dimensional Cube of u component of the data.
1415 cube_v: Cube
1416 2 dimensional Cube of v component of the data.
1417 filename: str
1418 Filename of the plot to write.
1419 title: str
1420 Plot title.
1421 """
1422 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1423 # Create a cube containing the magnitude of the vector field.
1424 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1425 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1426 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1427 cube_vec_mag.rename(
1428 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1429 )
1431 # Specify the color bar
1432 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1434 # Setup plot map projection, extent and coastlines and borderlines.
1435 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1437 if method == "contourf":
1438 # Filled contour plot of the field.
1439 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1440 elif method == "pcolormesh":
1441 try:
1442 vmin = min(levels)
1443 vmax = max(levels)
1444 except TypeError:
1445 vmin, vmax = None, None
1446 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1447 # if levels are defined.
1448 if norm is not None:
1449 vmin = None
1450 vmax = None
1451 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1452 else:
1453 raise ValueError(f"Unknown plotting method: {method}")
1455 # Check to see if transect, and if so, adjust y axis.
1456 if is_transect(cube_vec_mag):
1457 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1458 axes.invert_yaxis()
1459 axes.set_yscale("log")
1460 axes.set_ylim(1100, 100)
1461 # If both model_level_number and level_height exists, iplt can construct
1462 # plot as a function of height above orography (NOT sea level).
1463 elif {"model_level_number", "level_height"}.issubset(
1464 {coord.name() for coord in cube_vec_mag.coords()}
1465 ):
1466 axes.set_yscale("log")
1468 axes.set_title(
1469 f"{title}\n"
1470 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1471 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1472 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1473 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1474 fontsize=16,
1475 )
1477 else:
1478 # Add title.
1479 axes.set_title(title, fontsize=16)
1481 # Add watermark with min/max/mean. Currently not user togglable.
1482 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1483 cube_min, cube_max, cube_mean = calc_array_stats(cube_vec_mag.data)
1484 axes.annotate(
1485 f"Min: {cube_min:.3g} Max: {cube_max:.3g} Mean: {cube_mean:.3g}",
1486 xy=(0.05, -0.05),
1487 xycoords="axes fraction",
1488 xytext=(-5, 5),
1489 textcoords="offset points",
1490 ha="right",
1491 va="bottom",
1492 size=11,
1493 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1494 )
1496 # Add colour bar.
1497 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1498 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1499 # add ticks and tick_labels for every levels if less than 20 levels exist
1500 if levels is not None and len(levels) < 20:
1501 cbar.set_ticks(levels)
1502 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1504 # 30 barbs along the longest axis of the plot, or a barb per point for data
1505 # with less than 30 points.
1506 step = max(max(cube_u.shape) // 30, 1)
1507 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1509 # Save plot.
1510 _save_close_figure(fig, "vector", filename)
1513def _plot_and_save_histogram_series(
1514 cubes: iris.cube.Cube | iris.cube.CubeList,
1515 filename: str,
1516 title: str,
1517 vmin: float,
1518 vmax: float,
1519 **kwargs,
1520):
1521 """Plot and save a histogram series.
1523 Parameters
1524 ----------
1525 cubes: Cube or CubeList
1526 2 dimensional Cube or CubeList of the data to plot as histogram.
1527 filename: str
1528 Filename of the plot to write.
1529 title: str
1530 Plot title.
1531 vmin: float
1532 minimum for colorbar
1533 vmax: float
1534 maximum for colorbar
1535 """
1536 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1537 ax = plt.gca()
1539 model_colors_map = get_model_colors_map(cubes)
1541 # Set default that histograms will produce probability density function
1542 # at each bin (integral over range sums to 1).
1543 density = True
1545 for cube in iter_maybe(cubes):
1546 # Easier to check title (where var name originates)
1547 # than seeing if long names exist etc.
1548 # Exception case, where distribution better fits log scales/bins.
1549 if (
1550 ("surface_microphysical" in title)
1551 or ("rain accumulation" in title)
1552 or ("Rainfall rate Composite" in title)
1553 or ("Nimrod_5min" in title)
1554 ):
1555 if "amount" in title:
1556 # Compute histogram following Klingaman et al. (2017): ASoP
1557 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1558 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1559 density = False
1560 else:
1561 bins = 10.0 ** (
1562 np.arange(-10, 27, 1) / 10.0
1563 ) # Suggestion from RMED toolbox.
1564 bins = np.insert(bins, 0, 0)
1565 ax.set_yscale("log")
1566 vmin = bins[1]
1567 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1568 ax.set_xscale("log")
1569 elif "lightning" in title:
1570 bins = [0, 1, 2, 3, 4, 5]
1571 else:
1572 bins = np.linspace(vmin, vmax, 51)
1573 logger.debug(
1574 "Plotting histogram with %s bins %s - %s.",
1575 np.size(bins),
1576 np.min(bins),
1577 np.max(bins),
1578 )
1580 # Reshape cube data into a single array to allow for a single histogram.
1581 # Otherwise we plot xdim histograms stacked.
1582 cube_data_1d = (cube.data).flatten()
1584 label = None
1585 color = "black"
1586 if model_colors_map:
1587 label = cube.attributes.get("model_name")
1588 color = model_colors_map[label]
1589 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1591 # Compute area under curve.
1592 if (
1593 ("surface_microphysical" in title and "amount" in title)
1594 or ("rain_accumulation" in title)
1595 or ("Rainfall rate Composite" in title)
1596 or ("Nimrod_5min" in title)
1597 ):
1598 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1599 x = x * bin_mean / x.sum()
1600 x = x[1:]
1601 y = y[1:]
1603 ax.plot(
1604 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1605 )
1607 # Add some labels and tweak the style.
1608 ax.set_title(title, fontsize=16)
1609 ax.set_xlabel(
1610 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1611 )
1612 ax.set_ylabel("Normalised probability density", fontsize=14)
1613 if (
1614 ("surface_microphysical" in title and "amount" in title)
1615 or ("rain accumulation" in title)
1616 or ("Nimrod_5min" in title)
1617 ):
1618 ax.set_ylabel(
1619 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1620 )
1621 try:
1622 ax.set_xlim(vmin, vmax)
1623 except ValueError:
1624 pass
1625 ax.tick_params(axis="both", labelsize=12)
1627 # Overlay grid-lines onto histogram plot.
1628 ax.grid(linestyle="--", color="grey", linewidth=1)
1629 if model_colors_map:
1630 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1632 # Save plot.
1633 _save_close_figure(fig, "histogram", filename)
1636def _plot_and_save_postage_stamp_histogram_series(
1637 cube: iris.cube.Cube,
1638 filename: str,
1639 title: str,
1640 stamp_coordinate: str,
1641 vmin: float,
1642 vmax: float,
1643 **kwargs,
1644):
1645 """Plot and save postage (ensemble members) stamps for a histogram series.
1647 Parameters
1648 ----------
1649 cube: Cube
1650 2 dimensional Cube of the data to plot as histogram.
1651 filename: str
1652 Filename of the plot to write.
1653 title: str
1654 Plot title.
1655 stamp_coordinate: str
1656 Coordinate that becomes different plots.
1657 vmin: float
1658 minimum for pdf x-axis
1659 vmax: float
1660 maximum for pdf x-axis
1661 """
1662 # Use the smallest square grid that will fit the members.
1663 nmember = len(cube.coord(stamp_coordinate).points)
1664 grid_rows = int(math.sqrt(nmember))
1665 grid_size = math.ceil(nmember / grid_rows)
1667 fig = plt.figure(
1668 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1669 )
1670 # Make a subplot for each member.
1671 for member, subplot in zip(
1672 cube.slices_over(stamp_coordinate),
1673 range(1, grid_size * grid_rows + 1),
1674 strict=False,
1675 ):
1676 # Implicit interface is much easier here, due to needing to have the
1677 # cartopy GeoAxes generated.
1678 plt.subplot(grid_rows, grid_size, subplot)
1679 # Reshape cube data into a single array to allow for a single histogram.
1680 # Otherwise we plot xdim histograms stacked.
1681 member_data_1d = (member.data).flatten()
1682 plt.hist(member_data_1d, density=True, stacked=True)
1683 axes = plt.gca()
1684 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1685 axes.set_title(f"{mtitle}")
1686 axes.set_xlim(vmin, vmax)
1688 # Overall figure title.
1689 fig.suptitle(title, fontsize=16)
1691 # Save plot.
1692 _save_close_figure(fig, "histogram postage stamp", filename)
1695def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1696 cube: iris.cube.Cube,
1697 filename: str,
1698 title: str,
1699 stamp_coordinate: str,
1700 vmin: float,
1701 vmax: float,
1702 **kwargs,
1703):
1704 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1705 ax.set_title(title, fontsize=16)
1706 ax.set_xlim(vmin, vmax)
1707 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1708 ax.set_ylabel("normalised probability density", fontsize=14)
1709 # Loop over all slices along the stamp_coordinate
1710 for member in cube.slices_over(stamp_coordinate):
1711 # Flatten the member data to 1D
1712 member_data_1d = member.data.flatten()
1713 # Plot the histogram using plt.hist
1714 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1715 plt.hist(
1716 member_data_1d,
1717 density=True,
1718 stacked=True,
1719 label=f"{mtitle}",
1720 )
1722 # Add a legend
1723 ax.legend(fontsize=16)
1725 # Save plot.
1726 _save_close_figure(fig, "histogram postage stamp", filename)
1729def _plot_and_save_scatter_series(
1730 cubes: iris.cube.Cube | iris.cube.CubeList,
1731 filename: str,
1732 title: str,
1733 vmin: float,
1734 vmax: float,
1735 hexbin: bool,
1736 **kwargs,
1737):
1738 """Plot and save a scatter plot series.
1740 Parameters
1741 ----------
1742 cubes: Cube or CubeList
1743 2 dimensional Cube or CubeList of the data to plot as scatter.
1744 filename: str
1745 Filename of the plot to write.
1746 title: str
1747 Plot title.
1748 vmin: float
1749 minimum for colorbar
1750 vmax: float
1751 maximum for colorbar
1752 hexbin: bool
1753 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1754 Else scatter of all points, with potential to overplot many comparisons on same plot.
1755 """
1756 if hexbin:
1757 # Check cubes using same functionality as the difference operator.
1758 if len(cubes) != 2:
1759 raise ValueError(
1760 "Cubes should contain exactly 2 cubes for hexbin plotting."
1761 )
1762 title = title.replace("scatter", "hexbin")
1763 filename = filename.replace("scatter", "hexbin")
1765 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1766 ax = plt.gca()
1768 model_colors_map = get_model_colors_map(cubes)
1770 percentiles = np.arange(0, 100, 5)
1771 percentiles[0] = 1
1772 percentiles[-1] = 99
1773 quantiles = iris.cube.CubeList()
1775 # Loop through all output cubes for both data points and overplotting quantiles.
1776 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1777 for plottype in ["points", "quantiles"]:
1778 nplot = 0
1779 for cube in iter_maybe(cubes):
1780 label = None
1781 color = "black"
1782 if model_colors_map: 1782 ↛ 1787line 1782 didn't jump to line 1787 because the condition on line 1782 was always true
1783 label = cube.attributes.get("model_name")
1784 color = model_colors_map[label]
1786 # Plot all data points
1787 if plottype == "points":
1788 if nplot > 0:
1789 if hexbin:
1790 hb = plt.hexbin(
1791 cubes[0].data.flatten(),
1792 cube.data.flatten(),
1793 alpha=0.3,
1794 gridsize=100,
1795 mincnt=1,
1796 )
1797 else:
1798 plt.scatter(
1799 cubes[0].data.flatten(),
1800 cube.data.flatten(),
1801 color=color,
1802 marker="+",
1803 label=None,
1804 alpha=0.3,
1805 )
1807 elif plottype == "quantiles": 1807 ↛ 1826line 1807 didn't jump to line 1826 because the condition on line 1807 was always true
1808 # Construct Q-Q plot
1809 quantiles.append(
1810 cube.collapsed(
1811 cube.coords(dim_coords=True),
1812 iris.analysis.PERCENTILE,
1813 percent=percentiles,
1814 )
1815 )
1816 if nplot > 0:
1817 iplt.scatter(
1818 quantiles[0],
1819 quantiles[-1],
1820 color=color,
1821 marker="o",
1822 label=label,
1823 edgecolors="black",
1824 )
1826 nplot = nplot + 1
1828 # Add some labels and tweak the style.
1829 ax.set_title(title, fontsize=16)
1830 ax.set_xlabel(
1831 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1832 )
1833 ax.set_ylabel(
1834 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1835 )
1836 ax.tick_params(axis="both", labelsize=12)
1837 ax.autoscale()
1839 # Set 1:1 line and equal axes if scatter plot of common cube names
1840 nameA = iter_maybe(cubes)[0].name()
1841 nameB = iter_maybe(cubes)[1].name()
1842 if any(part in nameB.split("_") for part in nameA.split("_")): 1842 ↛ 1853line 1842 didn't jump to line 1853 because the condition on line 1842 was always true
1843 lims = [
1844 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1845 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1846 ]
1847 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1848 ax.set_aspect("equal")
1849 ax.set_xlim(lims)
1850 ax.set_ylim(lims)
1852 # Overlay grid-lines onto scatter plot.
1853 ax.grid(linestyle="--", color="grey", linewidth=1)
1854 if model_colors_map: 1854 ↛ 1858line 1854 didn't jump to line 1858 because the condition on line 1854 was always true
1855 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1857 # Add colorbar if hexbin output
1858 if hexbin:
1859 cb = plt.colorbar(
1860 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1861 )
1862 cb.set_label("Number of data points", size=12)
1864 # Save plot.
1865 _save_close_figure(fig, "scatter", filename)
1868def _spatial_plot(
1869 method: Literal["contourf", "pcolormesh", "scatter"],
1870 cube: iris.cube.Cube,
1871 filename: str | None,
1872 sequence_coordinate: str,
1873 stamp_coordinate: str,
1874 strict_title: bool = False,
1875 overlay_cube: iris.cube.Cube | None = None,
1876 contour_cube: iris.cube.Cube | None = None,
1877 point_cube: iris.cube.Cube | None = None,
1878 **kwargs,
1879):
1880 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1882 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1883 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1884 is present then postage stamp plots will be produced.
1886 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1887 be overplotted on the same figure.
1889 Parameters
1890 ----------
1891 method: "contourf" | "pcolormesh" | "scatter"
1892 The plotting method to use.
1893 Select choice of "contourf" or "pcolormesh" for gridded data.
1894 Use "scatter" for point-based data.
1895 cube: Cube
1896 Iris cube of the data to plot. It should have two spatial dimensions,
1897 such as lat and lon, and may also have a another two dimension to be
1898 plotted sequentially and/or as postage stamp plots.
1899 filename: str | None
1900 Name of the plot to write, used as a prefix for plot sequences. If None
1901 uses the recipe name.
1902 sequence_coordinate: str
1903 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1904 This coordinate must exist in the cube.
1905 stamp_coordinate: str
1906 Coordinate about which to plot postage stamp plots. Defaults to
1907 ``"realization"``.
1908 strict_title: bool, optional
1909 Logical switch that if set to True will ensure that the MODEL_NAME
1910 string is not prepended to the plot title. The default is False.
1911 overlay_cube: Cube | None, optional
1912 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1913 contour_cube: Cube | None, optional
1914 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1915 point_cube: Cube | None, optional
1916 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
1918 Raises
1919 ------
1920 ValueError
1921 If the cube doesn't have the right dimensions.
1922 TypeError
1923 If the cube isn't a single cube.
1924 """
1925 # Ensure we've got a single cube.
1926 cube = check_single_cube(cube)
1928 # Set title based on recipe metadata or use cube name
1929 recipe_title = get_recipe_metadata().get("title", cube.name())
1931 # Check if there is a valid stamp coordinate in cube dimensions.
1932 if stamp_coordinate == "realization": 1932 ↛ 1937line 1932 didn't jump to line 1937 because the condition on line 1932 was always true
1933 stamp_coordinate = check_stamp_coordinate(cube)
1935 # Make postage stamp plots if stamp_coordinate exists and has more than a
1936 # single point.
1937 plotting_func = _plot_and_save_spatial_plot
1938 try:
1939 if cube.coord(stamp_coordinate).shape[0] > 1:
1940 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1941 except iris.exceptions.CoordinateNotFoundError:
1942 pass
1944 # Produce a geographical scatter plot if the data have a
1945 # dimension called observation or model_obs_error
1946 if any(
1947 crd.var_name == "station"
1948 or crd.var_name == "Station_Name"
1949 or crd.var_name == "model_obs_error"
1950 for crd in cube.coords()
1951 ):
1952 plotting_func = _plot_and_save_spatial_plot
1953 method = "scatter"
1955 # Must have a sequence coordinate.
1956 try:
1957 cube.coord(sequence_coordinate)
1958 except iris.exceptions.CoordinateNotFoundError as err:
1959 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1961 # Create a plot for each value of the sequence coordinate.
1962 plot_index = []
1963 nplot = np.size(cube.coord(sequence_coordinate).points)
1965 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1966 # Set plot titles and filename
1967 seq_coord = cube_slice.coord(sequence_coordinate)
1969 if "model_name" in cube.attributes: 1969 ↛ 1970line 1969 didn't jump to line 1970 because the condition on line 1969 was never true
1970 model_name = cube.attributes["model_name"]
1971 else:
1972 model_name = None
1973 if strict_title: 1973 ↛ 1974line 1973 didn't jump to line 1974 because the condition on line 1973 was never true
1974 model_name = None
1976 plot_title, plot_filename = _set_title_and_filename(
1977 seq_coord, nplot, recipe_title, filename, model_name=model_name
1978 )
1980 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1981 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1982 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1983 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1985 # Do the actual plotting.
1986 plotting_func(
1987 cube_slice,
1988 filename=plot_filename,
1989 stamp_coordinate=stamp_coordinate,
1990 title=plot_title,
1991 method=method,
1992 overlay_cube=overlay_slice,
1993 contour_cube=contour_slice,
1994 point_cube=point_slice,
1995 **kwargs,
1996 )
1997 plot_index.append(plot_filename)
1999 # Add list of plots to plot metadata.
2000 complete_plot_index = _append_to_plot_index(plot_index)
2002 # Make a page to display the plots.
2003 _make_plot_html_page(complete_plot_index)
2006####################
2007# Public functions #
2008####################
2011def spatial_contour_plot(
2012 cube: iris.cube.Cube,
2013 filename: str | None = None,
2014 sequence_coordinate: str = "time",
2015 stamp_coordinate: str = "realization",
2016 **kwargs,
2017) -> iris.cube.Cube:
2018 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2020 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2021 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2022 is present then postage stamp plots will be produced.
2024 Parameters
2025 ----------
2026 cube: Cube
2027 Iris cube of the data to plot. It should have two spatial dimensions,
2028 such as lat and lon, and may also have a another two dimension to be
2029 plotted sequentially and/or as postage stamp plots.
2030 filename: str, optional
2031 Name of the plot to write, used as a prefix for plot sequences. Defaults
2032 to the recipe name.
2033 sequence_coordinate: str, optional
2034 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2035 This coordinate must exist in the cube.
2036 stamp_coordinate: str, optional
2037 Coordinate about which to plot postage stamp plots. Defaults to
2038 ``"realization"``.
2040 Returns
2041 -------
2042 Cube
2043 The original cube (so further operations can be applied).
2045 Raises
2046 ------
2047 ValueError
2048 If the cube doesn't have the right dimensions.
2049 TypeError
2050 If the cube isn't a single cube.
2051 """
2052 _spatial_plot(
2053 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2054 )
2055 return cube
2058def spatial_pcolormesh_plot(
2059 cubes: iris.cube.Cube | iris.cube.CubeList,
2060 filename: str | None = None,
2061 sequence_coordinate: str = "time",
2062 stamp_coordinate: str = "realization",
2063 strict_title: bool = False,
2064 **kwargs,
2065) -> iris.cube.Cube:
2066 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2068 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2069 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2070 is present then postage stamp plots will be produced.
2072 This function is significantly faster than ``spatial_contour_plot``,
2073 especially at high resolutions, and should be preferred unless contiguous
2074 contour areas are important.
2076 Parameters
2077 ----------
2078 cube: Cubes
2079 Iris cube or cubelist of the data to plot. Each cube should have two spatial dimensions,
2080 such as lat and lon, and may also have a another two dimension to be
2081 plotted sequentially and/or as postage stamp plots.
2082 filename: str, optional
2083 Name of the plot to write, used as a prefix for plot sequences. Defaults
2084 to the recipe name.
2085 sequence_coordinate: str, optional
2086 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2087 This coordinate must exist in the cube.
2088 stamp_coordinate: str, optional
2089 Coordinate about which to plot postage stamp plots. Defaults to
2090 ``"realization"``.
2091 strict_title: bool, optional
2092 Logical switch that if set to True will ensure that the MODEL_NAME
2093 string is not prepended to the plot title. The default is False.
2095 Returns
2096 -------
2097 Cubes
2098 The original cube/cubelist (so further operations can be applied).
2100 Raises
2101 ------
2102 ValueError
2103 If the cube doesn't have the right dimensions.
2104 """
2105 if isinstance(cubes, iris.cube.CubeList): 2105 ↛ 2106line 2105 didn't jump to line 2106 because the condition on line 2105 was never true
2106 for model_cube in cubes:
2107 _spatial_plot(
2108 "pcolormesh",
2109 model_cube,
2110 filename,
2111 sequence_coordinate,
2112 stamp_coordinate,
2113 strict_title,
2114 **kwargs,
2115 )
2116 elif isinstance(cubes, iris.cube.Cube): 2116 ↛ 2126line 2116 didn't jump to line 2126 because the condition on line 2116 was always true
2117 _spatial_plot(
2118 "pcolormesh",
2119 cubes,
2120 filename,
2121 sequence_coordinate,
2122 stamp_coordinate,
2123 strict_title,
2124 **kwargs,
2125 )
2126 return cubes
2129def spatial_multi_pcolormesh_plot(
2130 cube: iris.cube.Cube,
2131 overlay_cube: iris.cube.Cube | None = None,
2132 contour_cube: iris.cube.Cube | None = None,
2133 point_cube: iris.cube.Cube | None = None,
2134 filename: str | None = None,
2135 sequence_coordinate: str = "time",
2136 stamp_coordinate: str = "realization",
2137 **kwargs,
2138) -> iris.cube.Cube:
2139 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2141 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2142 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2143 is present then postage stamp plots will be produced.
2145 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2147 If specified, contours of a contour_cube can be overplotted on top of those.
2149 If specified, a spatial scatter map of point_cube can be overplotted.
2151 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2153 This function is significantly faster than ``spatial_contour_plot``,
2154 especially at high resolutions, and should be preferred unless contiguous
2155 contour areas are important.
2157 Parameters
2158 ----------
2159 cube: Cube
2160 Iris cube of the data to plot. It should have two spatial dimensions,
2161 such as lat and lon, and may also have two additional dimensions to be
2162 plotted sequentially and/or as postage stamp plots.
2163 overlay_cube: Cube, optional
2164 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2165 such as lat and lon, and may also have two additional dimensions to be
2166 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.
2167 If not provided, output plot generated without overlay cube.
2168 contour_cube: Cube, optional
2169 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,
2170 such as lat and lon, and may also have two additional dimensions to be
2171 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2172 point_cube: Cube, optional
2173 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
2174 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2175 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2176 postage stamp plots. If not provided, output plot generated without point-based layer.
2177 filename: str, optional
2178 Name of the plot to write, used as a prefix for plot sequences. Defaults
2179 to the recipe name.
2180 sequence_coordinate: str, optional
2181 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2182 This coordinate must exist in the cube.
2183 stamp_coordinate: str, optional
2184 Coordinate about which to plot postage stamp plots. Defaults to
2185 ``"realization"``.
2187 Returns
2188 -------
2189 Cube
2190 The original cube (so further operations can be applied).
2192 Raises
2193 ------
2194 ValueError
2195 If the cube doesn't have the right dimensions.
2196 TypeError
2197 If the cube isn't a single cube.
2198 """
2199 _spatial_plot(
2200 "pcolormesh",
2201 cube,
2202 filename,
2203 sequence_coordinate,
2204 stamp_coordinate,
2205 overlay_cube=overlay_cube,
2206 contour_cube=contour_cube,
2207 point_cube=point_cube,
2208 )
2209 return cube, overlay_cube, contour_cube, point_cube
2212# TODO: Expand function to handle ensemble data.
2213# line_coordinate: str, optional
2214# Coordinate about which to plot multiple lines. Defaults to
2215# ``"realization"``.
2216def plot_line_series(
2217 cube: iris.cube.Cube | iris.cube.CubeList,
2218 filename: str | None = None,
2219 series_coordinate: str = "time",
2220 sequence_coordinate: str = "time",
2221 # add the following for ensembles
2222 stamp_coordinate: str = "realization",
2223 single_plot: bool = False,
2224 **kwargs,
2225) -> iris.cube.Cube | iris.cube.CubeList:
2226 """Plot a line plot for the specified coordinate.
2228 The Cube or CubeList must be 1D.
2230 Parameters
2231 ----------
2232 iris.cube | iris.cube.CubeList
2233 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2234 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2235 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2236 filename: str, optional
2237 Name of the plot to write, used as a prefix for plot sequences. Defaults
2238 to the recipe name.
2239 series_coordinate: str, optional
2240 Coordinate about which to make a series. Defaults to ``"time"``. This
2241 coordinate must exist in the cube.
2243 Returns
2244 -------
2245 iris.cube.Cube | iris.cube.CubeList
2246 The original Cube or CubeList (so further operations can be applied).
2248 Raises
2249 ------
2250 ValueError
2251 If the cubes don't have the right dimensions.
2252 TypeError
2253 If the cube isn't a Cube or CubeList.
2254 """
2255 # Ensure we have a name for the plot file.
2256 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2258 num_models = get_num_models(cube)
2260 validate_cube_shape(cube, num_models)
2262 # Iterate over all cubes and extract coordinate to plot.
2263 cubes = iris.cube.CubeList(iter_maybe(cube))
2265 coords = []
2266 for model_cube in cubes:
2267 try:
2268 coords.append(model_cube.coord(series_coordinate))
2269 except iris.exceptions.CoordinateNotFoundError as err:
2270 raise ValueError(
2271 f"Cube must have a {series_coordinate} coordinate."
2272 ) from err
2273 # Count cube dimensions and exclude realization and
2274 # forecast_reference_time if they exist.
2275 ndim = model_cube.ndim
2277 if model_cube.coords("realization"): 2277 ↛ 2285line 2277 didn't jump to line 2285 because the condition on line 2277 was always true
2278 # returns coord dimension
2279 realization_dims = model_cube.coord_dims("realization")
2281 # Only subtract if realization is a dimension coordinate
2282 if realization_dims:
2283 ndim -= len(realization_dims)
2285 if model_cube.coords("forecast_reference_time"):
2286 frt_dims = model_cube.coord_dims("forecast_reference_time")
2288 # Only subtract if frt is a dimension coordinate
2289 if frt_dims: 2289 ↛ 2290line 2289 didn't jump to line 2290 because the condition on line 2289 was never true
2290 ndim -= len(frt_dims)
2292 if ndim > 2:
2293 raise ValueError(
2294 "Cube must be 1D or 2D (excluding any realization or forecast_reference_time dimensions)."
2295 )
2297 plot_index = []
2299 # Check if this is a spectral plot by looking for spectral coordinates
2300 is_spectral_plot = series_coordinate in [
2301 "frequency",
2302 "physical_wavenumber",
2303 "wavelength",
2304 ]
2306 if is_spectral_plot:
2307 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2308 # coordinate frequency/wavenumber.
2309 # If several power spectra are plotted with time as sequence_coordinate for the
2310 # time slider option.
2312 # Internal plotting function.
2313 plotting_func = _plot_and_save_line_power_spectrum_series
2315 for model_cube in cubes:
2316 try:
2317 model_cube.coord(sequence_coordinate)
2318 except iris.exceptions.CoordinateNotFoundError as err:
2319 raise ValueError(
2320 f"Cube must have a {sequence_coordinate} coordinate."
2321 ) from err
2323 if num_models == 1: 2323 ↛ 2338line 2323 didn't jump to line 2338 because the condition on line 2323 was always true
2324 # check for ensembles
2325 if ( 2325 ↛ 2329line 2325 didn't jump to line 2329 because the condition on line 2325 was never true
2326 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2327 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2328 ):
2329 if single_plot:
2330 # Plot spectra, mean and ensemble spread on 1 plot
2331 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2332 else:
2333 # Plot postage stamps
2334 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2335 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2336 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2337 else:
2338 all_points = sorted(
2339 set(
2340 itertools.chain.from_iterable(
2341 cb.coord(sequence_coordinate).points for cb in cubes
2342 )
2343 )
2344 )
2345 all_slices = list(
2346 itertools.chain.from_iterable(
2347 cb.slices_over(sequence_coordinate) for cb in cubes
2348 )
2349 )
2350 # Matched slices (matched by seq coord point; it may happen that
2351 # evaluated models do not cover the same seq coord range, hence matching
2352 # necessary)
2353 cube_iterables = [
2354 iris.cube.CubeList(
2355 s
2356 for s in all_slices
2357 if s.coord(sequence_coordinate).points[0] == point
2358 )
2359 for point in all_points
2360 ]
2361 nplot = len(all_points)
2363 # Create a plot for each value of the sequence coordinate. Allowing for
2364 # multiple cubes in a CubeList to be plotted in the same plot for similar
2365 # sequence values. Passing a CubeList into the internal plotting function
2366 # for similar values of the sequence coordinate. cube_slice can be an
2367 # iris.cube.Cube or an iris.cube.CubeList.
2369 for cube_slice in cube_iterables:
2370 # Normalize cube_slice to a list of cubes
2371 if isinstance(cube_slice, iris.cube.CubeList): 2371 ↛ 2372line 2371 didn't jump to line 2372 because the condition on line 2371 was never true
2372 cubes = list(cube_slice)
2373 elif isinstance(cube_slice, iris.cube.Cube): 2373 ↛ 2376line 2373 didn't jump to line 2376 because the condition on line 2373 was always true
2374 cubes = [cube_slice]
2375 else:
2376 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2378 # Use sequence value so multiple sequences can merge.
2379 seq_coord = cube_slice[0].coord(sequence_coordinate)
2380 plot_title, plot_filename = _set_title_and_filename(
2381 seq_coord, nplot, recipe_title, filename
2382 )
2384 # Format the coordinate value in a unit appropriate way.
2385 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2387 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2388 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2388 ↛ 2389line 2388 didn't jump to line 2389 because the condition on line 2388 was never true
2389 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2391 # Do the actual plotting.
2392 plotting_func(
2393 cube_slice,
2394 coords,
2395 stamp_coordinate,
2396 plot_filename,
2397 title,
2398 series_coordinate,
2399 )
2401 plot_index.append(plot_filename)
2402 else:
2403 # Format the title and filename using plotted series coordinate
2404 nplot = 1
2405 seq_coord = coords[0]
2406 plot_title, plot_filename = _set_title_and_filename(
2407 seq_coord, nplot, recipe_title, filename
2408 )
2410 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2411 if (
2412 "station" in [c.name() for c in cubes[0].coords()]
2413 and len(cubes[0].coord("station").points) > 1
2414 ):
2415 for station in cubes[0].coord("station").points:
2416 station_cubes = cubes.extract(iris.Constraint(station=station))
2417 station_name = station_cubes[0].coord("Station_Name").points[0]
2418 station_plotname = plot_filename.replace(
2419 ".png", "_" + station_name + ".png"
2420 )
2421 _plot_and_save_line_series(
2422 station_cubes,
2423 coords,
2424 "realization",
2425 station_plotname,
2426 f"{plot_title} {station_name}",
2427 )
2428 plot_index.append(station_plotname)
2430 else:
2431 # Do the actual plotting for all other series coordinate options.
2432 _plot_and_save_line_series(
2433 cubes, coords, stamp_coordinate, plot_filename, plot_title
2434 )
2436 plot_index.append(plot_filename)
2438 # append plot to list of plots
2439 complete_plot_index = _append_to_plot_index(plot_index)
2441 # Make a page to display the plots.
2442 _make_plot_html_page(complete_plot_index)
2444 return cube
2447def plot_vertical_line_series(
2448 cubes: iris.cube.Cube | iris.cube.CubeList,
2449 filename: str | None = None,
2450 series_coordinate: str = "model_level_number",
2451 sequence_coordinate: str = "time",
2452 # line_coordinate: str = "realization",
2453 **kwargs,
2454) -> iris.cube.Cube | iris.cube.CubeList:
2455 """Plot a line plot against a type of vertical coordinate.
2457 The Cube or CubeList must be 1D.
2459 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2460 then a sequence of plots will be produced.
2462 Parameters
2463 ----------
2464 iris.cube | iris.cube.CubeList
2465 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2466 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2467 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2468 filename: str, optional
2469 Name of the plot to write, used as a prefix for plot sequences. Defaults
2470 to the recipe name.
2471 series_coordinate: str, optional
2472 Coordinate to plot on the y-axis. Can be ``pressure`` or
2473 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2474 for LFRic. Defaults to ``model_level_number``.
2475 This coordinate must exist in the cube.
2476 sequence_coordinate: str, optional
2477 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2478 This coordinate must exist in the cube.
2480 Returns
2481 -------
2482 iris.cube.Cube | iris.cube.CubeList
2483 The original Cube or CubeList (so further operations can be applied).
2484 Plotted data.
2486 Raises
2487 ------
2488 ValueError
2489 If the cubes doesn't have the right dimensions.
2490 TypeError
2491 If the cube isn't a Cube or CubeList.
2492 """
2493 # Ensure we have a name for the plot file.
2494 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2496 cubes = iter_maybe(cubes)
2497 # Initialise empty list to hold all data from all cubes in a CubeList
2498 all_data = []
2500 # Store min/max ranges for x range.
2501 x_levels = []
2503 num_models = get_num_models(cubes)
2505 validate_cube_shape(cubes, num_models)
2507 # Iterate over all cubes in cube or CubeList and plot.
2508 coords = []
2509 for cube in cubes:
2510 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2511 try:
2512 coords.append(cube.coord(series_coordinate))
2513 except iris.exceptions.CoordinateNotFoundError as err:
2514 raise ValueError(
2515 f"Cube must have a {series_coordinate} coordinate."
2516 ) from err
2518 try:
2519 if cube.ndim > 1 or not cube.coords("realization"): 2519 ↛ 2527line 2519 didn't jump to line 2527 because the condition on line 2519 was always true
2520 cube.coord(sequence_coordinate)
2521 except iris.exceptions.CoordinateNotFoundError as err:
2522 raise ValueError(
2523 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2524 ) from err
2526 # Get minimum and maximum from levels information.
2527 _, levels, _ = colorbar_map_levels(cube, axis="x")
2528 if levels is not None: 2528 ↛ 2532line 2528 didn't jump to line 2532 because the condition on line 2528 was always true
2529 x_levels.append(min(levels))
2530 x_levels.append(max(levels))
2531 else:
2532 all_data.append(cube.data)
2534 if len(x_levels) == 0: 2534 ↛ 2536line 2534 didn't jump to line 2536 because the condition on line 2534 was never true
2535 # Combine all data into a single NumPy array
2536 combined_data = np.concatenate(all_data)
2538 # Set the lower and upper limit for the x-axis to ensure all plots have
2539 # same range. This needs to read the whole cube over the range of the
2540 # sequence and if applicable postage stamp coordinate.
2541 vmin = np.floor(combined_data.min())
2542 vmax = np.ceil(combined_data.max())
2543 else:
2544 vmin = min(x_levels)
2545 vmax = max(x_levels)
2547 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2548 # a single profile directly without iterating over a sequence.
2549 sequence_coords = [
2550 cube.coord(sequence_coordinate)
2551 for cube in cubes
2552 if cube.coords(sequence_coordinate)
2553 ]
2554 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2555 np.size(coord.points) > 1 for coord in sequence_coords
2556 )
2557 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2558 np.size(coord.points) == 1 for coord in sequence_coords
2559 )
2561 plot_index = []
2562 if has_sequence_coord: 2562 ↛ 2587line 2562 didn't jump to line 2587 because the condition on line 2562 was always true
2563 # Matching the slices (matching by seq coord point; it may happen that
2564 # evaluated models do not cover the same seq coord range, hence matching
2565 # necessary)
2566 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2567 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2568 for cubes_slice in cube_iterables:
2569 # Format the coordinate value in a unit appropriate way.
2570 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2571 plot_title, plot_filename = _set_title_and_filename(
2572 seq_coord, nplot, recipe_title, filename
2573 )
2575 # Do the actual plotting.
2576 _plot_and_save_vertical_line_series(
2577 cubes_slice,
2578 coords,
2579 "realization",
2580 plot_filename,
2581 series_coordinate,
2582 title=plot_title,
2583 vmin=vmin,
2584 vmax=vmax,
2585 )
2586 plot_index.append(plot_filename)
2587 elif has_scalar_sequence_coord:
2588 # Scalar sequence coordinate (typically aggregated time bounds):
2589 # make one plot and include sequence period in title/filename.
2590 plot_title, plot_filename = _set_title_and_filename(
2591 sequence_coords[0], 1, recipe_title, filename
2592 )
2594 _plot_and_save_vertical_line_series(
2595 cubes,
2596 coords,
2597 "realization",
2598 plot_filename,
2599 series_coordinate,
2600 title=plot_title,
2601 vmin=vmin,
2602 vmax=vmax,
2603 )
2604 plot_index.append(plot_filename)
2605 else:
2606 # 1D case: no sequence coordinate, plot a single profile.
2607 plot_title = recipe_title
2608 if filename:
2609 plot_filename = filename
2610 else:
2611 plot_filename = f"{slugify(plot_title)}.png"
2613 _plot_and_save_vertical_line_series(
2614 cubes,
2615 coords,
2616 "realization",
2617 plot_filename,
2618 series_coordinate,
2619 title=plot_title,
2620 vmin=vmin,
2621 vmax=vmax,
2622 )
2623 plot_index.append(plot_filename)
2625 # Add list of plots to plot metadata.
2626 complete_plot_index = _append_to_plot_index(plot_index)
2628 # Make a page to display the plots.
2629 _make_plot_html_page(complete_plot_index)
2631 return cubes
2634def qq_plot(
2635 cubes: iris.cube.CubeList,
2636 coordinates: list[str],
2637 percentiles: list[float],
2638 model_names: list[str],
2639 filename: str | None = None,
2640 one_to_one: bool = True,
2641 **kwargs,
2642) -> iris.cube.CubeList:
2643 """Plot a Quantile-Quantile plot between two models for common time points.
2645 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2646 collapsed within the operator over all specified coordinates such as
2647 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2648 ensemble members to ensure a 1D cube (array).
2650 Parameters
2651 ----------
2652 cubes: iris.cube.CubeList
2653 Two cubes of the same variable with different models.
2654 coordinate: list[str]
2655 The list of coordinates to collapse over. This list should be
2656 every coordinate within the cube to result in a 1D cube around
2657 the percentile coordinate.
2658 percent: list[float]
2659 A list of percentiles to appear in the plot.
2660 model_names: list[str]
2661 A list of model names to appear on the axis of the plot.
2662 filename: str, optional
2663 Filename of the plot to write.
2664 one_to_one: bool, optional
2665 If True a 1:1 line is plotted; if False it is not. Default is True.
2667 Raises
2668 ------
2669 ValueError
2670 When the cubes are not compatible.
2672 Notes
2673 -----
2674 The quantile-quantile plot is a variant on the scatter plot representing
2675 two datasets by their quantiles (percentiles) for common time points.
2676 This plot does not use a theoretical distribution to compare against, but
2677 compares percentiles of two datasets. This plot does
2678 not use all raw data points, but plots the selected percentiles (quantiles) of
2679 each variable instead for the two datasets, thereby normalising the data for a
2680 direct comparison between the selected percentiles of the two dataset distributions.
2682 Quantile-quantile plots are valuable for comparing against
2683 observations and other models. Identical percentiles between the variables
2684 will lie on the one-to-one line implying the values correspond well to each
2685 other. Where there is a deviation from the one-to-one line a range of
2686 possibilities exist depending on how and where the data is shifted (e.g.,
2687 Wilks 2011 [Wilks2011]_).
2689 For distributions above the one-to-one line the distribution is left-skewed;
2690 below is right-skewed. A distinct break implies a bimodal distribution, and
2691 closer values/values further apart at the tails imply poor representation of
2692 the extremes.
2694 """
2695 # Check cubes using same functionality as the difference operator.
2696 if len(cubes) != 2:
2697 raise ValueError("cubes should contain exactly 2 cubes.")
2698 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2699 other: Cube = cubes.extract_cube(
2700 iris.Constraint(
2701 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2702 )
2703 )
2705 # Get spatial coord names.
2706 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2707 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2709 # Ensure cubes to compare are on common differencing grid.
2710 # This is triggered if either
2711 # i) latitude and longitude shapes are not the same. Note grid points
2712 # are not compared directly as these can differ through rounding
2713 # errors.
2714 # ii) or variables are known to often sit on different grid staggering
2715 # in different models (e.g. cell center vs cell edge), as is the case
2716 # for UM and LFRic comparisons.
2717 # In future greater choice of regridding method might be applied depending
2718 # on variable type. Linear regridding can in general be appropriate for smooth
2719 # variables. Care should be taken with interpretation of differences
2720 # given this dependency on regridding.
2721 if (
2722 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2723 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2724 ) or (
2725 base.long_name
2726 in [
2727 "eastward_wind_at_10m",
2728 "northward_wind_at_10m",
2729 "northward_wind_at_cell_centres",
2730 "eastward_wind_at_cell_centres",
2731 "zonal_wind_at_pressure_levels",
2732 "meridional_wind_at_pressure_levels",
2733 "potential_vorticity_at_pressure_levels",
2734 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2735 ]
2736 ):
2737 logger.debug("Linear regridding base cube to other grid to compute differences")
2738 base = regrid_onto_cube(base, other, method="Linear")
2740 # Extract just common time points.
2741 base, other = _extract_common_time_points(base, other)
2743 # Equalise attributes so we can merge.
2744 fully_equalise_attributes([base, other])
2745 logger.debug("Base: %s\nOther: %s", base, other)
2747 # Collapse cubes.
2748 base = collapse(
2749 base,
2750 coordinate=coordinates,
2751 method="PERCENTILE",
2752 additional_percent=percentiles,
2753 )
2754 other = collapse(
2755 other,
2756 coordinate=coordinates,
2757 method="PERCENTILE",
2758 additional_percent=percentiles,
2759 )
2761 # Ensure we have a name for the plot file.
2762 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2763 title = f"{recipe_title}"
2765 if filename is None:
2766 filename = slugify(recipe_title)
2768 # Add file extension.
2769 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2771 # Do the actual plotting on a scatter plot
2772 _plot_and_save_scatter_plot(
2773 base, other, plot_filename, title, one_to_one, model_names
2774 )
2776 # Add list of plots to plot metadata.
2777 plot_index = _append_to_plot_index([plot_filename])
2779 # Make a page to display the plots.
2780 _make_plot_html_page(plot_index)
2782 return iris.cube.CubeList([base, other])
2785def hinton_plot(
2786 cubes: iris.cube.CubeList, base_name: str, other_name: str, magnitude: bool = False
2787) -> None:
2788 """
2789 Plot a Hinton style triangle/scorecard plot.
2791 This plot type can be useful for summarising high level information, such as comparing
2792 how 'skillful' two models are when verified against observations for a variety of metrics,
2793 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2794 than customisable by the user as input arguments; many have been designed to automatically
2795 scale the plot depending on the number of x and y components.
2797 Parameters
2798 ----------
2799 cubes: iris.cube.CubeList
2800 A iris cubelist, containing at least two cubes of a skill metric to plot (model vs obs). Can
2801 include multiple variables, in which the plot will automatically scale for. If cubes containing
2802 a name significance_<var> exist, containing a bool array, then it will also plot whether each
2803 triangle is significant by using a thick black outline. Each cube should be 1D, with
2804 forecast_period as the only dimension.
2805 base_name: str
2806 The name of the base model to use as the control in the Hinton plot, as a string.
2807 other_name: str
2808 The name of the other model to use as the test in the Hinton plot, as a string.
2809 magnitude: bool
2810 Option bool when if True, then plot the numerical value difference in two values under each
2811 triangle.
2812 """
2813 # Ensure we have a name for the plot file.
2814 recipe_title = get_recipe_metadata().get("title", "Hinton")
2815 title = f"{recipe_title}"
2816 filename = slugify(recipe_title)
2818 # Check that all cubes only have one dimension called forecast_period
2819 for cube in cubes:
2820 if len(cube.dim_coords) > 1:
2821 raise ValueError(f"Should only have one dimension coord, {cube}")
2822 if cube.dim_coords[0].name() != "forecast_period":
2823 raise ValueError(
2824 f"Single coord should be forecast_period, not {cube.dim_coords[0].name()}"
2825 )
2827 # Separate out base cubes and other cubes.
2828 base_cubes = iris.cube.CubeList()
2829 other_cubes = iris.cube.CubeList()
2830 for c in cubes:
2831 if c.attributes["model_name"] == base_name:
2832 base_cubes.append(c)
2833 elif c.attributes["model_name"] == other_name:
2834 other_cubes.append(c)
2836 # base cubes should be the same length as other cubes, otherwise one is missing a variable.
2837 if len(base_cubes) != len(other_cubes):
2838 raise ValueError(
2839 f"base cubes {base_cubes} are not same number as {other_cubes}"
2840 )
2842 # Find common variable names in the two groups.
2843 base_vars = {cube.long_name for cube in base_cubes if cube.long_name is not None}
2844 other_vars = {cube.long_name for cube in other_cubes if cube.long_name is not None}
2845 common_vars = sorted(base_vars & other_vars)
2847 # Iterate over each variable (row)
2848 rows = []
2849 for var in common_vars:
2850 # Extract cube with matching variable name
2851 base_cube = next(
2852 (c for c in base_cubes if c.long_name == var),
2853 None,
2854 )
2856 other_cube = next(
2857 (c for c in other_cubes if c.long_name == var),
2858 None,
2859 )
2861 # If we can't find a variable in both cubes, then skip
2862 if base_cube is None or other_cube is None: 2862 ↛ 2863line 2862 didn't jump to line 2863 because the condition on line 2862 was never true
2863 continue
2865 # Compute difference (1D array)
2866 # We can already make assumption both on same forecast_periods as checked
2867 # prior to computing metric.
2868 diff = other_cube.data - base_cube.data
2870 # See if there is a significance cube present, if not, set as None.
2871 sig_cube = next(
2872 (cube for cube in cubes if cube.long_name == f"significance_{var}"),
2873 None,
2874 )
2876 # Append row information.
2877 rows.append(
2878 {
2879 "name": var,
2880 "forecast_periods": base_cube.coord("forecast_period").points,
2881 "change": diff,
2882 "significance": sig_cube.data.astype(bool)
2883 if sig_cube is not None
2884 else None,
2885 }
2886 )
2888 # For each row, compute standardised anomalies
2889 for row in rows:
2890 change = np.asarray(row["change"])
2892 anoms = change - np.mean(change)
2894 scale = np.max(np.abs(anoms))
2896 if scale > 0:
2897 scaled = anoms / scale
2898 else:
2899 scaled = np.zeros_like(anoms)
2901 row["anoms"] = anoms
2902 row["scaled"] = scaled
2904 # Setup colors of triangles
2905 color_pos = "#7CAE00"
2906 color_neg = "#7B68EE"
2908 # Setup cell/text size ratios
2909 figsize = None
2910 cell_size_in = 1.5
2911 text_row_ratio = 0.25
2913 # Get the number of x and y elements
2914 ny = len(rows)
2915 nx = max(len(row["forecast_periods"]) for row in rows)
2917 # Build non-uniform y coordinates
2918 tri_height = 1.0
2919 txt_height = text_row_ratio
2921 tri_y = []
2922 txt_y = []
2923 y_edges = [0.0]
2925 y = 0.0
2926 for _j in range(ny):
2927 tri_y.append(y + tri_height / 2)
2928 y += tri_height
2929 y_edges.append(y)
2931 if magnitude:
2932 txt_y.append(y + txt_height / 2)
2933 y += txt_height
2934 y_edges.append(y)
2936 total_height = y
2938 # Dynamic figure size
2939 if figsize is None: 2939 ↛ 2944line 2939 didn't jump to line 2944 because the condition on line 2939 was always true
2940 width = nx * cell_size_in
2941 height = total_height * cell_size_in + 2
2942 figsize = (width, height)
2944 fig, ax = plt.subplots(figsize=figsize)
2946 # Setup axes and grid.
2947 ax.set_aspect("equal", adjustable="box")
2948 ax.set_xlim(-0.5, nx - 0.5)
2949 ax.set_ylim(0, total_height)
2951 longest_row = max(rows, key=lambda row: len(row["forecast_periods"]))
2953 ax.set_xticks(np.arange(nx))
2954 ax.set_xticklabels(
2955 longest_row["forecast_periods"],
2956 rotation=90,
2957 )
2959 ax.set_yticks(tri_y)
2960 ax.set_yticklabels([row["name"] for row in rows])
2962 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2963 ax.set_yticks(y_edges, minor=True)
2965 ax.set_axisbelow(True)
2966 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2967 ax.grid(False, which="major")
2968 ax.tick_params(which="minor", length=0)
2970 ax.invert_yaxis()
2972 # Compute marker scaling (fixed overlap)
2973 fig.canvas.draw()
2975 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2976 width_in, height_in = bbox.width, bbox.height
2978 cell_w = (width_in * fig.dpi) / nx
2979 cell_h = (height_in * fig.dpi) / total_height
2980 cell_pixels = min(cell_w, cell_h)
2982 max_marker_size = (0.6 * cell_pixels) ** 2
2984 text_fontsize = cell_pixels * 0.15
2986 # Plot triangles + text
2987 for j, row in enumerate(rows):
2988 scaled = row["scaled"]
2989 anoms = row["anoms"]
2990 signif = row["significance"]
2992 for i in range(len(scaled)):
2993 val = scaled[i]
2995 if np.isnan(val): 2995 ↛ 2996line 2995 didn't jump to line 2996 because the condition on line 2995 was never true
2996 continue
2998 if abs(val) < 0.01:
2999 continue
3001 size = max_marker_size * abs(val)
3003 if val >= 0:
3004 marker = "^"
3005 color = color_pos
3006 else:
3007 marker = "v"
3008 color = color_neg
3010 if signif is not None:
3011 sig = bool(signif[i])
3012 edgecolor = "black" if sig else "none"
3013 linewidth = 0.6 if sig else 0.0
3014 else:
3015 edgecolor = "none"
3016 linewidth = 0.0
3018 ax.scatter(
3019 i,
3020 tri_y[j],
3021 s=size,
3022 marker=marker,
3023 c=color,
3024 edgecolors=edgecolor,
3025 linewidths=linewidth,
3026 zorder=3,
3027 clip_on=True,
3028 )
3030 # Text row
3031 if magnitude:
3032 mag_val = anoms[i]
3034 if not np.isnan(mag_val): 3034 ↛ 2992line 3034 didn't jump to line 2992 because the condition on line 3034 was always true
3035 ax.text(
3036 i,
3037 txt_y[j],
3038 f"{mag_val:.1f}",
3039 ha="center",
3040 va="center",
3041 fontsize=text_fontsize,
3042 color="black",
3043 zorder=4,
3044 )
3046 ax.set_title(title)
3047 plt.tight_layout()
3049 # Save plot.
3050 _save_close_figure(fig, "hinton", filename)
3052 # Add file extension.
3053 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3055 # Add list of plots to plot metadata.
3056 plot_index = _append_to_plot_index([plot_filename])
3058 # Make a page to display the plots.
3059 _make_plot_html_page(plot_index)
3062def scatter_plot(
3063 cube_x: iris.cube.Cube | iris.cube.CubeList,
3064 cube_y: iris.cube.Cube | iris.cube.CubeList,
3065 filename: str | None = None,
3066 one_to_one: bool = True,
3067 **kwargs,
3068) -> iris.cube.CubeList:
3069 """Plot a scatter plot between two variables.
3071 Both cubes must be 1D.
3073 Parameters
3074 ----------
3075 cube_x: Cube | CubeList
3076 1 dimensional Cube of the data to plot on y-axis.
3077 cube_y: Cube | CubeList
3078 1 dimensional Cube of the data to plot on x-axis.
3079 filename: str, optional
3080 Filename of the plot to write.
3081 one_to_one: bool, optional
3082 If True a 1:1 line is plotted; if False it is not. Default is True.
3084 Returns
3085 -------
3086 cubes: CubeList
3087 CubeList of the original x and y cubes for further processing.
3089 Raises
3090 ------
3091 ValueError
3092 If the cube doesn't have the right dimensions and cubes not the same
3093 size.
3094 TypeError
3095 If the cube isn't a single cube.
3097 Notes
3098 -----
3099 Scatter plots are used for determining if there is a relationship between
3100 two variables. Positive relations have a slope going from bottom left to top
3101 right; Negative relations have a slope going from top left to bottom right.
3102 """
3103 # Iterate over all cubes in cube or CubeList and plot.
3104 for cube_iter in iter_maybe(cube_x):
3105 # Check cubes are correct shape.
3106 cube_iter = check_single_cube(cube_iter)
3107 if cube_iter.ndim > 1:
3108 raise ValueError("cube_x must be 1D.")
3110 # Iterate over all cubes in cube or CubeList and plot.
3111 for cube_iter in iter_maybe(cube_y):
3112 # Check cubes are correct shape.
3113 cube_iter = check_single_cube(cube_iter)
3114 if cube_iter.ndim > 1:
3115 raise ValueError("cube_y must be 1D.")
3117 # Ensure we have a name for the plot file.
3118 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
3119 title = f"{recipe_title}"
3121 if filename is None:
3122 filename = slugify(recipe_title)
3124 # Add file extension.
3125 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3127 # Do the actual plotting.
3128 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3130 # Add list of plots to plot metadata.
3131 plot_index = _append_to_plot_index([plot_filename])
3133 # Make a page to display the plots.
3134 _make_plot_html_page(plot_index)
3136 return iris.cube.CubeList([cube_x, cube_y])
3139def vector_plot(
3140 cube_u: iris.cube.Cube,
3141 cube_v: iris.cube.Cube,
3142 filename: str | None = None,
3143 sequence_coordinate: str = "time",
3144 **kwargs,
3145) -> iris.cube.CubeList:
3146 """Plot a vector plot based on the input u and v components."""
3147 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
3149 # Cubes must have a matching sequence coordinate.
3150 try:
3151 # Check that the u and v cubes have the same sequence coordinate.
3152 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3152 ↛ anywhereline 3152 didn't jump anywhere: it always raised an exception.
3153 raise ValueError("Coordinates do not match.")
3154 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3155 raise ValueError(
3156 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3157 ) from err
3159 # Create a plot for each value of the sequence coordinate.
3160 plot_index = []
3161 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3162 for cube_u_slice, cube_v_slice in zip(
3163 cube_u.slices_over(sequence_coordinate),
3164 cube_v.slices_over(sequence_coordinate),
3165 strict=True,
3166 ):
3167 # Format the coordinate value in a unit appropriate way.
3168 seq_coord = cube_u_slice.coord(sequence_coordinate)
3169 plot_title, plot_filename = _set_title_and_filename(
3170 seq_coord, nplot, recipe_title, filename
3171 )
3173 # Do the actual plotting.
3174 _plot_and_save_vector_plot(
3175 cube_u_slice,
3176 cube_v_slice,
3177 filename=plot_filename,
3178 title=plot_title,
3179 method="pcolormesh",
3180 )
3181 plot_index.append(plot_filename)
3183 # Add list of plots to plot metadata.
3184 complete_plot_index = _append_to_plot_index(plot_index)
3186 # Make a page to display the plots.
3187 _make_plot_html_page(complete_plot_index)
3189 return iris.cube.CubeList([cube_u, cube_v])
3192def plot_histogram_series(
3193 cubes: iris.cube.Cube | iris.cube.CubeList,
3194 filename: str | None = None,
3195 sequence_coordinate: str = "time",
3196 stamp_coordinate: str = "realization",
3197 single_plot: bool = False,
3198 **kwargs,
3199) -> iris.cube.Cube | iris.cube.CubeList:
3200 """Plot a histogram plot for each vertical level provided.
3202 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3203 is present then a sequence of plots will be produced using the time slider
3204 functionality to scroll through histograms against time. If a
3205 stamp_coordinate is present then postage stamp plots will be produced. If
3206 stamp_coordinate and single_plot is True, all postage stamp plots will be
3207 plotted in a single plot instead of separate postage stamp plots.
3209 Parameters
3210 ----------
3211 cubes: Cube | iris.cube.CubeList
3212 Iris cube or CubeList of the data to plot. It should have a single dimension other
3213 than the stamp coordinate.
3214 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3215 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3216 filename: str, optional
3217 Name of the plot to write, used as a prefix for plot sequences. Defaults
3218 to the recipe name.
3219 sequence_coordinate: str, optional
3220 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3221 This coordinate must exist in the cube and will be used for the time
3222 slider.
3223 stamp_coordinate: str, optional
3224 Coordinate about which to plot postage stamp plots. Defaults to
3225 ``"realization"``.
3226 single_plot: bool, optional
3227 If True, all postage stamp plots will be plotted in a single plot. If
3228 False, each postage stamp plot will be plotted separately. Is only valid
3229 if stamp_coordinate exists and has more than a single point.
3231 Returns
3232 -------
3233 iris.cube.Cube | iris.cube.CubeList
3234 The original Cube or CubeList (so further operations can be applied).
3235 Plotted data.
3237 Raises
3238 ------
3239 ValueError
3240 If the cube doesn't have the right dimensions.
3241 TypeError
3242 If the cube isn't a Cube or CubeList.
3243 """
3244 recipe_title = get_recipe_metadata().get("title", "Histogram")
3246 cubes = iter_maybe(cubes)
3248 # Internal plotting function.
3249 plotting_func = _plot_and_save_histogram_series
3251 num_models = get_num_models(cubes)
3253 validate_cube_shape(cubes, num_models)
3255 # If several histograms are plotted, check sequence_coordinate
3256 check_sequence_coordinate(cubes, sequence_coordinate)
3258 # Get axis minimum and maximum from levels information.
3259 # If no levels set, derive minima and maxima from data in CubeList.
3260 vmin, vmax = _set_axis_range(cubes)
3262 # Make postage stamp plots if stamp_coordinate exists and has more than a
3263 # single point. If single_plot is True:
3264 # -- all postage stamp plots will be plotted in a single plot instead of
3265 # separate postage stamp plots.
3266 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3267 # produced per single model only
3268 if num_models == 1:
3269 if ( 3269 ↛ 3273line 3269 didn't jump to line 3273 because the condition on line 3269 was never true
3270 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3271 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3272 ):
3273 if single_plot:
3274 plotting_func = (
3275 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3276 )
3277 else:
3278 plotting_func = _plot_and_save_postage_stamp_histogram_series
3279 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3280 else:
3281 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3283 plot_index = []
3284 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3285 # Create a plot for each value of the sequence coordinate. Allowing for
3286 # multiple cubes in a CubeList to be plotted in the same plot for similar
3287 # sequence values. Passing a CubeList into the internal plotting function
3288 # for similar values of the sequence coordinate. cube_slice can be an
3289 # iris.cube.Cube or an iris.cube.CubeList.
3290 for cube_slice in cube_iterables:
3291 single_cube = cube_slice
3292 if isinstance(cube_slice, iris.cube.CubeList):
3293 single_cube = cube_slice[0]
3295 # Ensure valid stamp coordinate in cube dimensions
3296 if stamp_coordinate == "realization": 3296 ↛ 3299line 3296 didn't jump to line 3299 because the condition on line 3296 was always true
3297 stamp_coordinate = check_stamp_coordinate(single_cube)
3298 # Set plot titles and filename, based on sequence coordinate
3299 seq_coord = single_cube.coord(sequence_coordinate)
3300 # Use time coordinate in title and filename if single histogram output.
3301 if sequence_coordinate == "realization" and nplot == 1: 3301 ↛ 3302line 3301 didn't jump to line 3302 because the condition on line 3301 was never true
3302 seq_coord = single_cube.coord("time")
3303 # Use station name in title and filename if model vs obs comparison
3304 if sequence_coordinate == "station": 3304 ↛ 3305line 3304 didn't jump to line 3305 because the condition on line 3304 was never true
3305 seq_coord = single_cube.coord("Station_Name")
3307 plot_title, plot_filename = _set_title_and_filename(
3308 seq_coord, nplot, recipe_title, filename
3309 )
3311 # Do the actual plotting.
3312 plotting_func(
3313 cube_slice,
3314 filename=plot_filename,
3315 stamp_coordinate=stamp_coordinate,
3316 title=plot_title,
3317 vmin=vmin,
3318 vmax=vmax,
3319 )
3320 plot_index.append(plot_filename)
3322 # Add list of plots to plot metadata.
3323 complete_plot_index = _append_to_plot_index(plot_index)
3325 # Make a page to display the plots.
3326 _make_plot_html_page(complete_plot_index)
3328 return cubes
3331def plot_scatter_series(
3332 cubes: iris.cube.Cube | iris.cube.CubeList,
3333 filename: str | None = None,
3334 sequence_coordinate: str = "time",
3335 stamp_coordinate: str = "realization",
3336 hexbin: bool = False,
3337 **kwargs,
3338) -> iris.cube.Cube | iris.cube.CubeList:
3339 """Plot a scatter plot for each sequence coordinate provided.
3341 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3342 is present then a sequence of plots will be produced using the time slider
3343 functionality to scroll through scatter against time. If a
3344 stamp_coordinate is present then postage stamp plots will be produced. If
3345 stamp_coordinate and single_plot is True, all postage stamp plots will be
3346 plotted in a single plot instead of separate postage stamp plots.
3348 Parameters
3349 ----------
3350 cubes: Cube | iris.cube.CubeList
3351 Iris cube or CubeList of the data to plot. It should have a single dimension other
3352 than the stamp coordinate.
3353 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3354 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3355 filename: str, optional
3356 Name of the plot to write, used as a prefix for plot sequences. Defaults
3357 to the recipe name.
3358 sequence_coordinate: str, optional
3359 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3360 This coordinate must exist in the cube and will be used for the time
3361 slider.
3362 stamp_coordinate: str, optional
3363 Coordinate about which to plot postage stamp plots. Defaults to
3364 ``"realization"``.
3365 hexbin: bool, optional
3366 If True, generate hexbin comparison plot.
3367 If False, generate point-by-point scatter plot.
3369 Returns
3370 -------
3371 iris.cube.Cube | iris.cube.CubeList
3372 The original Cube or CubeList (so further operations can be applied).
3373 Plotted data.
3375 Raises
3376 ------
3377 ValueError
3378 If the cube doesn't have the right dimensions.
3379 TypeError
3380 If the cube isn't a Cube or CubeList.
3381 """
3382 recipe_title = get_recipe_metadata().get("title", "Scatter")
3384 cubes = iter_maybe(cubes)
3386 # Internal plotting function.
3387 plotting_func = _plot_and_save_scatter_series
3389 num_models = get_num_models(cubes)
3391 validate_cube_shape(cubes, num_models)
3393 check_sequence_coordinate(cubes, sequence_coordinate)
3395 vmin, vmax = _set_axis_range(cubes)
3397 # Require >1 models to compare on scatter plot
3398 if num_models > 1:
3399 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3400 else:
3401 raise ValueError(
3402 "Scatter plot series requires multiple number of models in input data."
3403 )
3405 plot_index = []
3406 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3407 # Create a plot for each value of the sequence coordinate. Allowing for
3408 # multiple cubes in a CubeList to be plotted in the same plot for similar
3409 # sequence values. Passing a CubeList into the internal plotting function
3410 # for similar values of the sequence coordinate. cube_slice can be an
3411 # iris.cube.Cube or an iris.cube.CubeList.
3412 for cube_slice in cube_iterables:
3413 single_cube = cube_slice
3414 if isinstance(cube_slice, iris.cube.CubeList): 3414 ↛ 3418line 3414 didn't jump to line 3418 because the condition on line 3414 was always true
3415 single_cube = cube_slice[0]
3417 # Ensure valid stamp coordinate in cube dimensions
3418 if stamp_coordinate == "realization": 3418 ↛ 3421line 3418 didn't jump to line 3421 because the condition on line 3418 was always true
3419 stamp_coordinate = check_stamp_coordinate(single_cube)
3420 # Set plot titles and filename, based on sequence coordinate
3421 seq_coord = single_cube.coord(sequence_coordinate)
3422 # Use time coordinate in title and filename if single histogram output.
3423 if sequence_coordinate == "realization" and nplot == 1:
3424 seq_coord = single_cube.coord("time")
3425 # Use station name in title and filename if model vs obs comparison
3426 if sequence_coordinate == "station":
3427 seq_coord = single_cube.coord("Station_Name")
3429 plot_title, plot_filename = _set_title_and_filename(
3430 seq_coord, nplot, recipe_title, filename
3431 )
3433 # Do the actual plotting.
3434 plotting_func(
3435 cube_slice,
3436 filename=plot_filename,
3437 stamp_coordinate=stamp_coordinate,
3438 title=plot_title,
3439 vmin=vmin,
3440 vmax=vmax,
3441 hexbin=hexbin,
3442 )
3443 plot_index.append(plot_filename)
3445 # Add list of plots to plot metadata.
3446 complete_plot_index = _append_to_plot_index(plot_index)
3448 # Make a page to display the plots.
3449 _make_plot_html_page(complete_plot_index)
3451 return cubes
3454def _plot_and_save_postage_stamp_power_spectrum_series(
3455 cubes: iris.cube.Cube,
3456 coords: list[iris.coords.Coord],
3457 stamp_coordinate: str,
3458 filename: str,
3459 title: str,
3460 series_coordinate: str | None = None,
3461 **kwargs,
3462):
3463 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3465 Parameters
3466 ----------
3467 cubes: Cube or CubeList
3468 Cube or Cubelist of the power spectrum data.
3469 coords: list[Coord]
3470 Coordinates to plot on the x-axis, one per cube.
3471 stamp_coordinate: str
3472 Coordinate that becomes different plots.
3473 filename: str
3474 Filename of the plot to write.
3475 title: str
3476 Plot title.
3477 series_coordinate: str, optional
3478 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3480 """
3481 # Use the smallest square grid that will fit the members.
3482 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3484 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3485 model_colors_map = get_model_colors_map(cubes)
3486 # ax = plt.gca()
3487 # Make a subplot for each member.
3488 for member, subplot in zip(
3489 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3490 ):
3491 ax = plt.subplot(grid_size, grid_size, subplot)
3493 # Store min/max ranges.
3494 y_levels = []
3496 line_marker = None
3497 line_width = 1
3499 for cube in iter_maybe(member):
3500 xcoord = _select_series_coord(cube, series_coordinate)
3501 xname = xcoord.points
3503 yfield = cube.data # power spectrum
3504 label = None
3505 color = "black"
3506 if model_colors_map: 3506 ↛ 3507line 3506 didn't jump to line 3507 because the condition on line 3506 was never true
3507 label = cube.attributes.get("model_name")
3508 color = model_colors_map.get(label)
3510 if member.coord(stamp_coordinate).points == [0]:
3511 ax.plot(
3512 xname,
3513 yfield,
3514 color=color,
3515 marker=line_marker,
3516 ls="-",
3517 lw=line_width,
3518 label=f"{label} (control)"
3519 if len(cube.coord(stamp_coordinate).points) > 1
3520 else label,
3521 )
3522 # Label with member if part of an ensemble and not the control.
3523 else:
3524 ax.plot(
3525 xname,
3526 yfield,
3527 color=color,
3528 ls="-",
3529 lw=1.5,
3530 alpha=0.75,
3531 label=f"{label} (member)",
3532 )
3534 # Calculate the global min/max if multiple cubes are given.
3535 _, levels, _ = colorbar_map_levels(cube, axis="y")
3536 if levels is not None: 3536 ↛ 3537line 3536 didn't jump to line 3537 because the condition on line 3536 was never true
3537 y_levels.append(min(levels))
3538 y_levels.append(max(levels))
3540 # Add some labels and tweak the style.
3541 title = f"{title}"
3542 ax.set_title(title, fontsize=16)
3544 # Set appropriate x-axis label based on coordinate
3545 if series_coordinate == "wavelength" or ( 3545 ↛ 3548line 3545 didn't jump to line 3548 because the condition on line 3545 was never true
3546 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3547 ):
3548 ax.set_xlabel("Wavelength (km)", fontsize=14)
3549 elif series_coordinate == "physical_wavenumber" or ( 3549 ↛ 3554line 3549 didn't jump to line 3554 because the condition on line 3549 was always true
3550 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3551 ):
3552 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3553 else: # frequency or check units
3554 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3555 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3556 else:
3557 ax.set_xlabel("Wavenumber", fontsize=14)
3559 ax.set_ylabel("Power Spectral Density", fontsize=14)
3560 ax.tick_params(axis="both", labelsize=12)
3562 # Set log-log scale
3563 ax.set_xscale("log")
3564 ax.set_yscale("log")
3566 # Add gridlines
3567 ax.grid(linestyle="--", color="grey", linewidth=1)
3568 # Ientify unique labels for legend
3569 handles = list(
3570 {
3571 label: handle
3572 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3573 }.values()
3574 )
3575 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3577 ax = plt.gca()
3578 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3580 # Save plot.
3581 _save_close_figure(fig, "histogram postage stamp", filename)
3584def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3585 cubes: iris.cube.Cube,
3586 coords: list[iris.coords.Coord],
3587 stamp_coordinate: str,
3588 filename: str,
3589 title: str,
3590 series_coordinate: str | None = None,
3591 **kwargs,
3592):
3593 """Plot and save power spectra for ensemble members in single plot.
3595 Parameters
3596 ----------
3597 cubes: Cube or CubeList
3598 Cube or Cubelist of the power spectrum data.
3599 coords: list[Coord]
3600 Coordinates to plot on the x-axis, one per cube.
3601 stamp_coordinate: str
3602 Coordinate that becomes different plots.
3603 filename: str
3604 Filename of the plot to write.
3605 title: str
3606 Plot title.
3607 series_coordinate: str, optional
3608 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3610 """
3611 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3612 model_colors_map = get_model_colors_map(cubes)
3614 line_marker = None
3615 line_width = 1
3617 # Compute ensemble statistics to show spread
3618 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3619 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3620 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3622 xcoord_global = mean_cube.coord(series_coordinate)
3623 x_global = xcoord_global.points
3625 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3626 xcoord = _select_series_coord(member, series_coordinate)
3627 xname = xcoord.points
3629 yfield = member.data # power spectrum
3630 color = "black"
3631 if model_colors_map: 3631 ↛ 3635line 3631 didn't jump to line 3635 because the condition on line 3631 was always true
3632 label = member.attributes.get("model_name") if i == 0 else None
3633 color = model_colors_map.get(label)
3635 if member.coord(stamp_coordinate).points == [0]:
3636 ax.plot(
3637 xname,
3638 yfield,
3639 color=color,
3640 marker=line_marker,
3641 ls="-",
3642 lw=line_width,
3643 label=f"{label} (control)"
3644 if len(member.coord(stamp_coordinate).points) > 1
3645 else label,
3646 )
3647 # Label with member number if part of an ensemble and not the control.
3648 else:
3649 ax.plot(
3650 xname,
3651 yfield,
3652 color=color,
3653 ls="-",
3654 lw=1.5,
3655 alpha=0.75,
3656 label=label,
3657 )
3659 # Set appropriate x-axis label based on coordinate
3660 if series_coordinate == "wavelength" or ( 3660 ↛ 3663line 3660 didn't jump to line 3663 because the condition on line 3660 was never true
3661 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3662 ):
3663 ax.set_xlabel("Wavelength (km)", fontsize=14)
3664 elif series_coordinate == "physical_wavenumber" or ( 3664 ↛ 3669line 3664 didn't jump to line 3669 because the condition on line 3664 was always true
3665 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3666 ):
3667 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3668 else: # frequency or check units
3669 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3670 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3671 else:
3672 ax.set_xlabel("Wavenumber", fontsize=14)
3674 # Add ensemble spread shading
3675 ax.fill_between(
3676 x_global,
3677 min_cube.data,
3678 max_cube.data,
3679 color="grey",
3680 alpha=0.3,
3681 label="Ensemble spread",
3682 )
3684 # Add ensemble mean line
3685 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3687 ax.set_ylabel("Power Spectral Density", fontsize=14)
3688 ax.tick_params(axis="both", labelsize=12)
3690 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3691 # Set log-log scale
3692 ax.set_xscale("log")
3693 ax.set_yscale("log")
3695 # Add gridlines
3696 ax.grid(linestyle="--", color="grey", linewidth=1)
3697 # Identify unique labels for legend
3698 handles = list(
3699 {
3700 label: handle
3701 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3702 }.values()
3703 )
3704 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3706 # Figure title.
3707 ax.set_title(title, fontsize=16)
3709 # Save plot.
3710 _save_close_figure(fig, "power spectra postage stamp", filename)