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