Coverage for src/CSET/operators/plot.py: 78%
927 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-30 10:39 +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
24from typing import Literal
26import cartopy.crs as ccrs
27import cartopy.feature as cfeature
28import iris
29import iris.coords
30import iris.cube
31import iris.exceptions
32import iris.plot as iplt
33import matplotlib as mpl
34import matplotlib.pyplot as plt
35import numpy as np
36from cartopy.mpl.geoaxes import GeoAxes
37from iris.cube import Cube
38from markdown_it import MarkdownIt
39from mpl_toolkits.axes_grid1.inset_locator import inset_axes
41from CSET._common import (
42 filename_slugify,
43 get_recipe_metadata,
44 iter_maybe,
45 render_file,
46 slugify,
47)
48from CSET.operators._colormaps import (
49 colorbar_map_levels,
50 get_model_colors_map,
51)
52from CSET.operators._utils import (
53 check_sequence_coordinate,
54 check_single_cube,
55 check_stamp_coordinate,
56 fully_equalise_attributes,
57 get_cube_yxcoordname,
58 get_num_models,
59 is_transect,
60 slice_over_maybe,
61 validate_cube_shape,
62 validate_cubes_coords,
63)
64from CSET.operators.collapse import collapse
65from CSET.operators.misc import _extract_common_time_points
66from CSET.operators.regrid import regrid_onto_cube
68# Use a non-interactive plotting backend.
69mpl.use("agg")
72############################
73# Private helper functions #
74############################
77def _append_to_plot_index(plot_index: list) -> list:
78 """Add plots into the plot index, returning the complete plot index."""
79 with open("meta.json", "r+t", encoding="UTF-8") as fp:
80 fcntl.flock(fp, fcntl.LOCK_EX)
81 fp.seek(0)
82 meta = json.load(fp)
83 complete_plot_index = meta.get("plots", [])
84 complete_plot_index = complete_plot_index + plot_index
85 meta["plots"] = complete_plot_index
86 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
87 os.getenv("DO_CASE_AGGREGATION")
88 ):
89 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
90 fp.seek(0)
91 fp.truncate()
92 json.dump(meta, fp, indent=2)
93 return complete_plot_index
96def _make_plot_html_page(plots: list):
97 """Create a HTML page to display a plot image."""
98 # Debug check that plots actually contains some strings.
99 assert isinstance(plots[0], str)
101 # Load HTML template file.
102 operator_files = importlib.resources.files()
103 template_file = operator_files.joinpath("_plot_page_template.html")
105 # Get some metadata.
106 meta = get_recipe_metadata()
107 title = meta.get("title", "Untitled")
108 description = MarkdownIt().render(meta.get("description", "*No description.*"))
110 # Prepare template variables.
111 variables = {
112 "title": title,
113 "description": description,
114 "initial_plot": plots[0],
115 "plots": plots,
116 "title_slug": slugify(title),
117 }
119 # Render template.
120 html = render_file(template_file, **variables)
122 # Save completed HTML.
123 with open("index.html", "wt", encoding="UTF-8") as fp:
124 fp.write(html)
127def _setup_spatial_map(
128 cube: iris.cube.Cube,
129 figure,
130 cmap,
131 grid_size: tuple[int, int] | None = None,
132 subplot: int | None = None,
133):
134 """Define map projections, extent and add coastlines and borderlines for spatial plots.
136 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
137 is specified, and map extent defined based on the input data.
139 Parameters
140 ----------
141 cube: Cube
142 2 dimensional (lat and lon) Cube of the data to plot.
143 figure:
144 Matplotlib Figure object holding all plot elements.
145 cmap:
146 Matplotlib colormap.
147 grid_size: (int, int), optional
148 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
149 subplot: int, optional
150 Subplot index if multiple spatial subplots in figure.
152 Returns
153 -------
154 axes:
155 Matplotlib GeoAxes definition.
156 """
157 # Identify min/max plot bounds.
158 try:
159 lat_axis, lon_axis = get_cube_yxcoordname(cube)
160 x1 = np.min(cube.coord(lon_axis).points)
161 x2 = np.max(cube.coord(lon_axis).points)
162 y1 = np.min(cube.coord(lat_axis).points)
163 y2 = np.max(cube.coord(lat_axis).points)
165 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
166 if np.abs(x2 - x1) > 180.0:
167 x1 = x1 - 180.0
168 x2 = x2 - 180.0
169 logging.debug("Adjusting plot bounds to fit global extent.")
171 # Consider map projection orientation.
172 # Adapting orientation enables plotting across international dateline.
173 # Users can adapt the default central_longitude if alternative projections views.
174 if x2 > 180.0 or x1 < -180.0:
175 central_longitude = 180.0
176 else:
177 central_longitude = 0.0
179 # Define spatial map projection.
180 coord_system = cube.coord(lat_axis).coord_system
181 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
182 # Define rotated pole map projection for rotated pole inputs.
183 projection = ccrs.RotatedPole(
184 pole_longitude=coord_system.grid_north_pole_longitude,
185 pole_latitude=coord_system.grid_north_pole_latitude,
186 central_rotated_longitude=central_longitude,
187 )
188 crs = projection
189 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 189 ↛ 191line 189 didn't jump to line 191 because the condition on line 189 was never true
190 # Define Transverse Mercator projection for TM inputs.
191 projection = ccrs.TransverseMercator(
192 central_longitude=coord_system.longitude_of_central_meridian,
193 central_latitude=coord_system.latitude_of_projection_origin,
194 false_easting=coord_system.false_easting,
195 false_northing=coord_system.false_northing,
196 scale_factor=coord_system.scale_factor_at_central_meridian,
197 )
198 crs = projection
199 else:
200 # Define regular map projection for non-rotated pole inputs.
201 # Alternatives might include e.g. for global model outputs:
202 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
203 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
204 projection = ccrs.PlateCarree(central_longitude=central_longitude)
205 crs = ccrs.PlateCarree()
207 # Define axes for plot (or subplot) with required map projection.
208 if subplot is not None:
209 axes = figure.add_subplot(
210 grid_size[0], grid_size[1], subplot, projection=projection
211 )
212 else:
213 axes = figure.add_subplot(projection=projection)
215 # Add coastlines and borderlines if cube contains x and y map coordinates.
216 # Avoid adding lines for masked data or specific fixed ancillary spatial plots.
217 if iris.util.is_masked(cube.data) or any( 217 ↛ 220line 217 didn't jump to line 220 because the condition on line 217 was never true
218 name in cube.name() for name in ["land_", "orography", "altitude"]
219 ):
220 pass
221 else:
222 if cmap.name in ["viridis", "Greys"]:
223 coastcol = "magenta"
224 else:
225 coastcol = "black"
226 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
227 axes.coastlines(resolution="10m", color=coastcol)
228 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
230 # Add gridlines.
231 gl = axes.gridlines(
232 alpha=0.3,
233 draw_labels=True,
234 dms=False,
235 x_inline=False,
236 y_inline=False,
237 )
238 gl.top_labels = False
239 gl.right_labels = False
240 if subplot:
241 gl.bottom_labels = False
242 gl.left_labels = False
243 if subplot % grid_size[1] == 1:
244 gl.left_labels = True
245 if subplot > ((grid_size[0] - 1) * grid_size[1]): 245 ↛ 250line 245 didn't jump to line 250 because the condition on line 245 was always true
246 gl.bottom_labels = True
248 # If is lat/lon spatial map, fix extent to keep plot tight.
249 # Specifying crs within set_extent helps ensure only data region is shown.
250 if isinstance(coord_system, iris.coord_systems.GeogCS):
251 axes.set_extent([x1, x2, y1, y2], crs=crs)
253 except ValueError:
254 # Skip if not both x and y map coordinates.
255 axes = figure.gca()
256 pass
258 return axes
261def _get_plot_resolution() -> int:
262 """Get resolution of rasterised plots in pixels per inch."""
263 return get_recipe_metadata().get("plot_resolution", 100)
266def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
267 """Return title and filename based on start and end points or bounds."""
268 if use_bounds and seq_coord.has_bounds():
269 vals = seq_coord.bounds.flatten()
270 else:
271 vals = seq_coord.points
272 start = seq_coord.units.title(vals[0])
273 end = seq_coord.units.title(vals[-1])
275 if start == end:
276 sequence_title = f"\n [{start}]"
277 sequence_fname = f"_{filename_slugify(start)}"
278 else:
279 sequence_title = f"\n [{start} to {end}]"
280 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
282 # Do not include time if coord set to zero.
283 if (
284 seq_coord.units == "hours since 0001-01-01 00:00:00"
285 and vals[0] == 0
286 and vals[-1] == 0
287 ):
288 sequence_title = ""
289 sequence_fname = ""
291 return sequence_title, sequence_fname
294def _set_title_and_filename(
295 seq_coord: iris.coords.Coord,
296 nplot: int,
297 recipe_title: str,
298 filename: str,
299):
300 """Set plot title and filename based on cube coordinate.
302 Parameters
303 ----------
304 sequence_coordinate: iris.coords.Coord
305 Coordinate about which to make a plot sequence.
306 nplot: int
307 Number of output plots to generate - controls title/naming.
308 recipe_title: str
309 Default plot title, potentially to update.
310 filename: str
311 Input plot filename, potentially to update.
313 Returns
314 -------
315 plot_title: str
316 Output formatted plot title string, based on plotted data.
317 plot_filename: str
318 Output formatted plot filename string.
319 """
320 ndim = seq_coord.ndim
321 npoints = np.size(seq_coord.points)
322 sequence_title = ""
323 sequence_fname = ""
325 # Case 1: Multiple dimension sequence input - list number of aggregated cases
326 # (e.g. aggregation histogram plots)
327 if ndim > 1:
328 ncase = np.shape(seq_coord)[0]
329 sequence_title = f"\n [{ncase} cases]"
330 sequence_fname = f"_{ncase}cases"
332 # Case 2: Single dimension input
333 else:
334 # Single sequence point
335 if npoints == 1:
336 if nplot > 1:
337 # Default labels for sequence inputs
338 sequence_value = seq_coord.units.title(seq_coord.points[0])
339 sequence_title = f"\n [{sequence_value}]"
340 sequence_fname = f"_{filename_slugify(sequence_value)}"
341 else:
342 # Aggregated attribute available where input collapsed over aggregation
343 try:
344 ncase = seq_coord.attributes["number_reference_times"]
345 sequence_title = f"\n [{ncase} cases]"
346 sequence_fname = f"_{ncase}cases"
347 except KeyError:
348 sequence_title, sequence_fname = _get_start_end_strings(
349 seq_coord, use_bounds=seq_coord.has_bounds()
350 )
351 # Multiple sequence (e.g. time) points
352 else:
353 sequence_title, sequence_fname = _get_start_end_strings(
354 seq_coord, use_bounds=False
355 )
357 # Set plot title and filename
358 plot_title = f"{recipe_title}{sequence_title}"
360 # Set plot filename, defaulting to user input if provided.
361 if filename is None:
362 filename = slugify(recipe_title)
363 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
364 else:
365 if nplot > 1:
366 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
367 else:
368 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
370 return plot_title, plot_filename
373def _select_series_coord(cube, series_coordinate):
374 """Determine the grid coordinates to use to calculate grid spacing."""
375 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
376 if series_coordinate in spacing_coordinates: 376 ↛ 382line 376 didn't jump to line 382 because the condition on line 376 was always true
377 # Try the requested coordinate first then the fallbacks in order.
378 fallbacks = [series_coordinate] + [
379 c for c in spacing_coordinates if c != series_coordinate
380 ]
381 else:
382 fallbacks = {series_coordinate}
384 # Try each possible coordinate.
385 for coord in fallbacks:
386 try:
387 return cube.coord(coord)
388 except iris.exceptions.CoordinateNotFoundError:
389 logging.debug("Coordinate %s not found.", coord)
391 # If we get here, none of the fallback options were found.
392 raise iris.exceptions.CoordinateNotFoundError(
393 f"No valid coordinate found for '{series_coordinate}' "
394 f"or fallback options {fallbacks}"
395 )
398def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
399 """Control postage stamp plot output titles based on stamp coordinate."""
400 if stamp_coord.name() == "realization":
401 mtitle = "Member"
402 else:
403 mtitle = stamp_coord.name().capitalize()
405 if stamp_coord.name() == "time":
406 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
407 else:
408 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
410 return mtitle
413def _set_axis_range(cubes):
414 """Get minimum and maximum from levels information."""
415 levels = None
416 for cube in cubes: 416 ↛ 432line 416 didn't jump to line 432 because the loop on line 416 didn't complete
417 # First check if user-specified "auto" range variable.
418 # This maintains the value of levels as None, so proceed.
419 _, levels, _ = colorbar_map_levels(cube, axis="y")
420 if levels is None:
421 break
422 # If levels is changed, recheck to use the vmin,vmax or
423 # levels-based ranges for histogram plots.
424 _, levels, _ = colorbar_map_levels(cube)
425 logging.debug("levels: %s", levels)
426 if levels is not None: 426 ↛ 416line 426 didn't jump to line 416 because the condition on line 426 was always true
427 vmin = min(levels)
428 vmax = max(levels)
429 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
430 break
432 if levels is None:
433 vmin = min(cb.data.min() for cb in cubes)
434 vmax = max(cb.data.max() for cb in cubes)
436 return vmin, vmax
439def _find_matched_slices(cubes, sequence_coordinate):
440 """Identify matched cubes in CubeList by sequence_coordinate values.
442 Ensures common points are compared for multiple cube inputs.
443 """
444 all_points = sorted(
445 set(
446 itertools.chain.from_iterable(
447 cb.coord(sequence_coordinate).points for cb in cubes
448 )
449 )
450 )
451 all_slices = list(
452 itertools.chain.from_iterable(
453 cb.slices_over(sequence_coordinate) for cb in cubes
454 )
455 )
456 # Matched slices (matched by seq coord point; it may happen that
457 # evaluated models do not cover the same seq coord range, hence matching
458 # necessary)
459 cube_iterables = [
460 iris.cube.CubeList(
461 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
462 )
463 for point in all_points
464 ]
466 return cube_iterables
469def _plot_and_save_spatial_plot(
470 cube: iris.cube.Cube,
471 filename: str,
472 title: str,
473 method: Literal["contourf", "pcolormesh"],
474 overlay_cube: iris.cube.Cube | None = None,
475 contour_cube: iris.cube.Cube | None = None,
476 **kwargs,
477):
478 """Plot and save a spatial plot.
480 Parameters
481 ----------
482 cube: Cube
483 2 dimensional (lat and lon) Cube of the data to plot.
484 filename: str
485 Filename of the plot to write.
486 title: str
487 Plot title.
488 method: "contourf" | "pcolormesh"
489 The plotting method to use.
490 overlay_cube: Cube, optional
491 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
492 contour_cube: Cube, optional
493 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
494 """
495 # Setup plot details, size, resolution, etc.
496 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
498 # Specify the color bar
499 cmap, levels, norm = colorbar_map_levels(cube)
501 # If overplotting, set required colorbars
502 if overlay_cube:
503 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
504 if contour_cube:
505 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
507 # Setup plot map projection, extent and coastlines and borderlines.
508 axes = _setup_spatial_map(cube, fig, cmap)
510 # Set colorscale bounds
511 try:
512 vmin = min(levels)
513 vmax = max(levels)
514 except TypeError:
515 vmin, vmax = None, None
516 # Ensure to use norm and not vmin/vmax if levels are defined.
517 if norm is not None:
518 vmin = None
519 vmax = None
520 logging.debug("Plotting using defined levels.")
522 # Plot the field.
523 if method == "contourf":
524 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
525 elif method == "pcolormesh":
526 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
527 else:
528 raise ValueError(f"Unknown plotting method: {method}")
530 # Overplot overlay field, if required
531 if overlay_cube:
532 try:
533 over_vmin = min(over_levels)
534 over_vmax = max(over_levels)
535 except TypeError:
536 over_vmin, over_vmax = None, None
537 if over_norm is not None: 537 ↛ 538line 537 didn't jump to line 538 because the condition on line 537 was never true
538 over_vmin = None
539 over_vmax = None
540 overlay = iplt.pcolormesh(
541 overlay_cube,
542 cmap=over_cmap,
543 norm=over_norm,
544 alpha=0.8,
545 vmin=over_vmin,
546 vmax=over_vmax,
547 )
548 # Overplot contour field, if required, with contour labelling.
549 if contour_cube:
550 contour = iplt.contour(
551 contour_cube,
552 colors="darkgray",
553 levels=cntr_levels,
554 norm=cntr_norm,
555 alpha=0.5,
556 linestyles="--",
557 linewidths=1,
558 )
559 plt.clabel(contour)
561 # Check to see if transect, and if so, adjust y axis.
562 if is_transect(cube):
563 if "pressure" in [coord.name() for coord in cube.coords()]:
564 axes.invert_yaxis()
565 axes.set_yscale("log")
566 axes.set_ylim(1100, 100)
567 # If both model_level_number and level_height exists, iplt can construct
568 # plot as a function of height above orography (NOT sea level).
569 elif {"model_level_number", "level_height"}.issubset( 569 ↛ 574line 569 didn't jump to line 574 because the condition on line 569 was always true
570 {coord.name() for coord in cube.coords()}
571 ):
572 axes.set_yscale("log")
574 axes.set_title(
575 f"{title}\n"
576 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
577 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
578 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
579 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
580 fontsize=16,
581 )
583 # Inset code
584 axins = inset_axes(
585 axes,
586 width="20%",
587 height="20%",
588 loc="upper right",
589 axes_class=GeoAxes,
590 axes_kwargs={"map_projection": ccrs.PlateCarree()},
591 )
593 # Slightly transparent to reduce plot blocking.
594 axins.patch.set_alpha(0.4)
596 axins.coastlines(resolution="50m")
597 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
599 SLat, SLon, ELat, ELon = (
600 float(coord) for coord in cube.attributes["transect_coords"].split("_")
601 )
603 # Draw line between them
604 axins.plot(
605 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
606 )
608 # Plot points (note: lon, lat order for Cartopy)
609 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
610 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
612 lon_min, lon_max = sorted([SLon, ELon])
613 lat_min, lat_max = sorted([SLat, ELat])
615 # Midpoints
616 lon_mid = (lon_min + lon_max) / 2
617 lat_mid = (lat_min + lat_max) / 2
619 # Maximum half-range
620 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
621 if half_range == 0: # points identical → provide small default 621 ↛ 625line 621 didn't jump to line 625 because the condition on line 621 was always true
622 half_range = 1
624 # Set square extent
625 axins.set_extent(
626 [
627 lon_mid - half_range,
628 lon_mid + half_range,
629 lat_mid - half_range,
630 lat_mid + half_range,
631 ],
632 crs=ccrs.PlateCarree(),
633 )
635 # Ensure square aspect
636 axins.set_aspect("equal")
638 else:
639 # Add title.
640 axes.set_title(title, fontsize=16)
642 # Adjust padding if spatial plot or transect
643 if is_transect(cube):
644 yinfopad = -0.1
645 ycbarpad = 0.1
646 else:
647 yinfopad = 0.01
648 ycbarpad = 0.042
650 # Add watermark with min/max/mean. Currently not user togglable.
651 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
652 axes.annotate(
653 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
654 xy=(0.025, yinfopad),
655 xycoords="axes fraction",
656 xytext=(-5, 5),
657 textcoords="offset points",
658 ha="left",
659 va="bottom",
660 size=11,
661 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
662 )
664 # Add secondary colour bar for overlay_cube field if required.
665 if overlay_cube:
666 cbarB = fig.colorbar(
667 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
668 )
669 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
670 # add ticks and tick_labels for every levels if less than 20 levels exist
671 if over_levels is not None and len(over_levels) < 20: 671 ↛ 672line 671 didn't jump to line 672 because the condition on line 671 was never true
672 cbarB.set_ticks(over_levels)
673 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
674 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
675 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
676 logging.debug("Set secondary colorbar ticks and labels.")
678 # Add main colour bar.
679 cbar = fig.colorbar(
680 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
681 )
683 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
684 # add ticks and tick_labels for every levels if less than 20 levels exist
685 if levels is not None and len(levels) < 20:
686 cbar.set_ticks(levels)
687 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
688 if "rainfall" or "snowfall" or "visibility" in cube.name(): 688 ↛ 690line 688 didn't jump to line 690 because the condition on line 688 was always true
689 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
690 logging.debug("Set colorbar ticks and labels.")
692 # Save plot.
693 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
694 logging.info("Saved spatial plot to %s", filename)
695 plt.close(fig)
698def _plot_and_save_postage_stamp_spatial_plot(
699 cube: iris.cube.Cube,
700 filename: str,
701 stamp_coordinate: str,
702 title: str,
703 method: Literal["contourf", "pcolormesh"],
704 overlay_cube: iris.cube.Cube | None = None,
705 contour_cube: iris.cube.Cube | None = None,
706 **kwargs,
707):
708 """Plot postage stamp spatial plots from an ensemble.
710 Parameters
711 ----------
712 cube: Cube
713 Iris cube of data to be plotted. It must have the stamp coordinate.
714 filename: str
715 Filename of the plot to write.
716 stamp_coordinate: str
717 Coordinate that becomes different plots.
718 method: "contourf" | "pcolormesh"
719 The plotting method to use.
720 overlay_cube: Cube, optional
721 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
722 contour_cube: Cube, optional
723 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
725 Raises
726 ------
727 ValueError
728 If the cube doesn't have the right dimensions.
729 """
730 # Use the smallest square grid that will fit the members.
731 nmember = len(cube.coord(stamp_coordinate).points)
732 grid_rows = int(math.sqrt(nmember))
733 grid_size = math.ceil(nmember / grid_rows)
735 fig = plt.figure(
736 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
737 )
739 # Specify the color bar
740 cmap, levels, norm = colorbar_map_levels(cube)
741 # If overplotting, set required colorbars
742 if overlay_cube: 742 ↛ 743line 742 didn't jump to line 743 because the condition on line 742 was never true
743 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
744 if contour_cube: 744 ↛ 745line 744 didn't jump to line 745 because the condition on line 744 was never true
745 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
747 # Make a subplot for each member.
748 for member, subplot in zip(
749 cube.slices_over(stamp_coordinate),
750 range(1, grid_size * grid_rows + 1),
751 strict=False,
752 ):
753 # Setup subplot map projection, extent and coastlines and borderlines.
754 axes = _setup_spatial_map(
755 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
756 )
757 if method == "contourf":
758 # Filled contour plot of the field.
759 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
760 elif method == "pcolormesh":
761 if levels is not None:
762 vmin = min(levels)
763 vmax = max(levels)
764 else:
765 raise TypeError("Unknown vmin and vmax range.")
766 vmin, vmax = None, None
767 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
768 # if levels are defined.
769 if norm is not None: 769 ↛ 770line 769 didn't jump to line 770 because the condition on line 769 was never true
770 vmin = None
771 vmax = None
772 # pcolormesh plot of the field.
773 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
774 else:
775 raise ValueError(f"Unknown plotting method: {method}")
777 # Overplot overlay field, if required
778 if overlay_cube: 778 ↛ 779line 778 didn't jump to line 779 because the condition on line 778 was never true
779 try:
780 over_vmin = min(over_levels)
781 over_vmax = max(over_levels)
782 except TypeError:
783 over_vmin, over_vmax = None, None
784 if over_norm is not None:
785 over_vmin = None
786 over_vmax = None
787 iplt.pcolormesh(
788 overlay_cube[member.coord(stamp_coordinate).points[0]],
789 cmap=over_cmap,
790 norm=over_norm,
791 alpha=0.6,
792 vmin=over_vmin,
793 vmax=over_vmax,
794 )
795 # Overplot contour field, if required
796 if contour_cube: 796 ↛ 797line 796 didn't jump to line 797 because the condition on line 796 was never true
797 iplt.contour(
798 contour_cube[member.coord(stamp_coordinate).points[0]],
799 colors="darkgray",
800 levels=cntr_levels,
801 norm=cntr_norm,
802 alpha=0.6,
803 linestyles="--",
804 linewidths=1,
805 )
806 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
807 axes.set_title(f"{mtitle}")
809 # Put the shared colorbar in its own axes.
810 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
811 colorbar = fig.colorbar(
812 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
813 )
814 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
816 # Overall figure title.
817 fig.suptitle(title, fontsize=16)
819 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
820 logging.info("Saved contour postage stamp plot to %s", filename)
821 plt.close(fig)
824def _plot_and_save_line_series(
825 cubes: iris.cube.CubeList,
826 coords: list[iris.coords.Coord],
827 ensemble_coord: str,
828 filename: str,
829 title: str,
830 **kwargs,
831):
832 """Plot and save a 1D line series.
834 Parameters
835 ----------
836 cubes: Cube or CubeList
837 Cube or CubeList containing the cubes to plot on the y-axis.
838 coords: list[Coord]
839 Coordinates to plot on the x-axis, one per cube.
840 ensemble_coord: str
841 Ensemble coordinate in the cube.
842 filename: str
843 Filename of the plot to write.
844 title: str
845 Plot title.
846 """
847 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
849 model_colors_map = get_model_colors_map(cubes)
851 # Store min/max ranges.
852 y_levels = []
854 # Check match-up across sequence coords gives consistent sizes
855 validate_cubes_coords(cubes, coords)
857 for cube, coord in zip(cubes, coords, strict=True):
858 label = None
859 color = "black"
860 if model_colors_map:
861 label = cube.attributes.get("model_name")
862 color = model_colors_map.get(label)
863 for cube_slice in cube.slices_over(ensemble_coord):
864 # Label with (control) if part of an ensemble or not otherwise.
865 if cube_slice.coord(ensemble_coord).points == [0]:
866 iplt.plot(
867 coord,
868 cube_slice,
869 color=color,
870 marker="o",
871 ls="-",
872 lw=3,
873 label=f"{label} (control)"
874 if len(cube.coord(ensemble_coord).points) > 1
875 else label,
876 )
877 # Label with (perturbed) if part of an ensemble and not the control.
878 else:
879 iplt.plot(
880 coord,
881 cube_slice,
882 color=color,
883 ls="-",
884 lw=1.5,
885 alpha=0.75,
886 label=f"{label} (member)",
887 )
889 # Calculate the global min/max if multiple cubes are given.
890 _, levels, _ = colorbar_map_levels(cube, axis="y")
891 if levels is not None: 891 ↛ 892line 891 didn't jump to line 892 because the condition on line 891 was never true
892 y_levels.append(min(levels))
893 y_levels.append(max(levels))
895 # Get the current axes.
896 ax = plt.gca()
898 # Add some labels and tweak the style.
899 # check if cubes[0] works for single cube if not CubeList
900 if coords[0].name() == "time":
901 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
902 else:
903 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
904 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
905 ax.set_title(title, fontsize=16)
907 ax.ticklabel_format(axis="y", useOffset=False)
908 ax.tick_params(axis="x", labelrotation=15)
909 ax.tick_params(axis="both", labelsize=12)
911 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
912 if y_levels: 912 ↛ 913line 912 didn't jump to line 913 because the condition on line 912 was never true
913 ax.set_ylim(min(y_levels), max(y_levels))
914 # Add zero line.
915 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
916 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
917 logging.debug(
918 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
919 )
920 else:
921 ax.autoscale()
923 # Add gridlines
924 ax.grid(linestyle="--", color="grey", linewidth=1)
925 # Ientify unique labels for legend
926 handles = list(
927 {
928 label: handle
929 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
930 }.values()
931 )
932 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
934 # Save plot.
935 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
936 logging.info("Saved line plot to %s", filename)
937 plt.close(fig)
940def _plot_and_save_line_power_spectrum_series(
941 cubes: iris.cube.Cube | iris.cube.CubeList,
942 coords: list[iris.coords.Coord],
943 ensemble_coord: str,
944 filename: str,
945 title: str,
946 series_coordinate: str,
947 **kwargs,
948):
949 """Plot and save a 1D line series.
951 Parameters
952 ----------
953 cubes: Cube or CubeList
954 Cube or CubeList containing the cubes to plot on the y-axis.
955 coords: list[Coord]
956 Coordinates to plot on the x-axis, one per cube.
957 ensemble_coord: str
958 Ensemble coordinate in the cube.
959 filename: str
960 Filename of the plot to write.
961 title: str
962 Plot title.
963 series_coordinate: str
964 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
965 """
966 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
967 model_colors_map = get_model_colors_map(cubes)
968 ax = plt.gca()
970 # Store min/max ranges.
971 y_levels = []
973 line_marker = None
974 line_width = 1
976 for cube in iter_maybe(cubes):
977 # next 2 lines replace chunk of code.
978 xcoord = _select_series_coord(cube, series_coordinate)
979 xname = xcoord.points
981 yfield = cube.data # power spectrum
982 label = None
983 color = "black"
984 if model_colors_map: 984 ↛ 987line 984 didn't jump to line 987 because the condition on line 984 was always true
985 label = cube.attributes.get("model_name")
986 color = model_colors_map.get(label)
987 for cube_slice in cube.slices_over(ensemble_coord):
988 # Label with (control) if part of an ensemble or not otherwise.
989 if cube_slice.coord(ensemble_coord).points == [0]: 989 ↛ 1003line 989 didn't jump to line 1003 because the condition on line 989 was always true
990 ax.plot(
991 xname,
992 yfield,
993 color=color,
994 marker=line_marker,
995 ls="-",
996 lw=line_width,
997 label=f"{label} (control)"
998 if len(cube.coord(ensemble_coord).points) > 1
999 else label,
1000 )
1001 # Label with (perturbed) if part of an ensemble and not the control.
1002 else:
1003 ax.plot(
1004 xname,
1005 yfield,
1006 color=color,
1007 ls="-",
1008 lw=1.5,
1009 alpha=0.75,
1010 label=f"{label} (member)",
1011 )
1013 # Calculate the global min/max if multiple cubes are given.
1014 _, levels, _ = colorbar_map_levels(cube, axis="y")
1015 if levels is not None: 1015 ↛ 1016line 1015 didn't jump to line 1016 because the condition on line 1015 was never true
1016 y_levels.append(min(levels))
1017 y_levels.append(max(levels))
1019 # Add some labels and tweak the style.
1021 title = f"{title}"
1022 ax.set_title(title, fontsize=16)
1024 # Set appropriate x-axis label based on coordinate
1025 if series_coordinate == "wavelength" or ( 1025 ↛ 1028line 1025 didn't jump to line 1028 because the condition on line 1025 was never true
1026 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1027 ):
1028 ax.set_xlabel("Wavelength (km)", fontsize=14)
1029 elif series_coordinate == "physical_wavenumber" or ( 1029 ↛ 1032line 1029 didn't jump to line 1032 because the condition on line 1029 was never true
1030 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1031 ):
1032 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1033 else: # frequency or check units
1034 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1034 ↛ 1035line 1034 didn't jump to line 1035 because the condition on line 1034 was never true
1035 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1036 else:
1037 ax.set_xlabel("Wavenumber", fontsize=14)
1039 ax.set_ylabel("Power Spectral Density", fontsize=14)
1040 ax.tick_params(axis="both", labelsize=12)
1042 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1044 # Set log-log scale
1045 ax.set_xscale("log")
1046 ax.set_yscale("log")
1048 # Add gridlines
1049 ax.grid(linestyle="--", color="grey", linewidth=1)
1050 # Ientify unique labels for legend
1051 handles = list(
1052 {
1053 label: handle
1054 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1055 }.values()
1056 )
1057 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1059 # Save plot.
1060 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1061 logging.info("Saved line plot to %s", filename)
1062 plt.close(fig)
1065def _plot_and_save_vertical_line_series(
1066 cubes: iris.cube.CubeList,
1067 coords: list[iris.coords.Coord],
1068 ensemble_coord: str,
1069 filename: str,
1070 series_coordinate: str,
1071 title: str,
1072 vmin: float,
1073 vmax: float,
1074 **kwargs,
1075):
1076 """Plot and save a 1D line series in vertical.
1078 Parameters
1079 ----------
1080 cubes: CubeList
1081 1 dimensional Cube or CubeList of the data to plot on x-axis.
1082 coord: list[Coord]
1083 Coordinates to plot on the y-axis, one per cube.
1084 ensemble_coord: str
1085 Ensemble coordinate in the cube.
1086 filename: str
1087 Filename of the plot to write.
1088 series_coordinate: str
1089 Coordinate to use as vertical axis.
1090 title: str
1091 Plot title.
1092 vmin: float
1093 Minimum value for the x-axis.
1094 vmax: float
1095 Maximum value for the x-axis.
1096 """
1097 # plot the vertical pressure axis using log scale
1098 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1100 model_colors_map = get_model_colors_map(cubes)
1102 # Check match-up across sequence coords gives consistent sizes
1103 validate_cubes_coords(cubes, coords)
1105 for cube, coord in zip(cubes, coords, strict=True):
1106 label = None
1107 color = "black"
1108 if model_colors_map: 1108 ↛ 1109line 1108 didn't jump to line 1109 because the condition on line 1108 was never true
1109 label = cube.attributes.get("model_name")
1110 color = model_colors_map.get(label)
1112 for cube_slice in cube.slices_over(ensemble_coord):
1113 # If ensemble data given plot control member with (control)
1114 # unless single forecast.
1115 if cube_slice.coord(ensemble_coord).points == [0]:
1116 iplt.plot(
1117 cube_slice,
1118 coord,
1119 color=color,
1120 marker="o",
1121 ls="-",
1122 lw=3,
1123 label=f"{label} (control)"
1124 if len(cube.coord(ensemble_coord).points) > 1
1125 else label,
1126 )
1127 # If ensemble data given plot perturbed members with (perturbed).
1128 else:
1129 iplt.plot(
1130 cube_slice,
1131 coord,
1132 color=color,
1133 ls="-",
1134 lw=1.5,
1135 alpha=0.75,
1136 label=f"{label} (member)",
1137 )
1139 # Get the current axis
1140 ax = plt.gca()
1142 # Special handling for pressure level data.
1143 if series_coordinate == "pressure": 1143 ↛ 1165line 1143 didn't jump to line 1165 because the condition on line 1143 was always true
1144 # Invert y-axis and set to log scale.
1145 ax.invert_yaxis()
1146 ax.set_yscale("log")
1148 # Define y-ticks and labels for pressure log axis.
1149 y_tick_labels = [
1150 "1000",
1151 "850",
1152 "700",
1153 "500",
1154 "300",
1155 "200",
1156 "100",
1157 ]
1158 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1160 # Set y-axis limits and ticks.
1161 ax.set_ylim(1100, 100)
1163 # Test if series_coordinate is model level data. The UM data uses
1164 # model_level_number and lfric uses full_levels as coordinate.
1165 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1166 # Define y-ticks and labels for vertical axis.
1167 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1168 y_tick_labels = [str(int(i)) for i in y_ticks]
1169 ax.set_ylim(min(y_ticks), max(y_ticks))
1171 ax.set_yticks(y_ticks)
1172 ax.set_yticklabels(y_tick_labels)
1174 # Set x-axis limits.
1175 ax.set_xlim(vmin, vmax)
1176 # Mark y=0 if present in plot.
1177 if vmin < 0.0 and vmax > 0.0: 1177 ↛ 1178line 1177 didn't jump to line 1178 because the condition on line 1177 was never true
1178 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1180 # Add some labels and tweak the style.
1181 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1182 ax.set_xlabel(
1183 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1184 )
1185 ax.set_title(title, fontsize=16)
1186 ax.ticklabel_format(axis="x")
1187 ax.tick_params(axis="y")
1188 ax.tick_params(axis="both", labelsize=12)
1190 # Add gridlines
1191 ax.grid(linestyle="--", color="grey", linewidth=1)
1192 # Ientify unique labels for legend
1193 handles = list(
1194 {
1195 label: handle
1196 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1197 }.values()
1198 )
1199 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1201 # Save plot.
1202 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1203 logging.info("Saved line plot to %s", filename)
1204 plt.close(fig)
1207def _plot_and_save_scatter_plot(
1208 cube_x: iris.cube.Cube | iris.cube.CubeList,
1209 cube_y: iris.cube.Cube | iris.cube.CubeList,
1210 filename: str,
1211 title: str,
1212 one_to_one: bool,
1213 model_names: list[str] = None,
1214 **kwargs,
1215):
1216 """Plot and save a 2D scatter plot.
1218 Parameters
1219 ----------
1220 cube_x: Cube | CubeList
1221 1 dimensional Cube or CubeList of the data to plot on x-axis.
1222 cube_y: Cube | CubeList
1223 1 dimensional Cube or CubeList of the data to plot on y-axis.
1224 filename: str
1225 Filename of the plot to write.
1226 title: str
1227 Plot title.
1228 one_to_one: bool
1229 Whether a 1:1 line is plotted.
1230 """
1231 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1232 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1233 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1234 # over the pairs simultaneously.
1236 # Ensure cube_x and cube_y are iterable
1237 cube_x_iterable = iter_maybe(cube_x)
1238 cube_y_iterable = iter_maybe(cube_y)
1240 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1241 iplt.scatter(cube_x_iter, cube_y_iter)
1242 if one_to_one is True:
1243 plt.plot(
1244 [
1245 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1246 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1247 ],
1248 [
1249 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1250 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1251 ],
1252 "k",
1253 linestyle="--",
1254 )
1255 ax = plt.gca()
1257 # Add some labels and tweak the style.
1258 if model_names is None:
1259 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1260 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1261 else:
1262 # Add the model names, these should be order of base (x) and other (y).
1263 ax.set_xlabel(
1264 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1265 )
1266 ax.set_ylabel(
1267 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1268 )
1269 ax.set_title(title, fontsize=16)
1270 ax.ticklabel_format(axis="y", useOffset=False)
1271 ax.tick_params(axis="x", labelrotation=15)
1272 ax.tick_params(axis="both", labelsize=12)
1273 ax.autoscale()
1275 # Save plot.
1276 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1277 logging.info("Saved scatter plot to %s", filename)
1278 plt.close(fig)
1281def _plot_and_save_vector_plot(
1282 cube_u: iris.cube.Cube,
1283 cube_v: iris.cube.Cube,
1284 filename: str,
1285 title: str,
1286 method: Literal["contourf", "pcolormesh"],
1287 **kwargs,
1288):
1289 """Plot and save a 2D vector plot.
1291 Parameters
1292 ----------
1293 cube_u: Cube
1294 2 dimensional Cube of u component of the data.
1295 cube_v: Cube
1296 2 dimensional Cube of v component of the data.
1297 filename: str
1298 Filename of the plot to write.
1299 title: str
1300 Plot title.
1301 """
1302 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1304 # Create a cube containing the magnitude of the vector field.
1305 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1306 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1308 # Specify the color bar
1309 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1311 # Setup plot map projection, extent and coastlines and borderlines.
1312 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1314 if method == "contourf":
1315 # Filled contour plot of the field.
1316 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1317 elif method == "pcolormesh":
1318 try:
1319 vmin = min(levels)
1320 vmax = max(levels)
1321 except TypeError:
1322 vmin, vmax = None, None
1323 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1324 # if levels are defined.
1325 if norm is not None:
1326 vmin = None
1327 vmax = None
1328 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1329 else:
1330 raise ValueError(f"Unknown plotting method: {method}")
1332 # Check to see if transect, and if so, adjust y axis.
1333 if is_transect(cube_vec_mag):
1334 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1335 axes.invert_yaxis()
1336 axes.set_yscale("log")
1337 axes.set_ylim(1100, 100)
1338 # If both model_level_number and level_height exists, iplt can construct
1339 # plot as a function of height above orography (NOT sea level).
1340 elif {"model_level_number", "level_height"}.issubset(
1341 {coord.name() for coord in cube_vec_mag.coords()}
1342 ):
1343 axes.set_yscale("log")
1345 axes.set_title(
1346 f"{title}\n"
1347 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1348 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1349 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1350 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1351 fontsize=16,
1352 )
1354 else:
1355 # Add title.
1356 axes.set_title(title, fontsize=16)
1358 # Add watermark with min/max/mean. Currently not user togglable.
1359 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1360 axes.annotate(
1361 f"Min: {np.min(cube_vec_mag.data):.3g} Max: {np.max(cube_vec_mag.data):.3g} Mean: {np.mean(cube_vec_mag.data):.3g}",
1362 xy=(0.05, -0.05),
1363 xycoords="axes fraction",
1364 xytext=(-5, 5),
1365 textcoords="offset points",
1366 ha="right",
1367 va="bottom",
1368 size=11,
1369 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1370 )
1372 # Add colour bar.
1373 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1374 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1375 # add ticks and tick_labels for every levels if less than 20 levels exist
1376 if levels is not None and len(levels) < 20:
1377 cbar.set_ticks(levels)
1378 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1380 # 30 barbs along the longest axis of the plot, or a barb per point for data
1381 # with less than 30 points.
1382 step = max(max(cube_u.shape) // 30, 1)
1383 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1385 # Save plot.
1386 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1387 logging.info("Saved vector plot to %s", filename)
1388 plt.close(fig)
1391def _plot_and_save_histogram_series(
1392 cubes: iris.cube.Cube | iris.cube.CubeList,
1393 filename: str,
1394 title: str,
1395 vmin: float,
1396 vmax: float,
1397 **kwargs,
1398):
1399 """Plot and save a histogram series.
1401 Parameters
1402 ----------
1403 cubes: Cube or CubeList
1404 2 dimensional Cube or CubeList of the data to plot as histogram.
1405 filename: str
1406 Filename of the plot to write.
1407 title: str
1408 Plot title.
1409 vmin: float
1410 minimum for colorbar
1411 vmax: float
1412 maximum for colorbar
1413 """
1414 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1415 ax = plt.gca()
1417 model_colors_map = get_model_colors_map(cubes)
1419 # Set default that histograms will produce probability density function
1420 # at each bin (integral over range sums to 1).
1421 density = True
1423 for cube in iter_maybe(cubes):
1424 # Easier to check title (where var name originates)
1425 # than seeing if long names exist etc.
1426 # Exception case, where distribution better fits log scales/bins.
1427 if "surface_microphysical" in title:
1428 if "amount" in title: 1428 ↛ 1430line 1428 didn't jump to line 1430 because the condition on line 1428 was never true
1429 # Compute histogram following Klingaman et al. (2017): ASoP
1430 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1431 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1432 density = False
1433 else:
1434 bins = 10.0 ** (
1435 np.arange(-10, 27, 1) / 10.0
1436 ) # Suggestion from RMED toolbox.
1437 bins = np.insert(bins, 0, 0)
1438 ax.set_yscale("log")
1439 vmin = bins[1]
1440 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1441 ax.set_xscale("log")
1442 elif "lightning" in title:
1443 bins = [0, 1, 2, 3, 4, 5]
1444 else:
1445 bins = np.linspace(vmin, vmax, 51)
1446 logging.debug(
1447 "Plotting histogram with %s bins %s - %s.",
1448 np.size(bins),
1449 np.min(bins),
1450 np.max(bins),
1451 )
1453 # Reshape cube data into a single array to allow for a single histogram.
1454 # Otherwise we plot xdim histograms stacked.
1455 cube_data_1d = (cube.data).flatten()
1457 label = None
1458 color = "black"
1459 if model_colors_map: 1459 ↛ 1460line 1459 didn't jump to line 1460 because the condition on line 1459 was never true
1460 label = cube.attributes.get("model_name")
1461 color = model_colors_map[label]
1462 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1464 # Compute area under curve.
1465 if "surface_microphysical" in title and "amount" in title: 1465 ↛ 1466line 1465 didn't jump to line 1466 because the condition on line 1465 was never true
1466 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1467 x = x * bin_mean / x.sum()
1468 x = x[1:]
1469 y = y[1:]
1471 ax.plot(
1472 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1473 )
1475 # Add some labels and tweak the style.
1476 ax.set_title(title, fontsize=16)
1477 ax.set_xlabel(
1478 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1479 )
1480 ax.set_ylabel("Normalised probability density", fontsize=14)
1481 if "surface_microphysical" in title and "amount" in title: 1481 ↛ 1482line 1481 didn't jump to line 1482 because the condition on line 1481 was never true
1482 ax.set_ylabel(
1483 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1484 )
1485 ax.set_xlim(vmin, vmax)
1486 ax.tick_params(axis="both", labelsize=12)
1488 # Overlay grid-lines onto histogram plot.
1489 ax.grid(linestyle="--", color="grey", linewidth=1)
1490 if model_colors_map: 1490 ↛ 1491line 1490 didn't jump to line 1491 because the condition on line 1490 was never true
1491 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1493 # Save plot.
1494 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1495 logging.info("Saved histogram plot to %s", filename)
1496 plt.close(fig)
1499def _plot_and_save_postage_stamp_histogram_series(
1500 cube: iris.cube.Cube,
1501 filename: str,
1502 title: str,
1503 stamp_coordinate: str,
1504 vmin: float,
1505 vmax: float,
1506 **kwargs,
1507):
1508 """Plot and save postage (ensemble members) stamps for a histogram series.
1510 Parameters
1511 ----------
1512 cube: Cube
1513 2 dimensional Cube of the data to plot as histogram.
1514 filename: str
1515 Filename of the plot to write.
1516 title: str
1517 Plot title.
1518 stamp_coordinate: str
1519 Coordinate that becomes different plots.
1520 vmin: float
1521 minimum for pdf x-axis
1522 vmax: float
1523 maximum for pdf x-axis
1524 """
1525 # Use the smallest square grid that will fit the members.
1526 nmember = len(cube.coord(stamp_coordinate).points)
1527 grid_rows = int(math.sqrt(nmember))
1528 grid_size = math.ceil(nmember / grid_rows)
1530 fig = plt.figure(
1531 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1532 )
1533 # Make a subplot for each member.
1534 for member, subplot in zip(
1535 cube.slices_over(stamp_coordinate),
1536 range(1, grid_size * grid_rows + 1),
1537 strict=False,
1538 ):
1539 # Implicit interface is much easier here, due to needing to have the
1540 # cartopy GeoAxes generated.
1541 plt.subplot(grid_rows, grid_size, subplot)
1542 # Reshape cube data into a single array to allow for a single histogram.
1543 # Otherwise we plot xdim histograms stacked.
1544 member_data_1d = (member.data).flatten()
1545 plt.hist(member_data_1d, density=True, stacked=True)
1546 axes = plt.gca()
1547 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1548 axes.set_title(f"{mtitle}")
1549 axes.set_xlim(vmin, vmax)
1551 # Overall figure title.
1552 fig.suptitle(title, fontsize=16)
1554 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1555 logging.info("Saved histogram postage stamp plot to %s", filename)
1556 plt.close(fig)
1559def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1560 cube: iris.cube.Cube,
1561 filename: str,
1562 title: str,
1563 stamp_coordinate: str,
1564 vmin: float,
1565 vmax: float,
1566 **kwargs,
1567):
1568 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1569 ax.set_title(title, fontsize=16)
1570 ax.set_xlim(vmin, vmax)
1571 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1572 ax.set_ylabel("normalised probability density", fontsize=14)
1573 # Loop over all slices along the stamp_coordinate
1574 for member in cube.slices_over(stamp_coordinate):
1575 # Flatten the member data to 1D
1576 member_data_1d = member.data.flatten()
1577 # Plot the histogram using plt.hist
1578 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1579 plt.hist(
1580 member_data_1d,
1581 density=True,
1582 stacked=True,
1583 label=f"{mtitle}",
1584 )
1586 # Add a legend
1587 ax.legend(fontsize=16)
1589 # Save the figure to a file
1590 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1591 logging.info("Saved histogram postage stamp plot to %s", filename)
1593 # Close the figure
1594 plt.close(fig)
1597def _plot_and_save_scattermap_plot(
1598 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1599):
1600 """Plot and save a geographical scatter plot.
1602 Parameters
1603 ----------
1604 cube: Cube
1605 1 dimensional Cube of the data points with auxiliary latitude and
1606 longitude coordinates,
1607 filename: str
1608 Filename of the plot to write.
1609 title: str
1610 Plot title.
1611 projection: str
1612 Mapping projection to be used by cartopy.
1613 """
1614 # Setup plot details, size, resolution, etc.
1615 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1616 if projection is not None:
1617 # Apart from the default, the only projection we currently support is
1618 # a stereographic projection over the North Pole.
1619 if projection == "NP_Stereo":
1620 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1621 else:
1622 raise ValueError(f"Unknown projection: {projection}")
1623 else:
1624 axes = plt.axes(projection=ccrs.PlateCarree())
1626 # Scatter plot of the field. The marker size is chosen to give
1627 # symbols that decrease in size as the number of observations
1628 # increases, although the fraction of the figure covered by
1629 # symbols increases roughly as N^(1/2), disregarding overlaps,
1630 # and has been selected for the default figure size of (10, 10).
1631 # Should this be changed, the marker size should be adjusted in
1632 # proportion to the area of the figure.
1633 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1634 klon = None
1635 klat = None
1636 for kc in range(len(cube.aux_coords)):
1637 if cube.aux_coords[kc].standard_name == "latitude":
1638 klat = kc
1639 elif cube.aux_coords[kc].standard_name == "longitude":
1640 klon = kc
1641 scatter_map = iplt.scatter(
1642 cube.aux_coords[klon],
1643 cube.aux_coords[klat],
1644 c=cube.data[:],
1645 s=mrk_size,
1646 cmap="jet",
1647 edgecolors="k",
1648 )
1650 # Add coastlines and borderlines.
1651 try:
1652 axes.coastlines(resolution="10m")
1653 axes.add_feature(cfeature.BORDERS)
1654 except AttributeError:
1655 pass
1657 # Add title.
1658 axes.set_title(title, fontsize=16)
1660 # Add colour bar.
1661 cbar = fig.colorbar(scatter_map)
1662 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1664 # Save plot.
1665 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1666 logging.info("Saved geographical scatter plot to %s", filename)
1667 plt.close(fig)
1670def _spatial_plot(
1671 method: Literal["contourf", "pcolormesh"],
1672 cube: iris.cube.Cube,
1673 filename: str | None,
1674 sequence_coordinate: str,
1675 stamp_coordinate: str,
1676 overlay_cube: iris.cube.Cube | None = None,
1677 contour_cube: iris.cube.Cube | None = None,
1678 **kwargs,
1679):
1680 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1682 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1683 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1684 is present then postage stamp plots will be produced.
1686 If an overlay_cube and/or contour_cube are specified, multiple variables can
1687 be overplotted on the same figure.
1689 Parameters
1690 ----------
1691 method: "contourf" | "pcolormesh"
1692 The plotting method to use.
1693 cube: Cube
1694 Iris cube of the data to plot. It should have two spatial dimensions,
1695 such as lat and lon, and may also have a another two dimension to be
1696 plotted sequentially and/or as postage stamp plots.
1697 filename: str | None
1698 Name of the plot to write, used as a prefix for plot sequences. If None
1699 uses the recipe name.
1700 sequence_coordinate: str
1701 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1702 This coordinate must exist in the cube.
1703 stamp_coordinate: str
1704 Coordinate about which to plot postage stamp plots. Defaults to
1705 ``"realization"``.
1706 overlay_cube: Cube | None, optional
1707 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1708 contour_cube: Cube | None, optional
1709 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1711 Raises
1712 ------
1713 ValueError
1714 If the cube doesn't have the right dimensions.
1715 TypeError
1716 If the cube isn't a single cube.
1717 """
1718 recipe_title = get_recipe_metadata().get("title", "Untitled")
1720 # Ensure we've got a single cube.
1721 cube = check_single_cube(cube)
1723 # Check if there is a valid stamp coordinate in cube dimensions.
1724 if stamp_coordinate == "realization": 1724 ↛ 1729line 1724 didn't jump to line 1729 because the condition on line 1724 was always true
1725 stamp_coordinate = check_stamp_coordinate(cube)
1727 # Make postage stamp plots if stamp_coordinate exists and has more than a
1728 # single point.
1729 plotting_func = _plot_and_save_spatial_plot
1730 try:
1731 if cube.coord(stamp_coordinate).shape[0] > 1:
1732 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1733 except iris.exceptions.CoordinateNotFoundError:
1734 pass
1736 # Produce a geographical scatter plot if the data have a
1737 # dimension called observation or model_obs_error
1738 if any( 1738 ↛ 1742line 1738 didn't jump to line 1742 because the condition on line 1738 was never true
1739 crd.var_name == "station" or crd.var_name == "model_obs_error"
1740 for crd in cube.coords()
1741 ):
1742 plotting_func = _plot_and_save_scattermap_plot
1744 # Must have a sequence coordinate.
1745 try:
1746 cube.coord(sequence_coordinate)
1747 except iris.exceptions.CoordinateNotFoundError as err:
1748 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1750 # Create a plot for each value of the sequence coordinate.
1751 plot_index = []
1752 nplot = np.size(cube.coord(sequence_coordinate).points)
1754 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1755 # Set plot titles and filename
1756 seq_coord = cube_slice.coord(sequence_coordinate)
1757 plot_title, plot_filename = _set_title_and_filename(
1758 seq_coord, nplot, recipe_title, filename
1759 )
1761 # Extract sequence slice for overlay_cube and contour_cube if required.
1762 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1763 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1765 # Do the actual plotting.
1766 plotting_func(
1767 cube_slice,
1768 filename=plot_filename,
1769 stamp_coordinate=stamp_coordinate,
1770 title=plot_title,
1771 method=method,
1772 overlay_cube=overlay_slice,
1773 contour_cube=contour_slice,
1774 **kwargs,
1775 )
1776 plot_index.append(plot_filename)
1778 # Add list of plots to plot metadata.
1779 complete_plot_index = _append_to_plot_index(plot_index)
1781 # Make a page to display the plots.
1782 _make_plot_html_page(complete_plot_index)
1785####################
1786# Public functions #
1787####################
1790def spatial_contour_plot(
1791 cube: iris.cube.Cube,
1792 filename: str = None,
1793 sequence_coordinate: str = "time",
1794 stamp_coordinate: str = "realization",
1795 **kwargs,
1796) -> iris.cube.Cube:
1797 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1799 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1800 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1801 is present then postage stamp plots will be produced.
1803 Parameters
1804 ----------
1805 cube: Cube
1806 Iris cube of the data to plot. It should have two spatial dimensions,
1807 such as lat and lon, and may also have a another two dimension to be
1808 plotted sequentially and/or as postage stamp plots.
1809 filename: str, optional
1810 Name of the plot to write, used as a prefix for plot sequences. Defaults
1811 to the recipe name.
1812 sequence_coordinate: str, optional
1813 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1814 This coordinate must exist in the cube.
1815 stamp_coordinate: str, optional
1816 Coordinate about which to plot postage stamp plots. Defaults to
1817 ``"realization"``.
1819 Returns
1820 -------
1821 Cube
1822 The original cube (so further operations can be applied).
1824 Raises
1825 ------
1826 ValueError
1827 If the cube doesn't have the right dimensions.
1828 TypeError
1829 If the cube isn't a single cube.
1830 """
1831 _spatial_plot(
1832 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1833 )
1834 return cube
1837def spatial_pcolormesh_plot(
1838 cube: iris.cube.Cube,
1839 filename: str = None,
1840 sequence_coordinate: str = "time",
1841 stamp_coordinate: str = "realization",
1842 **kwargs,
1843) -> iris.cube.Cube:
1844 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1846 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1847 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1848 is present then postage stamp plots will be produced.
1850 This function is significantly faster than ``spatial_contour_plot``,
1851 especially at high resolutions, and should be preferred unless contiguous
1852 contour areas are important.
1854 Parameters
1855 ----------
1856 cube: Cube
1857 Iris cube of the data to plot. It should have two spatial dimensions,
1858 such as lat and lon, and may also have a another two dimension to be
1859 plotted sequentially and/or as postage stamp plots.
1860 filename: str, optional
1861 Name of the plot to write, used as a prefix for plot sequences. Defaults
1862 to the recipe name.
1863 sequence_coordinate: str, optional
1864 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1865 This coordinate must exist in the cube.
1866 stamp_coordinate: str, optional
1867 Coordinate about which to plot postage stamp plots. Defaults to
1868 ``"realization"``.
1870 Returns
1871 -------
1872 Cube
1873 The original cube (so further operations can be applied).
1875 Raises
1876 ------
1877 ValueError
1878 If the cube doesn't have the right dimensions.
1879 TypeError
1880 If the cube isn't a single cube.
1881 """
1882 _spatial_plot(
1883 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1884 )
1885 return cube
1888def spatial_multi_pcolormesh_plot(
1889 cube: iris.cube.Cube,
1890 overlay_cube: iris.cube.Cube,
1891 contour_cube: iris.cube.Cube,
1892 filename: str = None,
1893 sequence_coordinate: str = "time",
1894 stamp_coordinate: str = "realization",
1895 **kwargs,
1896) -> iris.cube.Cube:
1897 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1899 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1900 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1901 is present then postage stamp plots will be produced.
1903 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1905 If specified, contours of a contour_cube can be overplotted on top of those.
1907 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1909 This function is significantly faster than ``spatial_contour_plot``,
1910 especially at high resolutions, and should be preferred unless contiguous
1911 contour areas are important.
1913 Parameters
1914 ----------
1915 cube: Cube
1916 Iris cube of the data to plot. It should have two spatial dimensions,
1917 such as lat and lon, and may also have a another two dimension to be
1918 plotted sequentially and/or as postage stamp plots.
1919 overlay_cube: Cube
1920 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1921 such as lat and lon, and may also have a another two dimension to be
1922 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.
1923 contour_cube: Cube
1924 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,
1925 such as lat and lon, and may also have a another two dimension to be
1926 plotted sequentially and/or as postage stamp plots.
1927 filename: str, optional
1928 Name of the plot to write, used as a prefix for plot sequences. Defaults
1929 to the recipe name.
1930 sequence_coordinate: str, optional
1931 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1932 This coordinate must exist in the cube.
1933 stamp_coordinate: str, optional
1934 Coordinate about which to plot postage stamp plots. Defaults to
1935 ``"realization"``.
1937 Returns
1938 -------
1939 Cube
1940 The original cube (so further operations can be applied).
1942 Raises
1943 ------
1944 ValueError
1945 If the cube doesn't have the right dimensions.
1946 TypeError
1947 If the cube isn't a single cube.
1948 """
1949 _spatial_plot(
1950 "pcolormesh",
1951 cube,
1952 filename,
1953 sequence_coordinate,
1954 stamp_coordinate,
1955 overlay_cube=overlay_cube,
1956 contour_cube=contour_cube,
1957 )
1958 return cube, overlay_cube, contour_cube
1961# TODO: Expand function to handle ensemble data.
1962# line_coordinate: str, optional
1963# Coordinate about which to plot multiple lines. Defaults to
1964# ``"realization"``.
1965def plot_line_series(
1966 cube: iris.cube.Cube | iris.cube.CubeList,
1967 filename: str = None,
1968 series_coordinate: str = "time",
1969 sequence_coordinate: str = "time",
1970 # add the following for ensembles
1971 stamp_coordinate: str = "realization",
1972 single_plot: bool = False,
1973 **kwargs,
1974) -> iris.cube.Cube | iris.cube.CubeList:
1975 """Plot a line plot for the specified coordinate.
1977 The Cube or CubeList must be 1D.
1979 Parameters
1980 ----------
1981 iris.cube | iris.cube.CubeList
1982 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
1983 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
1984 We do not support different data such as temperature and humidity in the same CubeList for plotting.
1985 filename: str, optional
1986 Name of the plot to write, used as a prefix for plot sequences. Defaults
1987 to the recipe name.
1988 series_coordinate: str, optional
1989 Coordinate about which to make a series. Defaults to ``"time"``. This
1990 coordinate must exist in the cube.
1992 Returns
1993 -------
1994 iris.cube.Cube | iris.cube.CubeList
1995 The original Cube or CubeList (so further operations can be applied).
1996 plotted data.
1998 Raises
1999 ------
2000 ValueError
2001 If the cubes don't have the right dimensions.
2002 TypeError
2003 If the cube isn't a Cube or CubeList.
2004 """
2005 # Ensure we have a name for the plot file.
2006 recipe_title = get_recipe_metadata().get("title", "Untitled")
2008 num_models = get_num_models(cube)
2010 validate_cube_shape(cube, num_models)
2012 # Iterate over all cubes and extract coordinate to plot.
2013 cubes = iter_maybe(cube)
2014 coords = []
2015 for cube in cubes:
2016 try:
2017 coords.append(cube.coord(series_coordinate))
2018 except iris.exceptions.CoordinateNotFoundError as err:
2019 raise ValueError(
2020 f"Cube must have a {series_coordinate} coordinate."
2021 ) from err
2022 if cube.coords("realization"): 2022 ↛ 2026line 2022 didn't jump to line 2026 because the condition on line 2022 was always true
2023 if cube.ndim > 3: 2023 ↛ 2024line 2023 didn't jump to line 2024 because the condition on line 2023 was never true
2024 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2025 else:
2026 raise ValueError("Cube must have a realization coordinate.")
2028 plot_index = []
2030 # Check if this is a spectral plot by looking for spectral coordinates
2031 is_spectral_plot = series_coordinate in [
2032 "frequency",
2033 "physical_wavenumber",
2034 "wavelength",
2035 ]
2037 if is_spectral_plot:
2038 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2039 # coordinate frequency/wavenumber.
2040 # If several power spectra are plotted with time as sequence_coordinate for the
2041 # time slider option.
2043 # Internal plotting function.
2044 plotting_func = _plot_and_save_line_power_spectrum_series
2046 for cube in cubes:
2047 try:
2048 cube.coord(sequence_coordinate)
2049 except iris.exceptions.CoordinateNotFoundError as err:
2050 raise ValueError(
2051 f"Cube must have a {sequence_coordinate} coordinate."
2052 ) from err
2054 if num_models == 1: 2054 ↛ 2068line 2054 didn't jump to line 2068 because the condition on line 2054 was always true
2055 # check for ensembles
2056 if ( 2056 ↛ 2060line 2056 didn't jump to line 2060 because the condition on line 2056 was never true
2057 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2058 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2059 ):
2060 if single_plot:
2061 # Plot spectra, mean and ensemble spread on 1 plot
2062 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2063 else:
2064 # Plot postage stamps
2065 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2066 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2067 else:
2068 all_points = sorted(
2069 set(
2070 itertools.chain.from_iterable(
2071 cb.coord(sequence_coordinate).points for cb in cubes
2072 )
2073 )
2074 )
2075 all_slices = list(
2076 itertools.chain.from_iterable(
2077 cb.slices_over(sequence_coordinate) for cb in cubes
2078 )
2079 )
2080 # Matched slices (matched by seq coord point; it may happen that
2081 # evaluated models do not cover the same seq coord range, hence matching
2082 # necessary)
2083 cube_iterables = [
2084 iris.cube.CubeList(
2085 s
2086 for s in all_slices
2087 if s.coord(sequence_coordinate).points[0] == point
2088 )
2089 for point in all_points
2090 ]
2092 nplot = np.size(cube.coord(sequence_coordinate).points)
2094 # Create a plot for each value of the sequence coordinate. Allowing for
2095 # multiple cubes in a CubeList to be plotted in the same plot for similar
2096 # sequence values. Passing a CubeList into the internal plotting function
2097 # for similar values of the sequence coordinate. cube_slice can be an
2098 # iris.cube.Cube or an iris.cube.CubeList.
2100 for cube_slice in cube_iterables:
2101 # Normalize cube_slice to a list of cubes
2102 if isinstance(cube_slice, iris.cube.CubeList): 2102 ↛ 2103line 2102 didn't jump to line 2103 because the condition on line 2102 was never true
2103 cubes = list(cube_slice)
2104 elif isinstance(cube_slice, iris.cube.Cube): 2104 ↛ 2107line 2104 didn't jump to line 2107 because the condition on line 2104 was always true
2105 cubes = [cube_slice]
2106 else:
2107 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2109 # Use sequence value so multiple sequences can merge.
2110 seq_coord = cube_slice[0].coord(sequence_coordinate)
2111 plot_title, plot_filename = _set_title_and_filename(
2112 seq_coord, nplot, recipe_title, filename
2113 )
2115 # Format the coordinate value in a unit appropriate way.
2116 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2118 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2119 if nplot == 1 and seq_coord.has_bounds: 2119 ↛ 2124line 2119 didn't jump to line 2124 because the condition on line 2119 was always true
2120 if np.size(seq_coord.bounds) > 1: 2120 ↛ 2121line 2120 didn't jump to line 2121 because the condition on line 2120 was never true
2121 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2123 # Do the actual plotting.
2124 plotting_func(
2125 cube_slice,
2126 coords,
2127 stamp_coordinate,
2128 plot_filename,
2129 title,
2130 series_coordinate,
2131 )
2133 plot_index.append(plot_filename)
2134 else:
2135 # Format the title and filename using plotted series coordinate
2136 nplot = 1
2137 seq_coord = coords[0]
2138 plot_title, plot_filename = _set_title_and_filename(
2139 seq_coord, nplot, recipe_title, filename
2140 )
2141 # Do the actual plotting for all other series coordinate options.
2142 _plot_and_save_line_series(
2143 cubes, coords, stamp_coordinate, plot_filename, plot_title
2144 )
2146 plot_index.append(plot_filename)
2148 # append plot to list of plots
2149 complete_plot_index = _append_to_plot_index(plot_index)
2151 # Make a page to display the plots.
2152 _make_plot_html_page(complete_plot_index)
2154 return cube
2157def plot_vertical_line_series(
2158 cubes: iris.cube.Cube | iris.cube.CubeList,
2159 filename: str = None,
2160 series_coordinate: str = "model_level_number",
2161 sequence_coordinate: str = "time",
2162 # line_coordinate: str = "realization",
2163 **kwargs,
2164) -> iris.cube.Cube | iris.cube.CubeList:
2165 """Plot a line plot against a type of vertical coordinate.
2167 The Cube or CubeList must be 1D.
2169 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2170 then a sequence of plots will be produced.
2172 Parameters
2173 ----------
2174 iris.cube | iris.cube.CubeList
2175 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2176 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2177 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2178 filename: str, optional
2179 Name of the plot to write, used as a prefix for plot sequences. Defaults
2180 to the recipe name.
2181 series_coordinate: str, optional
2182 Coordinate to plot on the y-axis. Can be ``pressure`` or
2183 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2184 for LFRic. Defaults to ``model_level_number``.
2185 This coordinate must exist in the cube.
2186 sequence_coordinate: str, optional
2187 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2188 This coordinate must exist in the cube.
2190 Returns
2191 -------
2192 iris.cube.Cube | iris.cube.CubeList
2193 The original Cube or CubeList (so further operations can be applied).
2194 Plotted data.
2196 Raises
2197 ------
2198 ValueError
2199 If the cubes doesn't have the right dimensions.
2200 TypeError
2201 If the cube isn't a Cube or CubeList.
2202 """
2203 # Ensure we have a name for the plot file.
2204 recipe_title = get_recipe_metadata().get("title", "Untitled")
2206 cubes = iter_maybe(cubes)
2207 # Initialise empty list to hold all data from all cubes in a CubeList
2208 all_data = []
2210 # Store min/max ranges for x range.
2211 x_levels = []
2213 num_models = get_num_models(cubes)
2215 validate_cube_shape(cubes, num_models)
2217 # Iterate over all cubes in cube or CubeList and plot.
2218 coords = []
2219 for cube in cubes:
2220 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2221 try:
2222 coords.append(cube.coord(series_coordinate))
2223 except iris.exceptions.CoordinateNotFoundError as err:
2224 raise ValueError(
2225 f"Cube must have a {series_coordinate} coordinate."
2226 ) from err
2228 try:
2229 if cube.ndim > 1 or not cube.coords("realization"): 2229 ↛ 2237line 2229 didn't jump to line 2237 because the condition on line 2229 was always true
2230 cube.coord(sequence_coordinate)
2231 except iris.exceptions.CoordinateNotFoundError as err:
2232 raise ValueError(
2233 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2234 ) from err
2236 # Get minimum and maximum from levels information.
2237 _, levels, _ = colorbar_map_levels(cube, axis="x")
2238 if levels is not None: 2238 ↛ 2242line 2238 didn't jump to line 2242 because the condition on line 2238 was always true
2239 x_levels.append(min(levels))
2240 x_levels.append(max(levels))
2241 else:
2242 all_data.append(cube.data)
2244 if len(x_levels) == 0: 2244 ↛ 2246line 2244 didn't jump to line 2246 because the condition on line 2244 was never true
2245 # Combine all data into a single NumPy array
2246 combined_data = np.concatenate(all_data)
2248 # Set the lower and upper limit for the x-axis to ensure all plots have
2249 # same range. This needs to read the whole cube over the range of the
2250 # sequence and if applicable postage stamp coordinate.
2251 vmin = np.floor(combined_data.min())
2252 vmax = np.ceil(combined_data.max())
2253 else:
2254 vmin = min(x_levels)
2255 vmax = max(x_levels)
2257 # Matching the slices (matching by seq coord point; it may happen that
2258 # evaluated models do not cover the same seq coord range, hence matching
2259 # necessary)
2260 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2262 # Create a plot for each value of the sequence coordinate.
2263 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2264 # similar sequence values. Passing a CubeList into the internal plotting function
2265 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2266 # or an iris.cube.CubeList.
2267 plot_index = []
2268 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2269 for cubes_slice in cube_iterables:
2270 # Format the coordinate value in a unit appropriate way.
2271 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2272 plot_title, plot_filename = _set_title_and_filename(
2273 seq_coord, nplot, recipe_title, filename
2274 )
2276 # Do the actual plotting.
2277 _plot_and_save_vertical_line_series(
2278 cubes_slice,
2279 coords,
2280 "realization",
2281 plot_filename,
2282 series_coordinate,
2283 title=plot_title,
2284 vmin=vmin,
2285 vmax=vmax,
2286 )
2287 plot_index.append(plot_filename)
2289 # Add list of plots to plot metadata.
2290 complete_plot_index = _append_to_plot_index(plot_index)
2292 # Make a page to display the plots.
2293 _make_plot_html_page(complete_plot_index)
2295 return cubes
2298def qq_plot(
2299 cubes: iris.cube.CubeList,
2300 coordinates: list[str],
2301 percentiles: list[float],
2302 model_names: list[str],
2303 filename: str = None,
2304 one_to_one: bool = True,
2305 **kwargs,
2306) -> iris.cube.CubeList:
2307 """Plot a Quantile-Quantile plot between two models for common time points.
2309 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2310 collapsed within the operator over all specified coordinates such as
2311 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2312 ensemble members to ensure a 1D cube (array).
2314 Parameters
2315 ----------
2316 cubes: iris.cube.CubeList
2317 Two cubes of the same variable with different models.
2318 coordinate: list[str]
2319 The list of coordinates to collapse over. This list should be
2320 every coordinate within the cube to result in a 1D cube around
2321 the percentile coordinate.
2322 percent: list[float]
2323 A list of percentiles to appear in the plot.
2324 model_names: list[str]
2325 A list of model names to appear on the axis of the plot.
2326 filename: str, optional
2327 Filename of the plot to write.
2328 one_to_one: bool, optional
2329 If True a 1:1 line is plotted; if False it is not. Default is True.
2331 Raises
2332 ------
2333 ValueError
2334 When the cubes are not compatible.
2336 Notes
2337 -----
2338 The quantile-quantile plot is a variant on the scatter plot representing
2339 two datasets by their quantiles (percentiles) for common time points.
2340 This plot does not use a theoretical distribution to compare against, but
2341 compares percentiles of two datasets. This plot does
2342 not use all raw data points, but plots the selected percentiles (quantiles) of
2343 each variable instead for the two datasets, thereby normalising the data for a
2344 direct comparison between the selected percentiles of the two dataset distributions.
2346 Quantile-quantile plots are valuable for comparing against
2347 observations and other models. Identical percentiles between the variables
2348 will lie on the one-to-one line implying the values correspond well to each
2349 other. Where there is a deviation from the one-to-one line a range of
2350 possibilities exist depending on how and where the data is shifted (e.g.,
2351 Wilks 2011 [Wilks2011]_).
2353 For distributions above the one-to-one line the distribution is left-skewed;
2354 below is right-skewed. A distinct break implies a bimodal distribution, and
2355 closer values/values further apart at the tails imply poor representation of
2356 the extremes.
2358 References
2359 ----------
2360 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2361 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2362 """
2363 # Check cubes using same functionality as the difference operator.
2364 if len(cubes) != 2:
2365 raise ValueError("cubes should contain exactly 2 cubes.")
2366 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2367 other: Cube = cubes.extract_cube(
2368 iris.Constraint(
2369 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2370 )
2371 )
2373 # Get spatial coord names.
2374 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2375 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2377 # Ensure cubes to compare are on common differencing grid.
2378 # This is triggered if either
2379 # i) latitude and longitude shapes are not the same. Note grid points
2380 # are not compared directly as these can differ through rounding
2381 # errors.
2382 # ii) or variables are known to often sit on different grid staggering
2383 # in different models (e.g. cell center vs cell edge), as is the case
2384 # for UM and LFRic comparisons.
2385 # In future greater choice of regridding method might be applied depending
2386 # on variable type. Linear regridding can in general be appropriate for smooth
2387 # variables. Care should be taken with interpretation of differences
2388 # given this dependency on regridding.
2389 if (
2390 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2391 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2392 ) or (
2393 base.long_name
2394 in [
2395 "eastward_wind_at_10m",
2396 "northward_wind_at_10m",
2397 "northward_wind_at_cell_centres",
2398 "eastward_wind_at_cell_centres",
2399 "zonal_wind_at_pressure_levels",
2400 "meridional_wind_at_pressure_levels",
2401 "potential_vorticity_at_pressure_levels",
2402 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2403 ]
2404 ):
2405 logging.debug(
2406 "Linear regridding base cube to other grid to compute differences"
2407 )
2408 base = regrid_onto_cube(base, other, method="Linear")
2410 # Extract just common time points.
2411 base, other = _extract_common_time_points(base, other)
2413 # Equalise attributes so we can merge.
2414 fully_equalise_attributes([base, other])
2415 logging.debug("Base: %s\nOther: %s", base, other)
2417 # Collapse cubes.
2418 base = collapse(
2419 base,
2420 coordinate=coordinates,
2421 method="PERCENTILE",
2422 additional_percent=percentiles,
2423 )
2424 other = collapse(
2425 other,
2426 coordinate=coordinates,
2427 method="PERCENTILE",
2428 additional_percent=percentiles,
2429 )
2431 # Ensure we have a name for the plot file.
2432 recipe_title = get_recipe_metadata().get("title", "Untitled")
2433 title = f"{recipe_title}"
2435 if filename is None:
2436 filename = slugify(recipe_title)
2438 # Add file extension.
2439 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2441 # Do the actual plotting on a scatter plot
2442 _plot_and_save_scatter_plot(
2443 base, other, plot_filename, title, one_to_one, model_names
2444 )
2446 # Add list of plots to plot metadata.
2447 plot_index = _append_to_plot_index([plot_filename])
2449 # Make a page to display the plots.
2450 _make_plot_html_page(plot_index)
2452 return iris.cube.CubeList([base, other])
2455def scatter_plot(
2456 cube_x: iris.cube.Cube | iris.cube.CubeList,
2457 cube_y: iris.cube.Cube | iris.cube.CubeList,
2458 filename: str = None,
2459 one_to_one: bool = True,
2460 **kwargs,
2461) -> iris.cube.CubeList:
2462 """Plot a scatter plot between two variables.
2464 Both cubes must be 1D.
2466 Parameters
2467 ----------
2468 cube_x: Cube | CubeList
2469 1 dimensional Cube of the data to plot on y-axis.
2470 cube_y: Cube | CubeList
2471 1 dimensional Cube of the data to plot on x-axis.
2472 filename: str, optional
2473 Filename of the plot to write.
2474 one_to_one: bool, optional
2475 If True a 1:1 line is plotted; if False it is not. Default is True.
2477 Returns
2478 -------
2479 cubes: CubeList
2480 CubeList of the original x and y cubes for further processing.
2482 Raises
2483 ------
2484 ValueError
2485 If the cube doesn't have the right dimensions and cubes not the same
2486 size.
2487 TypeError
2488 If the cube isn't a single cube.
2490 Notes
2491 -----
2492 Scatter plots are used for determining if there is a relationship between
2493 two variables. Positive relations have a slope going from bottom left to top
2494 right; Negative relations have a slope going from top left to bottom right.
2495 """
2496 # Iterate over all cubes in cube or CubeList and plot.
2497 for cube_iter in iter_maybe(cube_x):
2498 # Check cubes are correct shape.
2499 cube_iter = check_single_cube(cube_iter)
2500 if cube_iter.ndim > 1:
2501 raise ValueError("cube_x must be 1D.")
2503 # Iterate over all cubes in cube or CubeList and plot.
2504 for cube_iter in iter_maybe(cube_y):
2505 # Check cubes are correct shape.
2506 cube_iter = check_single_cube(cube_iter)
2507 if cube_iter.ndim > 1:
2508 raise ValueError("cube_y must be 1D.")
2510 # Ensure we have a name for the plot file.
2511 recipe_title = get_recipe_metadata().get("title", "Untitled")
2512 title = f"{recipe_title}"
2514 if filename is None:
2515 filename = slugify(recipe_title)
2517 # Add file extension.
2518 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2520 # Do the actual plotting.
2521 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2523 # Add list of plots to plot metadata.
2524 plot_index = _append_to_plot_index([plot_filename])
2526 # Make a page to display the plots.
2527 _make_plot_html_page(plot_index)
2529 return iris.cube.CubeList([cube_x, cube_y])
2532def vector_plot(
2533 cube_u: iris.cube.Cube,
2534 cube_v: iris.cube.Cube,
2535 filename: str = None,
2536 sequence_coordinate: str = "time",
2537 **kwargs,
2538) -> iris.cube.CubeList:
2539 """Plot a vector plot based on the input u and v components."""
2540 recipe_title = get_recipe_metadata().get("title", "Untitled")
2542 # Cubes must have a matching sequence coordinate.
2543 try:
2544 # Check that the u and v cubes have the same sequence coordinate.
2545 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2545 ↛ anywhereline 2545 didn't jump anywhere: it always raised an exception.
2546 raise ValueError("Coordinates do not match.")
2547 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2548 raise ValueError(
2549 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2550 ) from err
2552 # Create a plot for each value of the sequence coordinate.
2553 plot_index = []
2554 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2555 for cube_u_slice, cube_v_slice in zip(
2556 cube_u.slices_over(sequence_coordinate),
2557 cube_v.slices_over(sequence_coordinate),
2558 strict=True,
2559 ):
2560 # Format the coordinate value in a unit appropriate way.
2561 seq_coord = cube_u_slice.coord(sequence_coordinate)
2562 plot_title, plot_filename = _set_title_and_filename(
2563 seq_coord, nplot, recipe_title, filename
2564 )
2566 # Do the actual plotting.
2567 _plot_and_save_vector_plot(
2568 cube_u_slice,
2569 cube_v_slice,
2570 filename=plot_filename,
2571 title=plot_title,
2572 method="contourf",
2573 )
2574 plot_index.append(plot_filename)
2576 # Add list of plots to plot metadata.
2577 complete_plot_index = _append_to_plot_index(plot_index)
2579 # Make a page to display the plots.
2580 _make_plot_html_page(complete_plot_index)
2582 return iris.cube.CubeList([cube_u, cube_v])
2585def plot_histogram_series(
2586 cubes: iris.cube.Cube | iris.cube.CubeList,
2587 filename: str = None,
2588 sequence_coordinate: str = "time",
2589 stamp_coordinate: str = "realization",
2590 single_plot: bool = False,
2591 **kwargs,
2592) -> iris.cube.Cube | iris.cube.CubeList:
2593 """Plot a histogram plot for each vertical level provided.
2595 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2596 is present then a sequence of plots will be produced using the time slider
2597 functionality to scroll through histograms against time. If a
2598 stamp_coordinate is present then postage stamp plots will be produced. If
2599 stamp_coordinate and single_plot is True, all postage stamp plots will be
2600 plotted in a single plot instead of separate postage stamp plots.
2602 Parameters
2603 ----------
2604 cubes: Cube | iris.cube.CubeList
2605 Iris cube or CubeList of the data to plot. It should have a single dimension other
2606 than the stamp coordinate.
2607 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2608 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2609 filename: str, optional
2610 Name of the plot to write, used as a prefix for plot sequences. Defaults
2611 to the recipe name.
2612 sequence_coordinate: str, optional
2613 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2614 This coordinate must exist in the cube and will be used for the time
2615 slider.
2616 stamp_coordinate: str, optional
2617 Coordinate about which to plot postage stamp plots. Defaults to
2618 ``"realization"``.
2619 single_plot: bool, optional
2620 If True, all postage stamp plots will be plotted in a single plot. If
2621 False, each postage stamp plot will be plotted separately. Is only valid
2622 if stamp_coordinate exists and has more than a single point.
2624 Returns
2625 -------
2626 iris.cube.Cube | iris.cube.CubeList
2627 The original Cube or CubeList (so further operations can be applied).
2628 Plotted data.
2630 Raises
2631 ------
2632 ValueError
2633 If the cube doesn't have the right dimensions.
2634 TypeError
2635 If the cube isn't a Cube or CubeList.
2636 """
2637 recipe_title = get_recipe_metadata().get("title", "Untitled")
2639 cubes = iter_maybe(cubes)
2640 # Ensure we have a name for the plot file.
2641 if filename is None:
2642 filename = slugify(recipe_title)
2644 # Internal plotting function.
2645 plotting_func = _plot_and_save_histogram_series
2647 num_models = get_num_models(cubes)
2649 validate_cube_shape(cubes, num_models)
2651 # If several histograms are plotted, check sequence_coordinate
2652 check_sequence_coordinate(cubes, sequence_coordinate)
2654 # Get axis minimum and maximum from levels information.
2655 # If no levels set, derive minima and maxima from data in CubeList.
2656 vmin, vmax = _set_axis_range(cubes)
2658 # Make postage stamp plots if stamp_coordinate exists and has more than a
2659 # single point. If single_plot is True:
2660 # -- all postage stamp plots will be plotted in a single plot instead of
2661 # separate postage stamp plots.
2662 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2663 # produced per single model only
2664 if num_models == 1: 2664 ↛ 2677line 2664 didn't jump to line 2677 because the condition on line 2664 was always true
2665 if ( 2665 ↛ 2669line 2665 didn't jump to line 2669 because the condition on line 2665 was never true
2666 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2667 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2668 ):
2669 if single_plot:
2670 plotting_func = (
2671 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2672 )
2673 else:
2674 plotting_func = _plot_and_save_postage_stamp_histogram_series
2675 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2676 else:
2677 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2679 plot_index = []
2680 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2681 # Create a plot for each value of the sequence coordinate. Allowing for
2682 # multiple cubes in a CubeList to be plotted in the same plot for similar
2683 # sequence values. Passing a CubeList into the internal plotting function
2684 # for similar values of the sequence coordinate. cube_slice can be an
2685 # iris.cube.Cube or an iris.cube.CubeList.
2686 for cube_slice in cube_iterables:
2687 single_cube = cube_slice
2688 if isinstance(cube_slice, iris.cube.CubeList): 2688 ↛ 2689line 2688 didn't jump to line 2689 because the condition on line 2688 was never true
2689 single_cube = cube_slice[0]
2691 # Ensure valid stamp coordinate in cube dimensions
2692 if stamp_coordinate == "realization": 2692 ↛ 2695line 2692 didn't jump to line 2695 because the condition on line 2692 was always true
2693 stamp_coordinate = check_stamp_coordinate(single_cube)
2694 # Set plot titles and filename, based on sequence coordinate
2695 seq_coord = single_cube.coord(sequence_coordinate)
2696 # Use time coordinate in title and filename if single histogram output.
2697 if sequence_coordinate == "realization" and nplot == 1: 2697 ↛ 2698line 2697 didn't jump to line 2698 because the condition on line 2697 was never true
2698 seq_coord = single_cube.coord("time")
2699 plot_title, plot_filename = _set_title_and_filename(
2700 seq_coord, nplot, recipe_title, filename
2701 )
2703 # Do the actual plotting.
2704 plotting_func(
2705 cube_slice,
2706 filename=plot_filename,
2707 stamp_coordinate=stamp_coordinate,
2708 title=plot_title,
2709 vmin=vmin,
2710 vmax=vmax,
2711 )
2712 plot_index.append(plot_filename)
2714 # Add list of plots to plot metadata.
2715 complete_plot_index = _append_to_plot_index(plot_index)
2717 # Make a page to display the plots.
2718 _make_plot_html_page(complete_plot_index)
2720 return cubes
2723def _plot_and_save_postage_stamp_power_spectrum_series(
2724 cubes: iris.cube.Cube,
2725 coords: list[iris.coords.Coord],
2726 stamp_coordinate: str,
2727 filename: str,
2728 title: str,
2729 series_coordinate: str = None,
2730 **kwargs,
2731):
2732 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2734 Parameters
2735 ----------
2736 cubes: Cube or CubeList
2737 Cube or Cubelist of the power spectrum data.
2738 coords: list[Coord]
2739 Coordinates to plot on the x-axis, one per cube.
2740 stamp_coordinate: str
2741 Coordinate that becomes different plots.
2742 filename: str
2743 Filename of the plot to write.
2744 title: str
2745 Plot title.
2746 series_coordinate: str, optional
2747 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2749 """
2750 # Use the smallest square grid that will fit the members.
2751 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
2753 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
2754 model_colors_map = get_model_colors_map(cubes)
2755 # ax = plt.gca()
2756 # Make a subplot for each member.
2757 for member, subplot in zip(
2758 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
2759 ):
2760 ax = plt.subplot(grid_size, grid_size, subplot)
2762 # Store min/max ranges.
2763 y_levels = []
2765 line_marker = None
2766 line_width = 1
2768 for cube in iter_maybe(member):
2769 xcoord = _select_series_coord(cube, series_coordinate)
2770 xname = xcoord.points
2772 yfield = cube.data # power spectrum
2773 label = None
2774 color = "black"
2775 if model_colors_map: 2775 ↛ 2776line 2775 didn't jump to line 2776 because the condition on line 2775 was never true
2776 label = cube.attributes.get("model_name")
2777 color = model_colors_map.get(label)
2779 if member.coord(stamp_coordinate).points == [0]:
2780 ax.plot(
2781 xname,
2782 yfield,
2783 color=color,
2784 marker=line_marker,
2785 ls="-",
2786 lw=line_width,
2787 label=f"{label} (control)"
2788 if len(cube.coord(stamp_coordinate).points) > 1
2789 else label,
2790 )
2791 # Label with member if part of an ensemble and not the control.
2792 else:
2793 ax.plot(
2794 xname,
2795 yfield,
2796 color=color,
2797 ls="-",
2798 lw=1.5,
2799 alpha=0.75,
2800 label=f"{label} (member)",
2801 )
2803 # Calculate the global min/max if multiple cubes are given.
2804 _, levels, _ = colorbar_map_levels(cube, axis="y")
2805 if levels is not None: 2805 ↛ 2806line 2805 didn't jump to line 2806 because the condition on line 2805 was never true
2806 y_levels.append(min(levels))
2807 y_levels.append(max(levels))
2809 # Add some labels and tweak the style.
2810 title = f"{title}"
2811 ax.set_title(title, fontsize=16)
2813 # Set appropriate x-axis label based on coordinate
2814 if series_coordinate == "wavelength" or ( 2814 ↛ 2817line 2814 didn't jump to line 2817 because the condition on line 2814 was never true
2815 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
2816 ):
2817 ax.set_xlabel("Wavelength (km)", fontsize=14)
2818 elif series_coordinate == "physical_wavenumber" or ( 2818 ↛ 2823line 2818 didn't jump to line 2823 because the condition on line 2818 was always true
2819 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
2820 ):
2821 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2822 else: # frequency or check units
2823 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
2824 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2825 else:
2826 ax.set_xlabel("Wavenumber", fontsize=14)
2828 ax.set_ylabel("Power Spectral Density", fontsize=14)
2829 ax.tick_params(axis="both", labelsize=12)
2831 # Set log-log scale
2832 ax.set_xscale("log")
2833 ax.set_yscale("log")
2835 # Add gridlines
2836 ax.grid(linestyle="--", color="grey", linewidth=1)
2837 # Ientify unique labels for legend
2838 handles = list(
2839 {
2840 label: handle
2841 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
2842 }.values()
2843 )
2844 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
2846 ax = plt.gca()
2847 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
2849 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
2850 logging.info("Saved histogram postage stamp plot to %s", filename)
2851 plt.close(fig)
2854def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
2855 cubes: iris.cube.Cube,
2856 coords: list[iris.coords.Coord],
2857 stamp_coordinate: str,
2858 filename: str,
2859 title: str,
2860 series_coordinate: str = None,
2861 **kwargs,
2862):
2863 """Plot and save power spectra for ensemble members in single plot.
2865 Parameters
2866 ----------
2867 cubes: Cube or CubeList
2868 Cube or Cubelist of the power spectrum data.
2869 coords: list[Coord]
2870 Coordinates to plot on the x-axis, one per cube.
2871 stamp_coordinate: str
2872 Coordinate that becomes different plots.
2873 filename: str
2874 Filename of the plot to write.
2875 title: str
2876 Plot title.
2877 series_coordinate: str, optional
2878 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2880 """
2881 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
2882 model_colors_map = get_model_colors_map(cubes)
2884 line_marker = None
2885 line_width = 1
2887 # Compute ensemble statistics to show spread
2888 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
2889 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
2890 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
2892 xcoord_global = mean_cube.coord(series_coordinate)
2893 x_global = xcoord_global.points
2895 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
2896 xcoord = _select_series_coord(member, series_coordinate)
2897 xname = xcoord.points
2899 yfield = member.data # power spectrum
2900 color = "black"
2901 if model_colors_map: 2901 ↛ 2905line 2901 didn't jump to line 2905 because the condition on line 2901 was always true
2902 label = member.attributes.get("model_name") if i == 0 else None
2903 color = model_colors_map.get(label)
2905 if member.coord(stamp_coordinate).points == [0]:
2906 ax.plot(
2907 xname,
2908 yfield,
2909 color=color,
2910 marker=line_marker,
2911 ls="-",
2912 lw=line_width,
2913 label=f"{label} (control)"
2914 if len(member.coord(stamp_coordinate).points) > 1
2915 else label,
2916 )
2917 # Label with member number if part of an ensemble and not the control.
2918 else:
2919 ax.plot(
2920 xname,
2921 yfield,
2922 color=color,
2923 ls="-",
2924 lw=1.5,
2925 alpha=0.75,
2926 label=label,
2927 )
2929 # Set appropriate x-axis label based on coordinate
2930 if series_coordinate == "wavelength" or ( 2930 ↛ 2933line 2930 didn't jump to line 2933 because the condition on line 2930 was never true
2931 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
2932 ):
2933 ax.set_xlabel("Wavelength (km)", fontsize=14)
2934 elif series_coordinate == "physical_wavenumber" or ( 2934 ↛ 2939line 2934 didn't jump to line 2939 because the condition on line 2934 was always true
2935 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
2936 ):
2937 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2938 else: # frequency or check units
2939 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
2940 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2941 else:
2942 ax.set_xlabel("Wavenumber", fontsize=14)
2944 # Add ensemble spread shading
2945 ax.fill_between(
2946 x_global,
2947 min_cube.data,
2948 max_cube.data,
2949 color="grey",
2950 alpha=0.3,
2951 label="Ensemble spread",
2952 )
2954 # Add ensemble mean line
2955 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
2957 ax.set_ylabel("Power Spectral Density", fontsize=14)
2958 ax.tick_params(axis="both", labelsize=12)
2960 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
2961 # Set log-log scale
2962 ax.set_xscale("log")
2963 ax.set_yscale("log")
2965 # Add gridlines
2966 ax.grid(linestyle="--", color="grey", linewidth=1)
2967 # Identify unique labels for legend
2968 handles = list(
2969 {
2970 label: handle
2971 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
2972 }.values()
2973 )
2974 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
2976 # Add a legend
2977 ax.legend(fontsize=16)
2979 # Figure title.
2980 ax.set_title(title, fontsize=16)
2982 # Save the figure to a file
2983 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
2985 # Close the figure
2986 plt.close(fig)