Coverage for src/CSET/operators/plot.py: 78%
929 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-02 14:04 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
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")
1303 # Create a cube containing the magnitude of the vector field.
1304 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1305 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1306 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1307 cube_vec_mag.rename(
1308 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1309 )
1311 # Specify the color bar
1312 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1314 # Setup plot map projection, extent and coastlines and borderlines.
1315 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1317 if method == "contourf":
1318 # Filled contour plot of the field.
1319 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1320 elif method == "pcolormesh":
1321 try:
1322 vmin = min(levels)
1323 vmax = max(levels)
1324 except TypeError:
1325 vmin, vmax = None, None
1326 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1327 # if levels are defined.
1328 if norm is not None:
1329 vmin = None
1330 vmax = None
1331 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1332 else:
1333 raise ValueError(f"Unknown plotting method: {method}")
1335 # Check to see if transect, and if so, adjust y axis.
1336 if is_transect(cube_vec_mag):
1337 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1338 axes.invert_yaxis()
1339 axes.set_yscale("log")
1340 axes.set_ylim(1100, 100)
1341 # If both model_level_number and level_height exists, iplt can construct
1342 # plot as a function of height above orography (NOT sea level).
1343 elif {"model_level_number", "level_height"}.issubset(
1344 {coord.name() for coord in cube_vec_mag.coords()}
1345 ):
1346 axes.set_yscale("log")
1348 axes.set_title(
1349 f"{title}\n"
1350 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1351 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1352 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1353 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1354 fontsize=16,
1355 )
1357 else:
1358 # Add title.
1359 axes.set_title(title, fontsize=16)
1361 # Add watermark with min/max/mean. Currently not user togglable.
1362 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1363 axes.annotate(
1364 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}",
1365 xy=(0.05, -0.05),
1366 xycoords="axes fraction",
1367 xytext=(-5, 5),
1368 textcoords="offset points",
1369 ha="right",
1370 va="bottom",
1371 size=11,
1372 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1373 )
1375 # Add colour bar.
1376 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1377 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1378 # add ticks and tick_labels for every levels if less than 20 levels exist
1379 if levels is not None and len(levels) < 20:
1380 cbar.set_ticks(levels)
1381 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1383 # 30 barbs along the longest axis of the plot, or a barb per point for data
1384 # with less than 30 points.
1385 step = max(max(cube_u.shape) // 30, 1)
1386 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1388 # Save plot.
1389 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1390 logging.info("Saved vector plot to %s", filename)
1391 plt.close(fig)
1394def _plot_and_save_histogram_series(
1395 cubes: iris.cube.Cube | iris.cube.CubeList,
1396 filename: str,
1397 title: str,
1398 vmin: float,
1399 vmax: float,
1400 **kwargs,
1401):
1402 """Plot and save a histogram series.
1404 Parameters
1405 ----------
1406 cubes: Cube or CubeList
1407 2 dimensional Cube or CubeList of the data to plot as histogram.
1408 filename: str
1409 Filename of the plot to write.
1410 title: str
1411 Plot title.
1412 vmin: float
1413 minimum for colorbar
1414 vmax: float
1415 maximum for colorbar
1416 """
1417 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1418 ax = plt.gca()
1420 model_colors_map = get_model_colors_map(cubes)
1422 # Set default that histograms will produce probability density function
1423 # at each bin (integral over range sums to 1).
1424 density = True
1426 for cube in iter_maybe(cubes):
1427 # Easier to check title (where var name originates)
1428 # than seeing if long names exist etc.
1429 # Exception case, where distribution better fits log scales/bins.
1430 if "surface_microphysical" in title:
1431 if "amount" in title: 1431 ↛ 1433line 1431 didn't jump to line 1433 because the condition on line 1431 was never true
1432 # Compute histogram following Klingaman et al. (2017): ASoP
1433 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1434 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1435 density = False
1436 else:
1437 bins = 10.0 ** (
1438 np.arange(-10, 27, 1) / 10.0
1439 ) # Suggestion from RMED toolbox.
1440 bins = np.insert(bins, 0, 0)
1441 ax.set_yscale("log")
1442 vmin = bins[1]
1443 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1444 ax.set_xscale("log")
1445 elif "lightning" in title:
1446 bins = [0, 1, 2, 3, 4, 5]
1447 else:
1448 bins = np.linspace(vmin, vmax, 51)
1449 logging.debug(
1450 "Plotting histogram with %s bins %s - %s.",
1451 np.size(bins),
1452 np.min(bins),
1453 np.max(bins),
1454 )
1456 # Reshape cube data into a single array to allow for a single histogram.
1457 # Otherwise we plot xdim histograms stacked.
1458 cube_data_1d = (cube.data).flatten()
1460 label = None
1461 color = "black"
1462 if model_colors_map: 1462 ↛ 1463line 1462 didn't jump to line 1463 because the condition on line 1462 was never true
1463 label = cube.attributes.get("model_name")
1464 color = model_colors_map[label]
1465 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1467 # Compute area under curve.
1468 if "surface_microphysical" in title and "amount" in title: 1468 ↛ 1469line 1468 didn't jump to line 1469 because the condition on line 1468 was never true
1469 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1470 x = x * bin_mean / x.sum()
1471 x = x[1:]
1472 y = y[1:]
1474 ax.plot(
1475 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1476 )
1478 # Add some labels and tweak the style.
1479 ax.set_title(title, fontsize=16)
1480 ax.set_xlabel(
1481 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1482 )
1483 ax.set_ylabel("Normalised probability density", fontsize=14)
1484 if "surface_microphysical" in title and "amount" in title: 1484 ↛ 1485line 1484 didn't jump to line 1485 because the condition on line 1484 was never true
1485 ax.set_ylabel(
1486 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1487 )
1488 ax.set_xlim(vmin, vmax)
1489 ax.tick_params(axis="both", labelsize=12)
1491 # Overlay grid-lines onto histogram plot.
1492 ax.grid(linestyle="--", color="grey", linewidth=1)
1493 if model_colors_map: 1493 ↛ 1494line 1493 didn't jump to line 1494 because the condition on line 1493 was never true
1494 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1496 # Save plot.
1497 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1498 logging.info("Saved histogram plot to %s", filename)
1499 plt.close(fig)
1502def _plot_and_save_postage_stamp_histogram_series(
1503 cube: iris.cube.Cube,
1504 filename: str,
1505 title: str,
1506 stamp_coordinate: str,
1507 vmin: float,
1508 vmax: float,
1509 **kwargs,
1510):
1511 """Plot and save postage (ensemble members) stamps for a histogram series.
1513 Parameters
1514 ----------
1515 cube: Cube
1516 2 dimensional Cube of the data to plot as histogram.
1517 filename: str
1518 Filename of the plot to write.
1519 title: str
1520 Plot title.
1521 stamp_coordinate: str
1522 Coordinate that becomes different plots.
1523 vmin: float
1524 minimum for pdf x-axis
1525 vmax: float
1526 maximum for pdf x-axis
1527 """
1528 # Use the smallest square grid that will fit the members.
1529 nmember = len(cube.coord(stamp_coordinate).points)
1530 grid_rows = int(math.sqrt(nmember))
1531 grid_size = math.ceil(nmember / grid_rows)
1533 fig = plt.figure(
1534 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1535 )
1536 # Make a subplot for each member.
1537 for member, subplot in zip(
1538 cube.slices_over(stamp_coordinate),
1539 range(1, grid_size * grid_rows + 1),
1540 strict=False,
1541 ):
1542 # Implicit interface is much easier here, due to needing to have the
1543 # cartopy GeoAxes generated.
1544 plt.subplot(grid_rows, grid_size, subplot)
1545 # Reshape cube data into a single array to allow for a single histogram.
1546 # Otherwise we plot xdim histograms stacked.
1547 member_data_1d = (member.data).flatten()
1548 plt.hist(member_data_1d, density=True, stacked=True)
1549 axes = plt.gca()
1550 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1551 axes.set_title(f"{mtitle}")
1552 axes.set_xlim(vmin, vmax)
1554 # Overall figure title.
1555 fig.suptitle(title, fontsize=16)
1557 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1558 logging.info("Saved histogram postage stamp plot to %s", filename)
1559 plt.close(fig)
1562def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1563 cube: iris.cube.Cube,
1564 filename: str,
1565 title: str,
1566 stamp_coordinate: str,
1567 vmin: float,
1568 vmax: float,
1569 **kwargs,
1570):
1571 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1572 ax.set_title(title, fontsize=16)
1573 ax.set_xlim(vmin, vmax)
1574 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1575 ax.set_ylabel("normalised probability density", fontsize=14)
1576 # Loop over all slices along the stamp_coordinate
1577 for member in cube.slices_over(stamp_coordinate):
1578 # Flatten the member data to 1D
1579 member_data_1d = member.data.flatten()
1580 # Plot the histogram using plt.hist
1581 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1582 plt.hist(
1583 member_data_1d,
1584 density=True,
1585 stacked=True,
1586 label=f"{mtitle}",
1587 )
1589 # Add a legend
1590 ax.legend(fontsize=16)
1592 # Save the figure to a file
1593 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1594 logging.info("Saved histogram postage stamp plot to %s", filename)
1596 # Close the figure
1597 plt.close(fig)
1600def _plot_and_save_scattermap_plot(
1601 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1602):
1603 """Plot and save a geographical scatter plot.
1605 Parameters
1606 ----------
1607 cube: Cube
1608 1 dimensional Cube of the data points with auxiliary latitude and
1609 longitude coordinates,
1610 filename: str
1611 Filename of the plot to write.
1612 title: str
1613 Plot title.
1614 projection: str
1615 Mapping projection to be used by cartopy.
1616 """
1617 # Setup plot details, size, resolution, etc.
1618 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1619 if projection is not None:
1620 # Apart from the default, the only projection we currently support is
1621 # a stereographic projection over the North Pole.
1622 if projection == "NP_Stereo":
1623 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1624 else:
1625 raise ValueError(f"Unknown projection: {projection}")
1626 else:
1627 axes = plt.axes(projection=ccrs.PlateCarree())
1629 # Scatter plot of the field. The marker size is chosen to give
1630 # symbols that decrease in size as the number of observations
1631 # increases, although the fraction of the figure covered by
1632 # symbols increases roughly as N^(1/2), disregarding overlaps,
1633 # and has been selected for the default figure size of (10, 10).
1634 # Should this be changed, the marker size should be adjusted in
1635 # proportion to the area of the figure.
1636 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1637 klon = None
1638 klat = None
1639 for kc in range(len(cube.aux_coords)):
1640 if cube.aux_coords[kc].standard_name == "latitude":
1641 klat = kc
1642 elif cube.aux_coords[kc].standard_name == "longitude":
1643 klon = kc
1644 scatter_map = iplt.scatter(
1645 cube.aux_coords[klon],
1646 cube.aux_coords[klat],
1647 c=cube.data[:],
1648 s=mrk_size,
1649 cmap="jet",
1650 edgecolors="k",
1651 )
1653 # Add coastlines and borderlines.
1654 try:
1655 axes.coastlines(resolution="10m")
1656 axes.add_feature(cfeature.BORDERS)
1657 except AttributeError:
1658 pass
1660 # Add title.
1661 axes.set_title(title, fontsize=16)
1663 # Add colour bar.
1664 cbar = fig.colorbar(scatter_map)
1665 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1667 # Save plot.
1668 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1669 logging.info("Saved geographical scatter plot to %s", filename)
1670 plt.close(fig)
1673def _spatial_plot(
1674 method: Literal["contourf", "pcolormesh"],
1675 cube: iris.cube.Cube,
1676 filename: str | None,
1677 sequence_coordinate: str,
1678 stamp_coordinate: str,
1679 overlay_cube: iris.cube.Cube | None = None,
1680 contour_cube: iris.cube.Cube | None = None,
1681 **kwargs,
1682):
1683 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1685 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1686 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1687 is present then postage stamp plots will be produced.
1689 If an overlay_cube and/or contour_cube are specified, multiple variables can
1690 be overplotted on the same figure.
1692 Parameters
1693 ----------
1694 method: "contourf" | "pcolormesh"
1695 The plotting method to use.
1696 cube: Cube
1697 Iris cube of the data to plot. It should have two spatial dimensions,
1698 such as lat and lon, and may also have a another two dimension to be
1699 plotted sequentially and/or as postage stamp plots.
1700 filename: str | None
1701 Name of the plot to write, used as a prefix for plot sequences. If None
1702 uses the recipe name.
1703 sequence_coordinate: str
1704 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1705 This coordinate must exist in the cube.
1706 stamp_coordinate: str
1707 Coordinate about which to plot postage stamp plots. Defaults to
1708 ``"realization"``.
1709 overlay_cube: Cube | None, optional
1710 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1711 contour_cube: Cube | None, optional
1712 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1714 Raises
1715 ------
1716 ValueError
1717 If the cube doesn't have the right dimensions.
1718 TypeError
1719 If the cube isn't a single cube.
1720 """
1721 recipe_title = get_recipe_metadata().get("title", "Untitled")
1723 # Ensure we've got a single cube.
1724 cube = check_single_cube(cube)
1726 # Check if there is a valid stamp coordinate in cube dimensions.
1727 if stamp_coordinate == "realization": 1727 ↛ 1732line 1727 didn't jump to line 1732 because the condition on line 1727 was always true
1728 stamp_coordinate = check_stamp_coordinate(cube)
1730 # Make postage stamp plots if stamp_coordinate exists and has more than a
1731 # single point.
1732 plotting_func = _plot_and_save_spatial_plot
1733 try:
1734 if cube.coord(stamp_coordinate).shape[0] > 1:
1735 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1736 except iris.exceptions.CoordinateNotFoundError:
1737 pass
1739 # Produce a geographical scatter plot if the data have a
1740 # dimension called observation or model_obs_error
1741 if any( 1741 ↛ 1745line 1741 didn't jump to line 1745 because the condition on line 1741 was never true
1742 crd.var_name == "station" or crd.var_name == "model_obs_error"
1743 for crd in cube.coords()
1744 ):
1745 plotting_func = _plot_and_save_scattermap_plot
1747 # Must have a sequence coordinate.
1748 try:
1749 cube.coord(sequence_coordinate)
1750 except iris.exceptions.CoordinateNotFoundError as err:
1751 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1753 # Create a plot for each value of the sequence coordinate.
1754 plot_index = []
1755 nplot = np.size(cube.coord(sequence_coordinate).points)
1757 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1758 # Set plot titles and filename
1759 seq_coord = cube_slice.coord(sequence_coordinate)
1760 plot_title, plot_filename = _set_title_and_filename(
1761 seq_coord, nplot, recipe_title, filename
1762 )
1764 # Extract sequence slice for overlay_cube and contour_cube if required.
1765 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1766 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1768 # Do the actual plotting.
1769 plotting_func(
1770 cube_slice,
1771 filename=plot_filename,
1772 stamp_coordinate=stamp_coordinate,
1773 title=plot_title,
1774 method=method,
1775 overlay_cube=overlay_slice,
1776 contour_cube=contour_slice,
1777 **kwargs,
1778 )
1779 plot_index.append(plot_filename)
1781 # Add list of plots to plot metadata.
1782 complete_plot_index = _append_to_plot_index(plot_index)
1784 # Make a page to display the plots.
1785 _make_plot_html_page(complete_plot_index)
1788####################
1789# Public functions #
1790####################
1793def spatial_contour_plot(
1794 cube: iris.cube.Cube,
1795 filename: str = None,
1796 sequence_coordinate: str = "time",
1797 stamp_coordinate: str = "realization",
1798 **kwargs,
1799) -> iris.cube.Cube:
1800 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1802 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1803 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1804 is present then postage stamp plots will be produced.
1806 Parameters
1807 ----------
1808 cube: Cube
1809 Iris cube of the data to plot. It should have two spatial dimensions,
1810 such as lat and lon, and may also have a another two dimension to be
1811 plotted sequentially and/or as postage stamp plots.
1812 filename: str, optional
1813 Name of the plot to write, used as a prefix for plot sequences. Defaults
1814 to the recipe name.
1815 sequence_coordinate: str, optional
1816 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1817 This coordinate must exist in the cube.
1818 stamp_coordinate: str, optional
1819 Coordinate about which to plot postage stamp plots. Defaults to
1820 ``"realization"``.
1822 Returns
1823 -------
1824 Cube
1825 The original cube (so further operations can be applied).
1827 Raises
1828 ------
1829 ValueError
1830 If the cube doesn't have the right dimensions.
1831 TypeError
1832 If the cube isn't a single cube.
1833 """
1834 _spatial_plot(
1835 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1836 )
1837 return cube
1840def spatial_pcolormesh_plot(
1841 cube: iris.cube.Cube,
1842 filename: str = None,
1843 sequence_coordinate: str = "time",
1844 stamp_coordinate: str = "realization",
1845 **kwargs,
1846) -> iris.cube.Cube:
1847 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1849 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1850 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1851 is present then postage stamp plots will be produced.
1853 This function is significantly faster than ``spatial_contour_plot``,
1854 especially at high resolutions, and should be preferred unless contiguous
1855 contour areas are important.
1857 Parameters
1858 ----------
1859 cube: Cube
1860 Iris cube of the data to plot. It should have two spatial dimensions,
1861 such as lat and lon, and may also have a another two dimension to be
1862 plotted sequentially and/or as postage stamp plots.
1863 filename: str, optional
1864 Name of the plot to write, used as a prefix for plot sequences. Defaults
1865 to the recipe name.
1866 sequence_coordinate: str, optional
1867 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1868 This coordinate must exist in the cube.
1869 stamp_coordinate: str, optional
1870 Coordinate about which to plot postage stamp plots. Defaults to
1871 ``"realization"``.
1873 Returns
1874 -------
1875 Cube
1876 The original cube (so further operations can be applied).
1878 Raises
1879 ------
1880 ValueError
1881 If the cube doesn't have the right dimensions.
1882 TypeError
1883 If the cube isn't a single cube.
1884 """
1885 _spatial_plot(
1886 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1887 )
1888 return cube
1891def spatial_multi_pcolormesh_plot(
1892 cube: iris.cube.Cube,
1893 overlay_cube: iris.cube.Cube | None = None,
1894 contour_cube: iris.cube.Cube | None = None,
1895 filename: str = None,
1896 sequence_coordinate: str = "time",
1897 stamp_coordinate: str = "realization",
1898 **kwargs,
1899) -> iris.cube.Cube:
1900 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1902 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1903 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1904 is present then postage stamp plots will be produced.
1906 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1908 If specified, contours of a contour_cube can be overplotted on top of those.
1910 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1912 This function is significantly faster than ``spatial_contour_plot``,
1913 especially at high resolutions, and should be preferred unless contiguous
1914 contour areas are important.
1916 Parameters
1917 ----------
1918 cube: Cube
1919 Iris cube of the data to plot. It should have two spatial dimensions,
1920 such as lat and lon, and may also have a another two dimension to be
1921 plotted sequentially and/or as postage stamp plots.
1922 overlay_cube: Cube, optional
1923 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1924 such as lat and lon, and may also have a another two dimension to be
1925 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.
1926 If not provided, output plot generated without overlay cube.
1927 contour_cube: Cube, optional
1928 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,
1929 such as lat and lon, and may also have a another two dimension to be
1930 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
1931 filename: str, optional
1932 Name of the plot to write, used as a prefix for plot sequences. Defaults
1933 to the recipe name.
1934 sequence_coordinate: str, optional
1935 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1936 This coordinate must exist in the cube.
1937 stamp_coordinate: str, optional
1938 Coordinate about which to plot postage stamp plots. Defaults to
1939 ``"realization"``.
1941 Returns
1942 -------
1943 Cube
1944 The original cube (so further operations can be applied).
1946 Raises
1947 ------
1948 ValueError
1949 If the cube doesn't have the right dimensions.
1950 TypeError
1951 If the cube isn't a single cube.
1952 """
1953 _spatial_plot(
1954 "pcolormesh",
1955 cube,
1956 filename,
1957 sequence_coordinate,
1958 stamp_coordinate,
1959 overlay_cube=overlay_cube,
1960 contour_cube=contour_cube,
1961 )
1962 return cube, overlay_cube, contour_cube
1965# TODO: Expand function to handle ensemble data.
1966# line_coordinate: str, optional
1967# Coordinate about which to plot multiple lines. Defaults to
1968# ``"realization"``.
1969def plot_line_series(
1970 cube: iris.cube.Cube | iris.cube.CubeList,
1971 filename: str = None,
1972 series_coordinate: str = "time",
1973 sequence_coordinate: str = "time",
1974 # add the following for ensembles
1975 stamp_coordinate: str = "realization",
1976 single_plot: bool = False,
1977 **kwargs,
1978) -> iris.cube.Cube | iris.cube.CubeList:
1979 """Plot a line plot for the specified coordinate.
1981 The Cube or CubeList must be 1D.
1983 Parameters
1984 ----------
1985 iris.cube | iris.cube.CubeList
1986 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
1987 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
1988 We do not support different data such as temperature and humidity in the same CubeList for plotting.
1989 filename: str, optional
1990 Name of the plot to write, used as a prefix for plot sequences. Defaults
1991 to the recipe name.
1992 series_coordinate: str, optional
1993 Coordinate about which to make a series. Defaults to ``"time"``. This
1994 coordinate must exist in the cube.
1996 Returns
1997 -------
1998 iris.cube.Cube | iris.cube.CubeList
1999 The original Cube or CubeList (so further operations can be applied).
2000 plotted data.
2002 Raises
2003 ------
2004 ValueError
2005 If the cubes don't have the right dimensions.
2006 TypeError
2007 If the cube isn't a Cube or CubeList.
2008 """
2009 # Ensure we have a name for the plot file.
2010 recipe_title = get_recipe_metadata().get("title", "Untitled")
2012 num_models = get_num_models(cube)
2014 validate_cube_shape(cube, num_models)
2016 # Iterate over all cubes and extract coordinate to plot.
2017 cubes = iter_maybe(cube)
2018 coords = []
2019 for cube in cubes:
2020 try:
2021 coords.append(cube.coord(series_coordinate))
2022 except iris.exceptions.CoordinateNotFoundError as err:
2023 raise ValueError(
2024 f"Cube must have a {series_coordinate} coordinate."
2025 ) from err
2026 if cube.coords("realization"): 2026 ↛ 2030line 2026 didn't jump to line 2030 because the condition on line 2026 was always true
2027 if cube.ndim > 3: 2027 ↛ 2028line 2027 didn't jump to line 2028 because the condition on line 2027 was never true
2028 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2029 else:
2030 raise ValueError("Cube must have a realization coordinate.")
2032 plot_index = []
2034 # Check if this is a spectral plot by looking for spectral coordinates
2035 is_spectral_plot = series_coordinate in [
2036 "frequency",
2037 "physical_wavenumber",
2038 "wavelength",
2039 ]
2041 if is_spectral_plot:
2042 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2043 # coordinate frequency/wavenumber.
2044 # If several power spectra are plotted with time as sequence_coordinate for the
2045 # time slider option.
2047 # Internal plotting function.
2048 plotting_func = _plot_and_save_line_power_spectrum_series
2050 for cube in cubes:
2051 try:
2052 cube.coord(sequence_coordinate)
2053 except iris.exceptions.CoordinateNotFoundError as err:
2054 raise ValueError(
2055 f"Cube must have a {sequence_coordinate} coordinate."
2056 ) from err
2058 if num_models == 1: 2058 ↛ 2072line 2058 didn't jump to line 2072 because the condition on line 2058 was always true
2059 # check for ensembles
2060 if ( 2060 ↛ 2064line 2060 didn't jump to line 2064 because the condition on line 2060 was never true
2061 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2062 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2063 ):
2064 if single_plot:
2065 # Plot spectra, mean and ensemble spread on 1 plot
2066 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2067 else:
2068 # Plot postage stamps
2069 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2070 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2071 else:
2072 all_points = sorted(
2073 set(
2074 itertools.chain.from_iterable(
2075 cb.coord(sequence_coordinate).points for cb in cubes
2076 )
2077 )
2078 )
2079 all_slices = list(
2080 itertools.chain.from_iterable(
2081 cb.slices_over(sequence_coordinate) for cb in cubes
2082 )
2083 )
2084 # Matched slices (matched by seq coord point; it may happen that
2085 # evaluated models do not cover the same seq coord range, hence matching
2086 # necessary)
2087 cube_iterables = [
2088 iris.cube.CubeList(
2089 s
2090 for s in all_slices
2091 if s.coord(sequence_coordinate).points[0] == point
2092 )
2093 for point in all_points
2094 ]
2096 nplot = np.size(cube.coord(sequence_coordinate).points)
2098 # Create a plot for each value of the sequence coordinate. Allowing for
2099 # multiple cubes in a CubeList to be plotted in the same plot for similar
2100 # sequence values. Passing a CubeList into the internal plotting function
2101 # for similar values of the sequence coordinate. cube_slice can be an
2102 # iris.cube.Cube or an iris.cube.CubeList.
2104 for cube_slice in cube_iterables:
2105 # Normalize cube_slice to a list of cubes
2106 if isinstance(cube_slice, iris.cube.CubeList): 2106 ↛ 2107line 2106 didn't jump to line 2107 because the condition on line 2106 was never true
2107 cubes = list(cube_slice)
2108 elif isinstance(cube_slice, iris.cube.Cube): 2108 ↛ 2111line 2108 didn't jump to line 2111 because the condition on line 2108 was always true
2109 cubes = [cube_slice]
2110 else:
2111 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2113 # Use sequence value so multiple sequences can merge.
2114 seq_coord = cube_slice[0].coord(sequence_coordinate)
2115 plot_title, plot_filename = _set_title_and_filename(
2116 seq_coord, nplot, recipe_title, filename
2117 )
2119 # Format the coordinate value in a unit appropriate way.
2120 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2122 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2123 if nplot == 1 and seq_coord.has_bounds: 2123 ↛ 2128line 2123 didn't jump to line 2128 because the condition on line 2123 was always true
2124 if np.size(seq_coord.bounds) > 1: 2124 ↛ 2125line 2124 didn't jump to line 2125 because the condition on line 2124 was never true
2125 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2127 # Do the actual plotting.
2128 plotting_func(
2129 cube_slice,
2130 coords,
2131 stamp_coordinate,
2132 plot_filename,
2133 title,
2134 series_coordinate,
2135 )
2137 plot_index.append(plot_filename)
2138 else:
2139 # Format the title and filename using plotted series coordinate
2140 nplot = 1
2141 seq_coord = coords[0]
2142 plot_title, plot_filename = _set_title_and_filename(
2143 seq_coord, nplot, recipe_title, filename
2144 )
2145 # Do the actual plotting for all other series coordinate options.
2146 _plot_and_save_line_series(
2147 cubes, coords, stamp_coordinate, plot_filename, plot_title
2148 )
2150 plot_index.append(plot_filename)
2152 # append plot to list of plots
2153 complete_plot_index = _append_to_plot_index(plot_index)
2155 # Make a page to display the plots.
2156 _make_plot_html_page(complete_plot_index)
2158 return cube
2161def plot_vertical_line_series(
2162 cubes: iris.cube.Cube | iris.cube.CubeList,
2163 filename: str = None,
2164 series_coordinate: str = "model_level_number",
2165 sequence_coordinate: str = "time",
2166 # line_coordinate: str = "realization",
2167 **kwargs,
2168) -> iris.cube.Cube | iris.cube.CubeList:
2169 """Plot a line plot against a type of vertical coordinate.
2171 The Cube or CubeList must be 1D.
2173 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2174 then a sequence of plots will be produced.
2176 Parameters
2177 ----------
2178 iris.cube | iris.cube.CubeList
2179 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2180 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2181 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2182 filename: str, optional
2183 Name of the plot to write, used as a prefix for plot sequences. Defaults
2184 to the recipe name.
2185 series_coordinate: str, optional
2186 Coordinate to plot on the y-axis. Can be ``pressure`` or
2187 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2188 for LFRic. Defaults to ``model_level_number``.
2189 This coordinate must exist in the cube.
2190 sequence_coordinate: str, optional
2191 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2192 This coordinate must exist in the cube.
2194 Returns
2195 -------
2196 iris.cube.Cube | iris.cube.CubeList
2197 The original Cube or CubeList (so further operations can be applied).
2198 Plotted data.
2200 Raises
2201 ------
2202 ValueError
2203 If the cubes doesn't have the right dimensions.
2204 TypeError
2205 If the cube isn't a Cube or CubeList.
2206 """
2207 # Ensure we have a name for the plot file.
2208 recipe_title = get_recipe_metadata().get("title", "Untitled")
2210 cubes = iter_maybe(cubes)
2211 # Initialise empty list to hold all data from all cubes in a CubeList
2212 all_data = []
2214 # Store min/max ranges for x range.
2215 x_levels = []
2217 num_models = get_num_models(cubes)
2219 validate_cube_shape(cubes, num_models)
2221 # Iterate over all cubes in cube or CubeList and plot.
2222 coords = []
2223 for cube in cubes:
2224 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2225 try:
2226 coords.append(cube.coord(series_coordinate))
2227 except iris.exceptions.CoordinateNotFoundError as err:
2228 raise ValueError(
2229 f"Cube must have a {series_coordinate} coordinate."
2230 ) from err
2232 try:
2233 if cube.ndim > 1 or not cube.coords("realization"): 2233 ↛ 2241line 2233 didn't jump to line 2241 because the condition on line 2233 was always true
2234 cube.coord(sequence_coordinate)
2235 except iris.exceptions.CoordinateNotFoundError as err:
2236 raise ValueError(
2237 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2238 ) from err
2240 # Get minimum and maximum from levels information.
2241 _, levels, _ = colorbar_map_levels(cube, axis="x")
2242 if levels is not None: 2242 ↛ 2246line 2242 didn't jump to line 2246 because the condition on line 2242 was always true
2243 x_levels.append(min(levels))
2244 x_levels.append(max(levels))
2245 else:
2246 all_data.append(cube.data)
2248 if len(x_levels) == 0: 2248 ↛ 2250line 2248 didn't jump to line 2250 because the condition on line 2248 was never true
2249 # Combine all data into a single NumPy array
2250 combined_data = np.concatenate(all_data)
2252 # Set the lower and upper limit for the x-axis to ensure all plots have
2253 # same range. This needs to read the whole cube over the range of the
2254 # sequence and if applicable postage stamp coordinate.
2255 vmin = np.floor(combined_data.min())
2256 vmax = np.ceil(combined_data.max())
2257 else:
2258 vmin = min(x_levels)
2259 vmax = max(x_levels)
2261 # Matching the slices (matching by seq coord point; it may happen that
2262 # evaluated models do not cover the same seq coord range, hence matching
2263 # necessary)
2264 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2266 # Create a plot for each value of the sequence coordinate.
2267 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2268 # similar sequence values. Passing a CubeList into the internal plotting function
2269 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2270 # or an iris.cube.CubeList.
2271 plot_index = []
2272 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2273 for cubes_slice in cube_iterables:
2274 # Format the coordinate value in a unit appropriate way.
2275 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2276 plot_title, plot_filename = _set_title_and_filename(
2277 seq_coord, nplot, recipe_title, filename
2278 )
2280 # Do the actual plotting.
2281 _plot_and_save_vertical_line_series(
2282 cubes_slice,
2283 coords,
2284 "realization",
2285 plot_filename,
2286 series_coordinate,
2287 title=plot_title,
2288 vmin=vmin,
2289 vmax=vmax,
2290 )
2291 plot_index.append(plot_filename)
2293 # Add list of plots to plot metadata.
2294 complete_plot_index = _append_to_plot_index(plot_index)
2296 # Make a page to display the plots.
2297 _make_plot_html_page(complete_plot_index)
2299 return cubes
2302def qq_plot(
2303 cubes: iris.cube.CubeList,
2304 coordinates: list[str],
2305 percentiles: list[float],
2306 model_names: list[str],
2307 filename: str = None,
2308 one_to_one: bool = True,
2309 **kwargs,
2310) -> iris.cube.CubeList:
2311 """Plot a Quantile-Quantile plot between two models for common time points.
2313 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2314 collapsed within the operator over all specified coordinates such as
2315 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2316 ensemble members to ensure a 1D cube (array).
2318 Parameters
2319 ----------
2320 cubes: iris.cube.CubeList
2321 Two cubes of the same variable with different models.
2322 coordinate: list[str]
2323 The list of coordinates to collapse over. This list should be
2324 every coordinate within the cube to result in a 1D cube around
2325 the percentile coordinate.
2326 percent: list[float]
2327 A list of percentiles to appear in the plot.
2328 model_names: list[str]
2329 A list of model names to appear on the axis of the plot.
2330 filename: str, optional
2331 Filename of the plot to write.
2332 one_to_one: bool, optional
2333 If True a 1:1 line is plotted; if False it is not. Default is True.
2335 Raises
2336 ------
2337 ValueError
2338 When the cubes are not compatible.
2340 Notes
2341 -----
2342 The quantile-quantile plot is a variant on the scatter plot representing
2343 two datasets by their quantiles (percentiles) for common time points.
2344 This plot does not use a theoretical distribution to compare against, but
2345 compares percentiles of two datasets. This plot does
2346 not use all raw data points, but plots the selected percentiles (quantiles) of
2347 each variable instead for the two datasets, thereby normalising the data for a
2348 direct comparison between the selected percentiles of the two dataset distributions.
2350 Quantile-quantile plots are valuable for comparing against
2351 observations and other models. Identical percentiles between the variables
2352 will lie on the one-to-one line implying the values correspond well to each
2353 other. Where there is a deviation from the one-to-one line a range of
2354 possibilities exist depending on how and where the data is shifted (e.g.,
2355 Wilks 2011 [Wilks2011]_).
2357 For distributions above the one-to-one line the distribution is left-skewed;
2358 below is right-skewed. A distinct break implies a bimodal distribution, and
2359 closer values/values further apart at the tails imply poor representation of
2360 the extremes.
2362 References
2363 ----------
2364 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2365 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2366 """
2367 # Check cubes using same functionality as the difference operator.
2368 if len(cubes) != 2:
2369 raise ValueError("cubes should contain exactly 2 cubes.")
2370 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2371 other: Cube = cubes.extract_cube(
2372 iris.Constraint(
2373 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2374 )
2375 )
2377 # Get spatial coord names.
2378 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2379 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2381 # Ensure cubes to compare are on common differencing grid.
2382 # This is triggered if either
2383 # i) latitude and longitude shapes are not the same. Note grid points
2384 # are not compared directly as these can differ through rounding
2385 # errors.
2386 # ii) or variables are known to often sit on different grid staggering
2387 # in different models (e.g. cell center vs cell edge), as is the case
2388 # for UM and LFRic comparisons.
2389 # In future greater choice of regridding method might be applied depending
2390 # on variable type. Linear regridding can in general be appropriate for smooth
2391 # variables. Care should be taken with interpretation of differences
2392 # given this dependency on regridding.
2393 if (
2394 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2395 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2396 ) or (
2397 base.long_name
2398 in [
2399 "eastward_wind_at_10m",
2400 "northward_wind_at_10m",
2401 "northward_wind_at_cell_centres",
2402 "eastward_wind_at_cell_centres",
2403 "zonal_wind_at_pressure_levels",
2404 "meridional_wind_at_pressure_levels",
2405 "potential_vorticity_at_pressure_levels",
2406 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2407 ]
2408 ):
2409 logging.debug(
2410 "Linear regridding base cube to other grid to compute differences"
2411 )
2412 base = regrid_onto_cube(base, other, method="Linear")
2414 # Extract just common time points.
2415 base, other = _extract_common_time_points(base, other)
2417 # Equalise attributes so we can merge.
2418 fully_equalise_attributes([base, other])
2419 logging.debug("Base: %s\nOther: %s", base, other)
2421 # Collapse cubes.
2422 base = collapse(
2423 base,
2424 coordinate=coordinates,
2425 method="PERCENTILE",
2426 additional_percent=percentiles,
2427 )
2428 other = collapse(
2429 other,
2430 coordinate=coordinates,
2431 method="PERCENTILE",
2432 additional_percent=percentiles,
2433 )
2435 # Ensure we have a name for the plot file.
2436 recipe_title = get_recipe_metadata().get("title", "Untitled")
2437 title = f"{recipe_title}"
2439 if filename is None:
2440 filename = slugify(recipe_title)
2442 # Add file extension.
2443 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2445 # Do the actual plotting on a scatter plot
2446 _plot_and_save_scatter_plot(
2447 base, other, plot_filename, title, one_to_one, model_names
2448 )
2450 # Add list of plots to plot metadata.
2451 plot_index = _append_to_plot_index([plot_filename])
2453 # Make a page to display the plots.
2454 _make_plot_html_page(plot_index)
2456 return iris.cube.CubeList([base, other])
2459def scatter_plot(
2460 cube_x: iris.cube.Cube | iris.cube.CubeList,
2461 cube_y: iris.cube.Cube | iris.cube.CubeList,
2462 filename: str = None,
2463 one_to_one: bool = True,
2464 **kwargs,
2465) -> iris.cube.CubeList:
2466 """Plot a scatter plot between two variables.
2468 Both cubes must be 1D.
2470 Parameters
2471 ----------
2472 cube_x: Cube | CubeList
2473 1 dimensional Cube of the data to plot on y-axis.
2474 cube_y: Cube | CubeList
2475 1 dimensional Cube of the data to plot on x-axis.
2476 filename: str, optional
2477 Filename of the plot to write.
2478 one_to_one: bool, optional
2479 If True a 1:1 line is plotted; if False it is not. Default is True.
2481 Returns
2482 -------
2483 cubes: CubeList
2484 CubeList of the original x and y cubes for further processing.
2486 Raises
2487 ------
2488 ValueError
2489 If the cube doesn't have the right dimensions and cubes not the same
2490 size.
2491 TypeError
2492 If the cube isn't a single cube.
2494 Notes
2495 -----
2496 Scatter plots are used for determining if there is a relationship between
2497 two variables. Positive relations have a slope going from bottom left to top
2498 right; Negative relations have a slope going from top left to bottom right.
2499 """
2500 # Iterate over all cubes in cube or CubeList and plot.
2501 for cube_iter in iter_maybe(cube_x):
2502 # Check cubes are correct shape.
2503 cube_iter = check_single_cube(cube_iter)
2504 if cube_iter.ndim > 1:
2505 raise ValueError("cube_x must be 1D.")
2507 # Iterate over all cubes in cube or CubeList and plot.
2508 for cube_iter in iter_maybe(cube_y):
2509 # Check cubes are correct shape.
2510 cube_iter = check_single_cube(cube_iter)
2511 if cube_iter.ndim > 1:
2512 raise ValueError("cube_y must be 1D.")
2514 # Ensure we have a name for the plot file.
2515 recipe_title = get_recipe_metadata().get("title", "Untitled")
2516 title = f"{recipe_title}"
2518 if filename is None:
2519 filename = slugify(recipe_title)
2521 # Add file extension.
2522 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2524 # Do the actual plotting.
2525 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2527 # Add list of plots to plot metadata.
2528 plot_index = _append_to_plot_index([plot_filename])
2530 # Make a page to display the plots.
2531 _make_plot_html_page(plot_index)
2533 return iris.cube.CubeList([cube_x, cube_y])
2536def vector_plot(
2537 cube_u: iris.cube.Cube,
2538 cube_v: iris.cube.Cube,
2539 filename: str = None,
2540 sequence_coordinate: str = "time",
2541 **kwargs,
2542) -> iris.cube.CubeList:
2543 """Plot a vector plot based on the input u and v components."""
2544 recipe_title = get_recipe_metadata().get("title", "Untitled")
2546 # Cubes must have a matching sequence coordinate.
2547 try:
2548 # Check that the u and v cubes have the same sequence coordinate.
2549 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2549 ↛ anywhereline 2549 didn't jump anywhere: it always raised an exception.
2550 raise ValueError("Coordinates do not match.")
2551 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2552 raise ValueError(
2553 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2554 ) from err
2556 # Create a plot for each value of the sequence coordinate.
2557 plot_index = []
2558 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2559 for cube_u_slice, cube_v_slice in zip(
2560 cube_u.slices_over(sequence_coordinate),
2561 cube_v.slices_over(sequence_coordinate),
2562 strict=True,
2563 ):
2564 # Format the coordinate value in a unit appropriate way.
2565 seq_coord = cube_u_slice.coord(sequence_coordinate)
2566 plot_title, plot_filename = _set_title_and_filename(
2567 seq_coord, nplot, recipe_title, filename
2568 )
2570 # Do the actual plotting.
2571 _plot_and_save_vector_plot(
2572 cube_u_slice,
2573 cube_v_slice,
2574 filename=plot_filename,
2575 title=plot_title,
2576 method="pcolormesh",
2577 )
2578 plot_index.append(plot_filename)
2580 # Add list of plots to plot metadata.
2581 complete_plot_index = _append_to_plot_index(plot_index)
2583 # Make a page to display the plots.
2584 _make_plot_html_page(complete_plot_index)
2586 return iris.cube.CubeList([cube_u, cube_v])
2589def plot_histogram_series(
2590 cubes: iris.cube.Cube | iris.cube.CubeList,
2591 filename: str = None,
2592 sequence_coordinate: str = "time",
2593 stamp_coordinate: str = "realization",
2594 single_plot: bool = False,
2595 **kwargs,
2596) -> iris.cube.Cube | iris.cube.CubeList:
2597 """Plot a histogram plot for each vertical level provided.
2599 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2600 is present then a sequence of plots will be produced using the time slider
2601 functionality to scroll through histograms against time. If a
2602 stamp_coordinate is present then postage stamp plots will be produced. If
2603 stamp_coordinate and single_plot is True, all postage stamp plots will be
2604 plotted in a single plot instead of separate postage stamp plots.
2606 Parameters
2607 ----------
2608 cubes: Cube | iris.cube.CubeList
2609 Iris cube or CubeList of the data to plot. It should have a single dimension other
2610 than the stamp coordinate.
2611 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2612 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2613 filename: str, optional
2614 Name of the plot to write, used as a prefix for plot sequences. Defaults
2615 to the recipe name.
2616 sequence_coordinate: str, optional
2617 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2618 This coordinate must exist in the cube and will be used for the time
2619 slider.
2620 stamp_coordinate: str, optional
2621 Coordinate about which to plot postage stamp plots. Defaults to
2622 ``"realization"``.
2623 single_plot: bool, optional
2624 If True, all postage stamp plots will be plotted in a single plot. If
2625 False, each postage stamp plot will be plotted separately. Is only valid
2626 if stamp_coordinate exists and has more than a single point.
2628 Returns
2629 -------
2630 iris.cube.Cube | iris.cube.CubeList
2631 The original Cube or CubeList (so further operations can be applied).
2632 Plotted data.
2634 Raises
2635 ------
2636 ValueError
2637 If the cube doesn't have the right dimensions.
2638 TypeError
2639 If the cube isn't a Cube or CubeList.
2640 """
2641 recipe_title = get_recipe_metadata().get("title", "Untitled")
2643 cubes = iter_maybe(cubes)
2644 # Ensure we have a name for the plot file.
2645 if filename is None:
2646 filename = slugify(recipe_title)
2648 # Internal plotting function.
2649 plotting_func = _plot_and_save_histogram_series
2651 num_models = get_num_models(cubes)
2653 validate_cube_shape(cubes, num_models)
2655 # If several histograms are plotted, check sequence_coordinate
2656 check_sequence_coordinate(cubes, sequence_coordinate)
2658 # Get axis minimum and maximum from levels information.
2659 # If no levels set, derive minima and maxima from data in CubeList.
2660 vmin, vmax = _set_axis_range(cubes)
2662 # Make postage stamp plots if stamp_coordinate exists and has more than a
2663 # single point. If single_plot is True:
2664 # -- all postage stamp plots will be plotted in a single plot instead of
2665 # separate postage stamp plots.
2666 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2667 # produced per single model only
2668 if num_models == 1: 2668 ↛ 2681line 2668 didn't jump to line 2681 because the condition on line 2668 was always true
2669 if ( 2669 ↛ 2673line 2669 didn't jump to line 2673 because the condition on line 2669 was never true
2670 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2671 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2672 ):
2673 if single_plot:
2674 plotting_func = (
2675 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2676 )
2677 else:
2678 plotting_func = _plot_and_save_postage_stamp_histogram_series
2679 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2680 else:
2681 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2683 plot_index = []
2684 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2685 # Create a plot for each value of the sequence coordinate. Allowing for
2686 # multiple cubes in a CubeList to be plotted in the same plot for similar
2687 # sequence values. Passing a CubeList into the internal plotting function
2688 # for similar values of the sequence coordinate. cube_slice can be an
2689 # iris.cube.Cube or an iris.cube.CubeList.
2690 for cube_slice in cube_iterables:
2691 single_cube = cube_slice
2692 if isinstance(cube_slice, iris.cube.CubeList): 2692 ↛ 2693line 2692 didn't jump to line 2693 because the condition on line 2692 was never true
2693 single_cube = cube_slice[0]
2695 # Ensure valid stamp coordinate in cube dimensions
2696 if stamp_coordinate == "realization": 2696 ↛ 2699line 2696 didn't jump to line 2699 because the condition on line 2696 was always true
2697 stamp_coordinate = check_stamp_coordinate(single_cube)
2698 # Set plot titles and filename, based on sequence coordinate
2699 seq_coord = single_cube.coord(sequence_coordinate)
2700 # Use time coordinate in title and filename if single histogram output.
2701 if sequence_coordinate == "realization" and nplot == 1: 2701 ↛ 2702line 2701 didn't jump to line 2702 because the condition on line 2701 was never true
2702 seq_coord = single_cube.coord("time")
2703 plot_title, plot_filename = _set_title_and_filename(
2704 seq_coord, nplot, recipe_title, filename
2705 )
2707 # Do the actual plotting.
2708 plotting_func(
2709 cube_slice,
2710 filename=plot_filename,
2711 stamp_coordinate=stamp_coordinate,
2712 title=plot_title,
2713 vmin=vmin,
2714 vmax=vmax,
2715 )
2716 plot_index.append(plot_filename)
2718 # Add list of plots to plot metadata.
2719 complete_plot_index = _append_to_plot_index(plot_index)
2721 # Make a page to display the plots.
2722 _make_plot_html_page(complete_plot_index)
2724 return cubes
2727def _plot_and_save_postage_stamp_power_spectrum_series(
2728 cubes: iris.cube.Cube,
2729 coords: list[iris.coords.Coord],
2730 stamp_coordinate: str,
2731 filename: str,
2732 title: str,
2733 series_coordinate: str = None,
2734 **kwargs,
2735):
2736 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2738 Parameters
2739 ----------
2740 cubes: Cube or CubeList
2741 Cube or Cubelist of the power spectrum data.
2742 coords: list[Coord]
2743 Coordinates to plot on the x-axis, one per cube.
2744 stamp_coordinate: str
2745 Coordinate that becomes different plots.
2746 filename: str
2747 Filename of the plot to write.
2748 title: str
2749 Plot title.
2750 series_coordinate: str, optional
2751 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2753 """
2754 # Use the smallest square grid that will fit the members.
2755 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
2757 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
2758 model_colors_map = get_model_colors_map(cubes)
2759 # ax = plt.gca()
2760 # Make a subplot for each member.
2761 for member, subplot in zip(
2762 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
2763 ):
2764 ax = plt.subplot(grid_size, grid_size, subplot)
2766 # Store min/max ranges.
2767 y_levels = []
2769 line_marker = None
2770 line_width = 1
2772 for cube in iter_maybe(member):
2773 xcoord = _select_series_coord(cube, series_coordinate)
2774 xname = xcoord.points
2776 yfield = cube.data # power spectrum
2777 label = None
2778 color = "black"
2779 if model_colors_map: 2779 ↛ 2780line 2779 didn't jump to line 2780 because the condition on line 2779 was never true
2780 label = cube.attributes.get("model_name")
2781 color = model_colors_map.get(label)
2783 if member.coord(stamp_coordinate).points == [0]:
2784 ax.plot(
2785 xname,
2786 yfield,
2787 color=color,
2788 marker=line_marker,
2789 ls="-",
2790 lw=line_width,
2791 label=f"{label} (control)"
2792 if len(cube.coord(stamp_coordinate).points) > 1
2793 else label,
2794 )
2795 # Label with member if part of an ensemble and not the control.
2796 else:
2797 ax.plot(
2798 xname,
2799 yfield,
2800 color=color,
2801 ls="-",
2802 lw=1.5,
2803 alpha=0.75,
2804 label=f"{label} (member)",
2805 )
2807 # Calculate the global min/max if multiple cubes are given.
2808 _, levels, _ = colorbar_map_levels(cube, axis="y")
2809 if levels is not None: 2809 ↛ 2810line 2809 didn't jump to line 2810 because the condition on line 2809 was never true
2810 y_levels.append(min(levels))
2811 y_levels.append(max(levels))
2813 # Add some labels and tweak the style.
2814 title = f"{title}"
2815 ax.set_title(title, fontsize=16)
2817 # Set appropriate x-axis label based on coordinate
2818 if series_coordinate == "wavelength" or ( 2818 ↛ 2821line 2818 didn't jump to line 2821 because the condition on line 2818 was never true
2819 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
2820 ):
2821 ax.set_xlabel("Wavelength (km)", fontsize=14)
2822 elif series_coordinate == "physical_wavenumber" or ( 2822 ↛ 2827line 2822 didn't jump to line 2827 because the condition on line 2822 was always true
2823 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
2824 ):
2825 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2826 else: # frequency or check units
2827 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
2828 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2829 else:
2830 ax.set_xlabel("Wavenumber", fontsize=14)
2832 ax.set_ylabel("Power Spectral Density", fontsize=14)
2833 ax.tick_params(axis="both", labelsize=12)
2835 # Set log-log scale
2836 ax.set_xscale("log")
2837 ax.set_yscale("log")
2839 # Add gridlines
2840 ax.grid(linestyle="--", color="grey", linewidth=1)
2841 # Ientify unique labels for legend
2842 handles = list(
2843 {
2844 label: handle
2845 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
2846 }.values()
2847 )
2848 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
2850 ax = plt.gca()
2851 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
2853 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
2854 logging.info("Saved histogram postage stamp plot to %s", filename)
2855 plt.close(fig)
2858def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
2859 cubes: iris.cube.Cube,
2860 coords: list[iris.coords.Coord],
2861 stamp_coordinate: str,
2862 filename: str,
2863 title: str,
2864 series_coordinate: str = None,
2865 **kwargs,
2866):
2867 """Plot and save power spectra for ensemble members in single plot.
2869 Parameters
2870 ----------
2871 cubes: Cube or CubeList
2872 Cube or Cubelist of the power spectrum data.
2873 coords: list[Coord]
2874 Coordinates to plot on the x-axis, one per cube.
2875 stamp_coordinate: str
2876 Coordinate that becomes different plots.
2877 filename: str
2878 Filename of the plot to write.
2879 title: str
2880 Plot title.
2881 series_coordinate: str, optional
2882 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2884 """
2885 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
2886 model_colors_map = get_model_colors_map(cubes)
2888 line_marker = None
2889 line_width = 1
2891 # Compute ensemble statistics to show spread
2892 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
2893 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
2894 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
2896 xcoord_global = mean_cube.coord(series_coordinate)
2897 x_global = xcoord_global.points
2899 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
2900 xcoord = _select_series_coord(member, series_coordinate)
2901 xname = xcoord.points
2903 yfield = member.data # power spectrum
2904 color = "black"
2905 if model_colors_map: 2905 ↛ 2909line 2905 didn't jump to line 2909 because the condition on line 2905 was always true
2906 label = member.attributes.get("model_name") if i == 0 else None
2907 color = model_colors_map.get(label)
2909 if member.coord(stamp_coordinate).points == [0]:
2910 ax.plot(
2911 xname,
2912 yfield,
2913 color=color,
2914 marker=line_marker,
2915 ls="-",
2916 lw=line_width,
2917 label=f"{label} (control)"
2918 if len(member.coord(stamp_coordinate).points) > 1
2919 else label,
2920 )
2921 # Label with member number if part of an ensemble and not the control.
2922 else:
2923 ax.plot(
2924 xname,
2925 yfield,
2926 color=color,
2927 ls="-",
2928 lw=1.5,
2929 alpha=0.75,
2930 label=label,
2931 )
2933 # Set appropriate x-axis label based on coordinate
2934 if series_coordinate == "wavelength" or ( 2934 ↛ 2937line 2934 didn't jump to line 2937 because the condition on line 2934 was never true
2935 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
2936 ):
2937 ax.set_xlabel("Wavelength (km)", fontsize=14)
2938 elif series_coordinate == "physical_wavenumber" or ( 2938 ↛ 2943line 2938 didn't jump to line 2943 because the condition on line 2938 was always true
2939 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
2940 ):
2941 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2942 else: # frequency or check units
2943 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
2944 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
2945 else:
2946 ax.set_xlabel("Wavenumber", fontsize=14)
2948 # Add ensemble spread shading
2949 ax.fill_between(
2950 x_global,
2951 min_cube.data,
2952 max_cube.data,
2953 color="grey",
2954 alpha=0.3,
2955 label="Ensemble spread",
2956 )
2958 # Add ensemble mean line
2959 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
2961 ax.set_ylabel("Power Spectral Density", fontsize=14)
2962 ax.tick_params(axis="both", labelsize=12)
2964 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
2965 # Set log-log scale
2966 ax.set_xscale("log")
2967 ax.set_yscale("log")
2969 # Add gridlines
2970 ax.grid(linestyle="--", color="grey", linewidth=1)
2971 # Identify unique labels for legend
2972 handles = list(
2973 {
2974 label: handle
2975 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
2976 }.values()
2977 )
2978 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
2980 # Add a legend
2981 ax.legend(fontsize=16)
2983 # Figure title.
2984 ax.set_title(title, fontsize=16)
2986 # Save the figure to a file
2987 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
2989 # Close the figure
2990 plt.close(fig)