Coverage for src/CSET/operators/plot.py: 78%
927 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-18 11:14 +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 | None = None,
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, optional
964 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
965 """
966 # xn = coords[0].name() # x-axis (e.g. frequency)
968 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
969 model_colors_map = get_model_colors_map(cubes)
970 ax = plt.gca()
972 # Store min/max ranges.
973 y_levels = []
975 line_marker = None
976 line_width = 1
978 # for cube, coord in zip(cubes, coords, strict=True):
979 for cube in iter_maybe(cubes):
980 # next 2 lines replace chunk of code.
981 xcoord = _select_series_coord(cube, series_coordinate)
982 xname = xcoord.points
984 yfield = cube.data # power spectrum
985 label = None
986 color = "black"
987 if model_colors_map: 987 ↛ 990line 987 didn't jump to line 990 because the condition on line 987 was always true
988 label = cube.attributes.get("model_name")
989 color = model_colors_map.get(label)
990 for cube_slice in cube.slices_over(ensemble_coord):
991 # Label with (control) if part of an ensemble or not otherwise.
992 if cube_slice.coord(ensemble_coord).points == [0]: 992 ↛ 1006line 992 didn't jump to line 1006 because the condition on line 992 was always true
993 ax.plot(
994 xname,
995 yfield,
996 color=color,
997 marker=line_marker,
998 ls="-",
999 lw=line_width,
1000 label=f"{label} (control)"
1001 if len(cube.coord(ensemble_coord).points) > 1
1002 else label,
1003 )
1004 # Label with (perturbed) if part of an ensemble and not the control.
1005 else:
1006 ax.plot(
1007 xname,
1008 yfield,
1009 color=color,
1010 ls="-",
1011 lw=1.5,
1012 alpha=0.75,
1013 label=f"{label} (member)",
1014 )
1016 # Calculate the global min/max if multiple cubes are given.
1017 _, levels, _ = colorbar_map_levels(cube, axis="y")
1018 if levels is not None: 1018 ↛ 1019line 1018 didn't jump to line 1019 because the condition on line 1018 was never true
1019 y_levels.append(min(levels))
1020 y_levels.append(max(levels))
1022 # Add some labels and tweak the style.
1023 # check if cubes[0] works for single cube if not CubeList
1025 title = f"{title}"
1026 ax.set_title(title, fontsize=16)
1028 # Set appropriate x-axis label based on coordinate
1029 if series_coordinate == "wavelength" 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 == "wavelength"
1031 ):
1032 ax.set_xlabel("Wavelength (km)", fontsize=14)
1033 elif series_coordinate == "physical_wavenumber" or ( 1033 ↛ 1036line 1033 didn't jump to line 1036 because the condition on line 1033 was never true
1034 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1035 ):
1036 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1037 else: # frequency or check units
1038 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1038 ↛ 1039line 1038 didn't jump to line 1039 because the condition on line 1038 was never true
1039 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1040 else:
1041 ax.set_xlabel("Wavenumber", fontsize=14)
1043 ax.set_ylabel("Power Spectral Density", fontsize=14)
1044 ax.tick_params(axis="both", labelsize=12)
1046 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1048 # Set log-log scale
1049 ax.set_xscale("log")
1050 ax.set_yscale("log")
1052 # Add gridlines
1053 ax.grid(linestyle="--", color="grey", linewidth=1)
1054 # Ientify unique labels for legend
1055 handles = list(
1056 {
1057 label: handle
1058 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1059 }.values()
1060 )
1061 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1063 # Save plot.
1064 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1065 logging.info("Saved line plot to %s", filename)
1066 plt.close(fig)
1069def _plot_and_save_vertical_line_series(
1070 cubes: iris.cube.CubeList,
1071 coords: list[iris.coords.Coord],
1072 ensemble_coord: str,
1073 filename: str,
1074 series_coordinate: str,
1075 title: str,
1076 vmin: float,
1077 vmax: float,
1078 **kwargs,
1079):
1080 """Plot and save a 1D line series in vertical.
1082 Parameters
1083 ----------
1084 cubes: CubeList
1085 1 dimensional Cube or CubeList of the data to plot on x-axis.
1086 coord: list[Coord]
1087 Coordinates to plot on the y-axis, one per cube.
1088 ensemble_coord: str
1089 Ensemble coordinate in the cube.
1090 filename: str
1091 Filename of the plot to write.
1092 series_coordinate: str
1093 Coordinate to use as vertical axis.
1094 title: str
1095 Plot title.
1096 vmin: float
1097 Minimum value for the x-axis.
1098 vmax: float
1099 Maximum value for the x-axis.
1100 """
1101 # plot the vertical pressure axis using log scale
1102 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1104 model_colors_map = get_model_colors_map(cubes)
1106 # Check match-up across sequence coords gives consistent sizes
1107 validate_cubes_coords(cubes, coords)
1109 for cube, coord in zip(cubes, coords, strict=True):
1110 label = None
1111 color = "black"
1112 if model_colors_map: 1112 ↛ 1113line 1112 didn't jump to line 1113 because the condition on line 1112 was never true
1113 label = cube.attributes.get("model_name")
1114 color = model_colors_map.get(label)
1116 for cube_slice in cube.slices_over(ensemble_coord):
1117 # If ensemble data given plot control member with (control)
1118 # unless single forecast.
1119 if cube_slice.coord(ensemble_coord).points == [0]:
1120 iplt.plot(
1121 cube_slice,
1122 coord,
1123 color=color,
1124 marker="o",
1125 ls="-",
1126 lw=3,
1127 label=f"{label} (control)"
1128 if len(cube.coord(ensemble_coord).points) > 1
1129 else label,
1130 )
1131 # If ensemble data given plot perturbed members with (perturbed).
1132 else:
1133 iplt.plot(
1134 cube_slice,
1135 coord,
1136 color=color,
1137 ls="-",
1138 lw=1.5,
1139 alpha=0.75,
1140 label=f"{label} (member)",
1141 )
1143 # Get the current axis
1144 ax = plt.gca()
1146 # Special handling for pressure level data.
1147 if series_coordinate == "pressure": 1147 ↛ 1169line 1147 didn't jump to line 1169 because the condition on line 1147 was always true
1148 # Invert y-axis and set to log scale.
1149 ax.invert_yaxis()
1150 ax.set_yscale("log")
1152 # Define y-ticks and labels for pressure log axis.
1153 y_tick_labels = [
1154 "1000",
1155 "850",
1156 "700",
1157 "500",
1158 "300",
1159 "200",
1160 "100",
1161 ]
1162 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1164 # Set y-axis limits and ticks.
1165 ax.set_ylim(1100, 100)
1167 # Test if series_coordinate is model level data. The UM data uses
1168 # model_level_number and lfric uses full_levels as coordinate.
1169 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1170 # Define y-ticks and labels for vertical axis.
1171 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1172 y_tick_labels = [str(int(i)) for i in y_ticks]
1173 ax.set_ylim(min(y_ticks), max(y_ticks))
1175 ax.set_yticks(y_ticks)
1176 ax.set_yticklabels(y_tick_labels)
1178 # Set x-axis limits.
1179 ax.set_xlim(vmin, vmax)
1180 # Mark y=0 if present in plot.
1181 if vmin < 0.0 and vmax > 0.0: 1181 ↛ 1182line 1181 didn't jump to line 1182 because the condition on line 1181 was never true
1182 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1184 # Add some labels and tweak the style.
1185 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1186 ax.set_xlabel(
1187 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1188 )
1189 ax.set_title(title, fontsize=16)
1190 ax.ticklabel_format(axis="x")
1191 ax.tick_params(axis="y")
1192 ax.tick_params(axis="both", labelsize=12)
1194 # Add gridlines
1195 ax.grid(linestyle="--", color="grey", linewidth=1)
1196 # Ientify unique labels for legend
1197 handles = list(
1198 {
1199 label: handle
1200 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1201 }.values()
1202 )
1203 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1205 # Save plot.
1206 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1207 logging.info("Saved line plot to %s", filename)
1208 plt.close(fig)
1211def _plot_and_save_scatter_plot(
1212 cube_x: iris.cube.Cube | iris.cube.CubeList,
1213 cube_y: iris.cube.Cube | iris.cube.CubeList,
1214 filename: str,
1215 title: str,
1216 one_to_one: bool,
1217 model_names: list[str] = None,
1218 **kwargs,
1219):
1220 """Plot and save a 2D scatter plot.
1222 Parameters
1223 ----------
1224 cube_x: Cube | CubeList
1225 1 dimensional Cube or CubeList of the data to plot on x-axis.
1226 cube_y: Cube | CubeList
1227 1 dimensional Cube or CubeList of the data to plot on y-axis.
1228 filename: str
1229 Filename of the plot to write.
1230 title: str
1231 Plot title.
1232 one_to_one: bool
1233 Whether a 1:1 line is plotted.
1234 """
1235 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1236 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1237 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1238 # over the pairs simultaneously.
1240 # Ensure cube_x and cube_y are iterable
1241 cube_x_iterable = iter_maybe(cube_x)
1242 cube_y_iterable = iter_maybe(cube_y)
1244 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1245 iplt.scatter(cube_x_iter, cube_y_iter)
1246 if one_to_one is True:
1247 plt.plot(
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 [
1253 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1254 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1255 ],
1256 "k",
1257 linestyle="--",
1258 )
1259 ax = plt.gca()
1261 # Add some labels and tweak the style.
1262 if model_names is None:
1263 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1264 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1265 else:
1266 # Add the model names, these should be order of base (x) and other (y).
1267 ax.set_xlabel(
1268 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1269 )
1270 ax.set_ylabel(
1271 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1272 )
1273 ax.set_title(title, fontsize=16)
1274 ax.ticklabel_format(axis="y", useOffset=False)
1275 ax.tick_params(axis="x", labelrotation=15)
1276 ax.tick_params(axis="both", labelsize=12)
1277 ax.autoscale()
1279 # Save plot.
1280 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1281 logging.info("Saved scatter plot to %s", filename)
1282 plt.close(fig)
1285def _plot_and_save_vector_plot(
1286 cube_u: iris.cube.Cube,
1287 cube_v: iris.cube.Cube,
1288 filename: str,
1289 title: str,
1290 method: Literal["contourf", "pcolormesh"],
1291 **kwargs,
1292):
1293 """Plot and save a 2D vector plot.
1295 Parameters
1296 ----------
1297 cube_u: Cube
1298 2 dimensional Cube of u component of the data.
1299 cube_v: Cube
1300 2 dimensional Cube of v component of the data.
1301 filename: str
1302 Filename of the plot to write.
1303 title: str
1304 Plot title.
1305 """
1306 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1308 # Create a cube containing the magnitude of the vector field.
1309 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1310 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1312 # Specify the color bar
1313 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1315 # Setup plot map projection, extent and coastlines and borderlines.
1316 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1318 if method == "contourf":
1319 # Filled contour plot of the field.
1320 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1321 elif method == "pcolormesh":
1322 try:
1323 vmin = min(levels)
1324 vmax = max(levels)
1325 except TypeError:
1326 vmin, vmax = None, None
1327 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1328 # if levels are defined.
1329 if norm is not None:
1330 vmin = None
1331 vmax = None
1332 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1333 else:
1334 raise ValueError(f"Unknown plotting method: {method}")
1336 # Check to see if transect, and if so, adjust y axis.
1337 if is_transect(cube_vec_mag):
1338 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1339 axes.invert_yaxis()
1340 axes.set_yscale("log")
1341 axes.set_ylim(1100, 100)
1342 # If both model_level_number and level_height exists, iplt can construct
1343 # plot as a function of height above orography (NOT sea level).
1344 elif {"model_level_number", "level_height"}.issubset(
1345 {coord.name() for coord in cube_vec_mag.coords()}
1346 ):
1347 axes.set_yscale("log")
1349 axes.set_title(
1350 f"{title}\n"
1351 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1352 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1353 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1354 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1355 fontsize=16,
1356 )
1358 else:
1359 # Add title.
1360 axes.set_title(title, fontsize=16)
1362 # Add watermark with min/max/mean. Currently not user togglable.
1363 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1364 axes.annotate(
1365 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}",
1366 xy=(0.05, -0.05),
1367 xycoords="axes fraction",
1368 xytext=(-5, 5),
1369 textcoords="offset points",
1370 ha="right",
1371 va="bottom",
1372 size=11,
1373 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1374 )
1376 # Add colour bar.
1377 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1378 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1379 # add ticks and tick_labels for every levels if less than 20 levels exist
1380 if levels is not None and len(levels) < 20:
1381 cbar.set_ticks(levels)
1382 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1384 # 30 barbs along the longest axis of the plot, or a barb per point for data
1385 # with less than 30 points.
1386 step = max(max(cube_u.shape) // 30, 1)
1387 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1389 # Save plot.
1390 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1391 logging.info("Saved vector plot to %s", filename)
1392 plt.close(fig)
1395def _plot_and_save_histogram_series(
1396 cubes: iris.cube.Cube | iris.cube.CubeList,
1397 filename: str,
1398 title: str,
1399 vmin: float,
1400 vmax: float,
1401 **kwargs,
1402):
1403 """Plot and save a histogram series.
1405 Parameters
1406 ----------
1407 cubes: Cube or CubeList
1408 2 dimensional Cube or CubeList of the data to plot as histogram.
1409 filename: str
1410 Filename of the plot to write.
1411 title: str
1412 Plot title.
1413 vmin: float
1414 minimum for colorbar
1415 vmax: float
1416 maximum for colorbar
1417 """
1418 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1419 ax = plt.gca()
1421 model_colors_map = get_model_colors_map(cubes)
1423 # Set default that histograms will produce probability density function
1424 # at each bin (integral over range sums to 1).
1425 density = True
1427 for cube in iter_maybe(cubes):
1428 # Easier to check title (where var name originates)
1429 # than seeing if long names exist etc.
1430 # Exception case, where distribution better fits log scales/bins.
1431 if "surface_microphysical" in title:
1432 if "amount" in title: 1432 ↛ 1434line 1432 didn't jump to line 1434 because the condition on line 1432 was never true
1433 # Compute histogram following Klingaman et al. (2017): ASoP
1434 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1435 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1436 density = False
1437 else:
1438 bins = 10.0 ** (
1439 np.arange(-10, 27, 1) / 10.0
1440 ) # Suggestion from RMED toolbox.
1441 bins = np.insert(bins, 0, 0)
1442 ax.set_yscale("log")
1443 vmin = bins[1]
1444 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1445 ax.set_xscale("log")
1446 elif "lightning" in title:
1447 bins = [0, 1, 2, 3, 4, 5]
1448 else:
1449 bins = np.linspace(vmin, vmax, 51)
1450 logging.debug(
1451 "Plotting histogram with %s bins %s - %s.",
1452 np.size(bins),
1453 np.min(bins),
1454 np.max(bins),
1455 )
1457 # Reshape cube data into a single array to allow for a single histogram.
1458 # Otherwise we plot xdim histograms stacked.
1459 cube_data_1d = (cube.data).flatten()
1461 label = None
1462 color = "black"
1463 if model_colors_map: 1463 ↛ 1464line 1463 didn't jump to line 1464 because the condition on line 1463 was never true
1464 label = cube.attributes.get("model_name")
1465 color = model_colors_map[label]
1466 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1468 # Compute area under curve.
1469 if "surface_microphysical" in title and "amount" in title: 1469 ↛ 1470line 1469 didn't jump to line 1470 because the condition on line 1469 was never true
1470 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1471 x = x * bin_mean / x.sum()
1472 x = x[1:]
1473 y = y[1:]
1475 ax.plot(
1476 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1477 )
1479 # Add some labels and tweak the style.
1480 ax.set_title(title, fontsize=16)
1481 ax.set_xlabel(
1482 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1483 )
1484 ax.set_ylabel("Normalised probability density", fontsize=14)
1485 if "surface_microphysical" in title and "amount" in title: 1485 ↛ 1486line 1485 didn't jump to line 1486 because the condition on line 1485 was never true
1486 ax.set_ylabel(
1487 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1488 )
1489 ax.set_xlim(vmin, vmax)
1490 ax.tick_params(axis="both", labelsize=12)
1492 # Overlay grid-lines onto histogram plot.
1493 ax.grid(linestyle="--", color="grey", linewidth=1)
1494 if model_colors_map: 1494 ↛ 1495line 1494 didn't jump to line 1495 because the condition on line 1494 was never true
1495 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1497 # Save plot.
1498 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1499 logging.info("Saved histogram plot to %s", filename)
1500 plt.close(fig)
1503def _plot_and_save_postage_stamp_histogram_series(
1504 cube: iris.cube.Cube,
1505 filename: str,
1506 title: str,
1507 stamp_coordinate: str,
1508 vmin: float,
1509 vmax: float,
1510 **kwargs,
1511):
1512 """Plot and save postage (ensemble members) stamps for a histogram series.
1514 Parameters
1515 ----------
1516 cube: Cube
1517 2 dimensional Cube of the data to plot as histogram.
1518 filename: str
1519 Filename of the plot to write.
1520 title: str
1521 Plot title.
1522 stamp_coordinate: str
1523 Coordinate that becomes different plots.
1524 vmin: float
1525 minimum for pdf x-axis
1526 vmax: float
1527 maximum for pdf x-axis
1528 """
1529 # Use the smallest square grid that will fit the members.
1530 nmember = len(cube.coord(stamp_coordinate).points)
1531 grid_rows = int(math.sqrt(nmember))
1532 grid_size = math.ceil(nmember / grid_rows)
1534 fig = plt.figure(
1535 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1536 )
1537 # Make a subplot for each member.
1538 for member, subplot in zip(
1539 cube.slices_over(stamp_coordinate),
1540 range(1, grid_size * grid_rows + 1),
1541 strict=False,
1542 ):
1543 # Implicit interface is much easier here, due to needing to have the
1544 # cartopy GeoAxes generated.
1545 plt.subplot(grid_rows, grid_size, subplot)
1546 # Reshape cube data into a single array to allow for a single histogram.
1547 # Otherwise we plot xdim histograms stacked.
1548 member_data_1d = (member.data).flatten()
1549 plt.hist(member_data_1d, density=True, stacked=True)
1550 axes = plt.gca()
1551 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1552 axes.set_title(f"{mtitle}")
1553 axes.set_xlim(vmin, vmax)
1555 # Overall figure title.
1556 fig.suptitle(title, fontsize=16)
1558 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1559 logging.info("Saved histogram postage stamp plot to %s", filename)
1560 plt.close(fig)
1563def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1564 cube: iris.cube.Cube,
1565 filename: str,
1566 title: str,
1567 stamp_coordinate: str,
1568 vmin: float,
1569 vmax: float,
1570 **kwargs,
1571):
1572 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1573 ax.set_title(title, fontsize=16)
1574 ax.set_xlim(vmin, vmax)
1575 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1576 ax.set_ylabel("normalised probability density", fontsize=14)
1577 # Loop over all slices along the stamp_coordinate
1578 for member in cube.slices_over(stamp_coordinate):
1579 # Flatten the member data to 1D
1580 member_data_1d = member.data.flatten()
1581 # Plot the histogram using plt.hist
1582 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1583 plt.hist(
1584 member_data_1d,
1585 density=True,
1586 stacked=True,
1587 label=f"{mtitle}",
1588 )
1590 # Add a legend
1591 ax.legend(fontsize=16)
1593 # Save the figure to a file
1594 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1595 logging.info("Saved histogram postage stamp plot to %s", filename)
1597 # Close the figure
1598 plt.close(fig)
1601def _plot_and_save_scattermap_plot(
1602 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1603):
1604 """Plot and save a geographical scatter plot.
1606 Parameters
1607 ----------
1608 cube: Cube
1609 1 dimensional Cube of the data points with auxiliary latitude and
1610 longitude coordinates,
1611 filename: str
1612 Filename of the plot to write.
1613 title: str
1614 Plot title.
1615 projection: str
1616 Mapping projection to be used by cartopy.
1617 """
1618 # Setup plot details, size, resolution, etc.
1619 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1620 if projection is not None:
1621 # Apart from the default, the only projection we currently support is
1622 # a stereographic projection over the North Pole.
1623 if projection == "NP_Stereo":
1624 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1625 else:
1626 raise ValueError(f"Unknown projection: {projection}")
1627 else:
1628 axes = plt.axes(projection=ccrs.PlateCarree())
1630 # Scatter plot of the field. The marker size is chosen to give
1631 # symbols that decrease in size as the number of observations
1632 # increases, although the fraction of the figure covered by
1633 # symbols increases roughly as N^(1/2), disregarding overlaps,
1634 # and has been selected for the default figure size of (10, 10).
1635 # Should this be changed, the marker size should be adjusted in
1636 # proportion to the area of the figure.
1637 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1638 klon = None
1639 klat = None
1640 for kc in range(len(cube.aux_coords)):
1641 if cube.aux_coords[kc].standard_name == "latitude":
1642 klat = kc
1643 elif cube.aux_coords[kc].standard_name == "longitude":
1644 klon = kc
1645 scatter_map = iplt.scatter(
1646 cube.aux_coords[klon],
1647 cube.aux_coords[klat],
1648 c=cube.data[:],
1649 s=mrk_size,
1650 cmap="jet",
1651 edgecolors="k",
1652 )
1654 # Add coastlines and borderlines.
1655 try:
1656 axes.coastlines(resolution="10m")
1657 axes.add_feature(cfeature.BORDERS)
1658 except AttributeError:
1659 pass
1661 # Add title.
1662 axes.set_title(title, fontsize=16)
1664 # Add colour bar.
1665 cbar = fig.colorbar(scatter_map)
1666 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1668 # Save plot.
1669 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1670 logging.info("Saved geographical scatter plot to %s", filename)
1671 plt.close(fig)
1674def _spatial_plot(
1675 method: Literal["contourf", "pcolormesh"],
1676 cube: iris.cube.Cube,
1677 filename: str | None,
1678 sequence_coordinate: str,
1679 stamp_coordinate: str,
1680 overlay_cube: iris.cube.Cube | None = None,
1681 contour_cube: iris.cube.Cube | None = None,
1682 **kwargs,
1683):
1684 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1686 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1687 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1688 is present then postage stamp plots will be produced.
1690 If an overlay_cube and/or contour_cube are specified, multiple variables can
1691 be overplotted on the same figure.
1693 Parameters
1694 ----------
1695 method: "contourf" | "pcolormesh"
1696 The plotting method to use.
1697 cube: Cube
1698 Iris cube of the data to plot. It should have two spatial dimensions,
1699 such as lat and lon, and may also have a another two dimension to be
1700 plotted sequentially and/or as postage stamp plots.
1701 filename: str | None
1702 Name of the plot to write, used as a prefix for plot sequences. If None
1703 uses the recipe name.
1704 sequence_coordinate: str
1705 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1706 This coordinate must exist in the cube.
1707 stamp_coordinate: str
1708 Coordinate about which to plot postage stamp plots. Defaults to
1709 ``"realization"``.
1710 overlay_cube: Cube | None, optional
1711 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1712 contour_cube: Cube | None, optional
1713 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1715 Raises
1716 ------
1717 ValueError
1718 If the cube doesn't have the right dimensions.
1719 TypeError
1720 If the cube isn't a single cube.
1721 """
1722 recipe_title = get_recipe_metadata().get("title", "Untitled")
1724 # Ensure we've got a single cube.
1725 cube = check_single_cube(cube)
1727 # Check if there is a valid stamp coordinate in cube dimensions.
1728 if stamp_coordinate == "realization": 1728 ↛ 1733line 1728 didn't jump to line 1733 because the condition on line 1728 was always true
1729 stamp_coordinate = check_stamp_coordinate(cube)
1731 # Make postage stamp plots if stamp_coordinate exists and has more than a
1732 # single point.
1733 plotting_func = _plot_and_save_spatial_plot
1734 try:
1735 if cube.coord(stamp_coordinate).shape[0] > 1:
1736 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1737 except iris.exceptions.CoordinateNotFoundError:
1738 pass
1740 # Produce a geographical scatter plot if the data have a
1741 # dimension called observation or model_obs_error
1742 if any( 1742 ↛ 1746line 1742 didn't jump to line 1746 because the condition on line 1742 was never true
1743 crd.var_name == "station" or crd.var_name == "model_obs_error"
1744 for crd in cube.coords()
1745 ):
1746 plotting_func = _plot_and_save_scattermap_plot
1748 # Must have a sequence coordinate.
1749 try:
1750 cube.coord(sequence_coordinate)
1751 except iris.exceptions.CoordinateNotFoundError as err:
1752 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1754 # Create a plot for each value of the sequence coordinate.
1755 plot_index = []
1756 nplot = np.size(cube.coord(sequence_coordinate).points)
1758 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1759 # Set plot titles and filename
1760 seq_coord = cube_slice.coord(sequence_coordinate)
1761 plot_title, plot_filename = _set_title_and_filename(
1762 seq_coord, nplot, recipe_title, filename
1763 )
1765 # Extract sequence slice for overlay_cube and contour_cube if required.
1766 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1767 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1769 # Do the actual plotting.
1770 plotting_func(
1771 cube_slice,
1772 filename=plot_filename,
1773 stamp_coordinate=stamp_coordinate,
1774 title=plot_title,
1775 method=method,
1776 overlay_cube=overlay_slice,
1777 contour_cube=contour_slice,
1778 **kwargs,
1779 )
1780 plot_index.append(plot_filename)
1782 # Add list of plots to plot metadata.
1783 complete_plot_index = _append_to_plot_index(plot_index)
1785 # Make a page to display the plots.
1786 _make_plot_html_page(complete_plot_index)
1789####################
1790# Public functions #
1791####################
1794def spatial_contour_plot(
1795 cube: iris.cube.Cube,
1796 filename: str = None,
1797 sequence_coordinate: str = "time",
1798 stamp_coordinate: str = "realization",
1799 **kwargs,
1800) -> iris.cube.Cube:
1801 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1803 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1804 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1805 is present then postage stamp plots will be produced.
1807 Parameters
1808 ----------
1809 cube: Cube
1810 Iris cube of the data to plot. It should have two spatial dimensions,
1811 such as lat and lon, and may also have a another two dimension to be
1812 plotted sequentially and/or as postage stamp plots.
1813 filename: str, optional
1814 Name of the plot to write, used as a prefix for plot sequences. Defaults
1815 to the recipe name.
1816 sequence_coordinate: str, optional
1817 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1818 This coordinate must exist in the cube.
1819 stamp_coordinate: str, optional
1820 Coordinate about which to plot postage stamp plots. Defaults to
1821 ``"realization"``.
1823 Returns
1824 -------
1825 Cube
1826 The original cube (so further operations can be applied).
1828 Raises
1829 ------
1830 ValueError
1831 If the cube doesn't have the right dimensions.
1832 TypeError
1833 If the cube isn't a single cube.
1834 """
1835 _spatial_plot(
1836 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1837 )
1838 return cube
1841def spatial_pcolormesh_plot(
1842 cube: iris.cube.Cube,
1843 filename: str = None,
1844 sequence_coordinate: str = "time",
1845 stamp_coordinate: str = "realization",
1846 **kwargs,
1847) -> iris.cube.Cube:
1848 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1850 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1851 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1852 is present then postage stamp plots will be produced.
1854 This function is significantly faster than ``spatial_contour_plot``,
1855 especially at high resolutions, and should be preferred unless contiguous
1856 contour areas are important.
1858 Parameters
1859 ----------
1860 cube: Cube
1861 Iris cube of the data to plot. It should have two spatial dimensions,
1862 such as lat and lon, and may also have a another two dimension to be
1863 plotted sequentially and/or as postage stamp plots.
1864 filename: str, optional
1865 Name of the plot to write, used as a prefix for plot sequences. Defaults
1866 to the recipe name.
1867 sequence_coordinate: str, optional
1868 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1869 This coordinate must exist in the cube.
1870 stamp_coordinate: str, optional
1871 Coordinate about which to plot postage stamp plots. Defaults to
1872 ``"realization"``.
1874 Returns
1875 -------
1876 Cube
1877 The original cube (so further operations can be applied).
1879 Raises
1880 ------
1881 ValueError
1882 If the cube doesn't have the right dimensions.
1883 TypeError
1884 If the cube isn't a single cube.
1885 """
1886 _spatial_plot(
1887 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1888 )
1889 return cube
1892def spatial_multi_pcolormesh_plot(
1893 cube: iris.cube.Cube,
1894 overlay_cube: iris.cube.Cube,
1895 contour_cube: iris.cube.Cube,
1896 filename: str = None,
1897 sequence_coordinate: str = "time",
1898 stamp_coordinate: str = "realization",
1899 **kwargs,
1900) -> iris.cube.Cube:
1901 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1903 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1904 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1905 is present then postage stamp plots will be produced.
1907 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1909 If specified, contours of a contour_cube can be overplotted on top of those.
1911 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1913 This function is significantly faster than ``spatial_contour_plot``,
1914 especially at high resolutions, and should be preferred unless contiguous
1915 contour areas are important.
1917 Parameters
1918 ----------
1919 cube: Cube
1920 Iris cube of the data to plot. It should have two spatial dimensions,
1921 such as lat and lon, and may also have a another two dimension to be
1922 plotted sequentially and/or as postage stamp plots.
1923 overlay_cube: Cube
1924 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1925 such as lat and lon, and may also have a another two dimension to be
1926 plotted sequentially and/or as postage stamp plots. This is likely to be a masked cube in order not to hide the underlying basis cube.
1927 contour_cube: Cube
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.
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="contourf",
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)