Coverage for src/CSET/operators/plot.py: 80%
1026 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-24 13:03 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24import sys
25from typing import Literal
27import cartopy.crs as ccrs
28import cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.pyplot as plt
36import numpy as np
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._colormaps import (
50 colorbar_map_levels,
51 get_model_colors_map,
52)
53from CSET.operators._utils import (
54 check_sequence_coordinate,
55 check_single_cube,
56 check_stamp_coordinate,
57 fully_equalise_attributes,
58 get_cube_yxcoordname,
59 get_num_models,
60 is_transect,
61 slice_over_maybe,
62 validate_cube_shape,
63 validate_cubes_coords,
64)
65from CSET.operators.collapse import collapse
66from CSET.operators.misc import _extract_common_time_points
67from CSET.operators.regrid import regrid_onto_cube
69# Use a non-interactive plotting backend.
70mpl.use("agg")
73############################
74# Private helper functions #
75############################
78def in_sphinx_gallery():
79 """Test if running plot code in sphinx-gallery context."""
80 return "sphinx_gallery" in sys.modules
83def _append_to_plot_index(plot_index: list) -> list:
84 """Add plots into the plot index, returning the complete plot index."""
85 with open("meta.json", "r+t", encoding="UTF-8") as fp:
86 fcntl.flock(fp, fcntl.LOCK_EX)
87 fp.seek(0)
88 meta = json.load(fp)
89 complete_plot_index = meta.get("plots", [])
90 complete_plot_index = complete_plot_index + plot_index
91 meta["plots"] = complete_plot_index
92 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
93 os.getenv("DO_CASE_AGGREGATION")
94 ):
95 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
96 fp.seek(0)
97 fp.truncate()
98 json.dump(meta, fp, indent=2)
99 return complete_plot_index
102def _make_plot_html_page(plots: list):
103 """Create a HTML page to display a plot image."""
104 # Debug check that plots actually contains some strings.
105 assert isinstance(plots[0], str)
107 # Load HTML template file.
108 operator_files = importlib.resources.files()
109 template_file = operator_files.joinpath("_plot_page_template.html")
111 # Get some metadata.
112 meta = get_recipe_metadata()
113 title = meta.get("title", "Untitled")
114 description = MarkdownIt().render(meta.get("description", "*No description.*"))
116 # Prepare template variables.
117 variables = {
118 "title": title,
119 "description": description,
120 "initial_plot": plots[0],
121 "plots": plots,
122 "title_slug": slugify(title),
123 }
125 # Render template.
126 html = render_file(template_file, **variables)
128 # Save completed HTML.
129 with open("index.html", "wt", encoding="UTF-8") as fp:
130 fp.write(html)
133def _setup_spatial_map(
134 cube: iris.cube.Cube,
135 figure,
136 cmap,
137 grid_size: tuple[int, int] | None = None,
138 subplot: int | None = None,
139):
140 """Define map projections, extent and add coastlines and borderlines for spatial plots.
142 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
143 is specified, and map extent defined based on the input data.
145 Parameters
146 ----------
147 cube: Cube
148 2 dimensional (lat and lon) Cube of the data to plot.
149 figure:
150 Matplotlib Figure object holding all plot elements.
151 cmap:
152 Matplotlib colormap.
153 grid_size: (int, int), optional
154 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
155 subplot: int, optional
156 Subplot index if multiple spatial subplots in figure.
158 Returns
159 -------
160 axes:
161 Matplotlib GeoAxes definition.
162 """
163 # Identify min/max plot bounds.
164 try:
165 lat_axis, lon_axis = get_cube_yxcoordname(cube)
166 xmin = np.nanmin(cube.coord(lon_axis).points)
167 xmax = np.nanmax(cube.coord(lon_axis).points)
168 ymin = np.nanmin(cube.coord(lat_axis).points)
169 ymax = np.nanmax(cube.coord(lat_axis).points)
171 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
172 if np.abs(xmax - xmin) > 180.0:
173 xmin = xmin - 180.0
174 xmax = xmax - 180.0
175 logging.debug("Adjusting plot bounds to fit global extent.")
177 # Consider map projection orientation.
178 # Adapting orientation enables plotting across international dateline.
179 # Users can adapt the default central_longitude if alternative projections views.
180 if xmax > 180.0 or xmin < -180.0:
181 central_longitude = 180.0
182 else:
183 central_longitude = 0.0
185 # Define spatial map projection.
186 coord_system = cube.coord(lat_axis).coord_system
187 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
188 # Define rotated pole map projection for rotated pole inputs.
189 projection = ccrs.RotatedPole(
190 pole_longitude=coord_system.grid_north_pole_longitude,
191 pole_latitude=coord_system.grid_north_pole_latitude,
192 central_rotated_longitude=central_longitude,
193 )
194 crs = projection
195 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 195 ↛ 197line 195 didn't jump to line 197 because the condition on line 195 was never true
196 # Define Transverse Mercator projection for TM inputs.
197 projection = ccrs.TransverseMercator(
198 central_longitude=coord_system.longitude_of_central_meridian,
199 central_latitude=coord_system.latitude_of_projection_origin,
200 false_easting=coord_system.false_easting,
201 false_northing=coord_system.false_northing,
202 scale_factor=coord_system.scale_factor_at_central_meridian,
203 )
204 crs = projection
205 else:
206 # Assume polar projection for regional grids encompassing N. Pole
207 if ymin > 20.0 and ymax > 80.0:
208 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
209 elif ymin < -80.0 and ymax < -20.0:
210 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
211 # Define regular map projection for non-rotated pole inputs.
212 # Alternatives might include e.g. for global model outputs:
213 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
214 # projection = ccrs.NearsidePerspective(
215 # central_longitude=180.0,
216 # central_latitude=0,
217 # satellite_height=35785831,
218 # )
219 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
220 else:
221 projection = ccrs.PlateCarree(central_longitude=central_longitude)
222 crs = ccrs.PlateCarree()
224 # Define axes for plot (or subplot) with required map projection.
225 if subplot is not None:
226 axes = figure.add_subplot(
227 grid_size[0], grid_size[1], subplot, projection=projection
228 )
229 else:
230 axes = figure.add_subplot(projection=projection)
232 # Add coastlines and borderlines if cube contains x and y map coordinates.
233 # Avoid adding lines for 2D masked data or specific fixed ancillary spatial plots.
234 if (cube.ndim > 1 and iris.util.is_masked(cube.data)) or any(
235 name in cube.name() for name in ["land_", "orography", "altitude"]
236 ):
237 pass
238 else:
239 if cmap.name in ["viridis", "Greys"]:
240 coastcol = "magenta"
241 else:
242 coastcol = "black"
243 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
244 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
245 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
247 # Add gridlines.
248 gl = axes.gridlines(
249 alpha=0.3,
250 draw_labels=True,
251 dms=False,
252 x_inline=False,
253 y_inline=False,
254 )
255 gl.top_labels = False
256 gl.right_labels = False
257 if subplot:
258 gl.bottom_labels = False
259 gl.left_labels = False
260 if subplot % grid_size[1] == 1:
261 gl.left_labels = True
262 if subplot > ((grid_size[0] - 1) * grid_size[1]): 262 ↛ 267line 262 didn't jump to line 267 because the condition on line 262 was always true
263 gl.bottom_labels = True
265 # If is lat/lon spatial map, fix extent to keep plot tight.
266 # Specifying crs within set_extent helps ensure only data region is shown.
267 if isinstance(coord_system, iris.coord_systems.GeogCS):
268 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
270 except ValueError:
271 # Skip if not both x and y map coordinates.
272 axes = figure.gca()
273 pass
275 return axes
278def _get_plot_resolution() -> int:
279 """Get resolution of rasterised plots in pixels per inch."""
280 return get_recipe_metadata().get("plot_resolution", 100)
283def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
284 """Return title and filename based on start and end points or bounds."""
285 if use_bounds and seq_coord.has_bounds():
286 vals = seq_coord.bounds.flatten()
287 else:
288 vals = seq_coord.points
289 start = seq_coord.units.title(vals[0])
290 end = seq_coord.units.title(vals[-1])
292 if start == end:
293 sequence_title = f"\n [{start}]"
294 sequence_fname = f"_{filename_slugify(start)}"
295 else:
296 sequence_title = f"\n [{start} to {end}]"
297 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
299 # Do not include time if coord set to zero.
300 if (
301 seq_coord.units == "hours since 0001-01-01 00:00:00"
302 and vals[0] == 0
303 and vals[-1] == 0
304 ):
305 sequence_title = ""
306 sequence_fname = ""
308 return sequence_title, sequence_fname
311def _set_title_and_filename(
312 seq_coord: iris.coords.Coord,
313 nplot: int,
314 recipe_title: str,
315 filename: str,
316):
317 """Set plot title and filename based on cube coordinate.
319 Parameters
320 ----------
321 sequence_coordinate: iris.coords.Coord
322 Coordinate about which to make a plot sequence.
323 nplot: int
324 Number of output plots to generate - controls title/naming.
325 recipe_title: str
326 Default plot title, potentially to update.
327 filename: str
328 Input plot filename, potentially to update.
330 Returns
331 -------
332 plot_title: str
333 Output formatted plot title string, based on plotted data.
334 plot_filename: str
335 Output formatted plot filename string.
336 """
337 ndim = seq_coord.ndim
338 npoints = np.size(seq_coord.points)
339 sequence_title = ""
340 sequence_fname = ""
342 # Case 1: Multiple dimension sequence input - list number of aggregated cases
343 # (e.g. aggregation histogram plots)
344 if ndim > 1:
345 ncase = np.shape(seq_coord)[0]
346 sequence_title = f"\n [{ncase} cases]"
347 sequence_fname = f"_{ncase}cases"
349 # Case 2: Single dimension input
350 else:
351 # Single sequence point
352 if npoints == 1:
353 if nplot > 1:
354 # Default labels for sequence inputs
355 sequence_value = seq_coord.units.title(seq_coord.points[0])
356 sequence_title = f"\n [{sequence_value}]"
357 sequence_fname = f"_{filename_slugify(sequence_value)}"
358 else:
359 # Aggregated attribute available where input collapsed over aggregation
360 try:
361 ncase = seq_coord.attributes["number_reference_times"]
362 sequence_title = f"\n [{ncase} cases]"
363 sequence_fname = f"_{ncase}cases"
364 except KeyError:
365 sequence_title, sequence_fname = _get_start_end_strings(
366 seq_coord, use_bounds=seq_coord.has_bounds()
367 )
368 # Multiple sequence (e.g. time) points
369 else:
370 sequence_title, sequence_fname = _get_start_end_strings(
371 seq_coord, use_bounds=False
372 )
374 # Set plot title and filename
375 plot_title = f"{recipe_title}{sequence_title}"
377 # Set plot filename, defaulting to user input if provided.
378 if filename is None:
379 filename = slugify(recipe_title)
380 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
381 else:
382 if nplot > 1:
383 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
384 else:
385 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
387 return plot_title, plot_filename
390def _select_series_coord(cube, series_coordinate):
391 """Determine the grid coordinates to use to calculate grid spacing."""
392 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
393 if series_coordinate in spacing_coordinates: 393 ↛ 399line 393 didn't jump to line 399 because the condition on line 393 was always true
394 # Try the requested coordinate first then the fallbacks in order.
395 fallbacks = [series_coordinate] + [
396 c for c in spacing_coordinates if c != series_coordinate
397 ]
398 else:
399 fallbacks = {series_coordinate}
401 # Try each possible coordinate.
402 for coord in fallbacks:
403 try:
404 return cube.coord(coord)
405 except iris.exceptions.CoordinateNotFoundError:
406 logging.debug("Coordinate %s not found.", coord)
408 # If we get here, none of the fallback options were found.
409 raise iris.exceptions.CoordinateNotFoundError(
410 f"No valid coordinate found for '{series_coordinate}' "
411 f"or fallback options {fallbacks}"
412 )
415def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
416 """Control postage stamp plot output titles based on stamp coordinate."""
417 if stamp_coord.name() == "realization":
418 mtitle = "Member"
419 else:
420 mtitle = stamp_coord.name().capitalize()
422 if stamp_coord.name() == "time":
423 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
424 else:
425 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
427 return mtitle
430def _set_axis_range(cubes):
431 """Get minimum and maximum from levels information."""
432 levels = None
433 for cube in cubes: 433 ↛ 449line 433 didn't jump to line 449 because the loop on line 433 didn't complete
434 # First check if user-specified "auto" range variable.
435 # This maintains the value of levels as None, so proceed.
436 _, levels, _ = colorbar_map_levels(cube, axis="y")
437 if levels is None:
438 break
439 # If levels is changed, recheck to use the vmin,vmax or
440 # levels-based ranges for histogram plots.
441 _, levels, _ = colorbar_map_levels(cube)
442 logging.debug("levels: %s", levels)
443 if levels is not None: 443 ↛ 433line 443 didn't jump to line 433 because the condition on line 443 was always true
444 vmin = min(levels)
445 vmax = max(levels)
446 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
447 break
449 if levels is None:
450 vmin = min(cb.data.min() for cb in cubes)
451 vmax = max(cb.data.max() for cb in cubes)
453 return vmin, vmax
456def _find_matched_slices(cubes, sequence_coordinate):
457 """Identify matched cubes in CubeList by sequence_coordinate values.
459 Ensures common points are compared for multiple cube inputs.
460 """
461 all_points = sorted(
462 set(
463 itertools.chain.from_iterable(
464 cb.coord(sequence_coordinate).points for cb in cubes
465 )
466 )
467 )
468 all_slices = list(
469 itertools.chain.from_iterable(
470 cb.slices_over(sequence_coordinate) for cb in cubes
471 )
472 )
473 # Matched slices (matched by seq coord point; it may happen that
474 # evaluated models do not cover the same seq coord range, hence matching
475 # necessary)
476 cube_iterables = [
477 iris.cube.CubeList(
478 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
479 )
480 for point in all_points
481 ]
483 return cube_iterables
486def _plot_and_save_spatial_plot(
487 cube: iris.cube.Cube,
488 filename: str,
489 title: str,
490 method: Literal["contourf", "pcolormesh", "scatter"],
491 overlay_cube: iris.cube.Cube | None = None,
492 contour_cube: iris.cube.Cube | None = None,
493 point_cube: iris.cube.Cube | None = None,
494 **kwargs,
495):
496 """Plot and save a spatial plot.
498 Parameters
499 ----------
500 cube: Cube
501 2 dimensional (lat and lon) Cube of the data to plot.
502 filename: str
503 Filename of the plot to write.
504 title: str
505 Plot title.
506 method: "contourf" | "pcolormesh" | "scatter"
507 The plotting method to use
508 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
509 overlay_cube: Cube, optional
510 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
511 contour_cube: Cube, optional
512 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
513 point_cube: Cube, optional
514 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
515 """
516 # Setup plot details, size, resolution, etc.
517 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
519 # Specify the color bar
520 cmap, levels, norm = colorbar_map_levels(cube)
522 # If overplotting, set required colorbars
523 if overlay_cube:
524 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
525 if contour_cube:
526 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
528 # Setup plot map projection, extent and coastlines and borderlines.
529 axes = _setup_spatial_map(cube, fig, cmap)
531 # Set colorscale bounds
532 try:
533 vmin = min(levels)
534 vmax = max(levels)
535 except TypeError:
536 vmin, vmax = None, None
537 # Ensure to use norm and not vmin/vmax if levels are defined.
538 if norm is not None:
539 vmin = None
540 vmax = None
541 logging.debug("Plotting using defined levels.")
543 # Plot the field.
544 if method == "contourf":
545 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
546 elif method == "pcolormesh":
547 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
548 elif method == "scatter": 548 ↛ 556line 548 didn't jump to line 556 because the condition on line 548 was never true
549 # Scatter plot of the field. The marker size is chosen to give
550 # symbols that decrease in size as the number of data points
551 # increases, although the fraction of the figure covered by
552 # symbols increases roughly as N^(1/2), disregarding overlaps,
553 # and has been selected for the default figure size of (10, 10).
554 # Should this be changed, the marker size should be adjusted in
555 # proportion to the area of the figure.
556 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
557 lat_axis, lon_axis = get_cube_yxcoordname(cube)
558 plot = iplt.scatter(
559 cube.coord(lon_axis),
560 cube.coord(lat_axis),
561 c=cube.data[:],
562 s=mrk_size,
563 cmap=cmap,
564 edgecolors="k",
565 norm=norm,
566 vmin=vmin,
567 vmax=vmax,
568 )
569 else:
570 raise ValueError(f"Unknown plotting method: {method}")
572 # Overplot overlay field, if required
573 if overlay_cube:
574 try:
575 over_vmin = min(over_levels)
576 over_vmax = max(over_levels)
577 except TypeError:
578 over_vmin, over_vmax = None, None
579 if over_norm is not None: 579 ↛ 580line 579 didn't jump to line 580 because the condition on line 579 was never true
580 over_vmin = None
581 over_vmax = None
582 overlay = iplt.pcolormesh(
583 overlay_cube,
584 cmap=over_cmap,
585 norm=over_norm,
586 alpha=0.8,
587 vmin=over_vmin,
588 vmax=over_vmax,
589 )
590 # Overplot contour field, if required, with contour labelling.
591 if contour_cube:
592 contour = iplt.contour(
593 contour_cube,
594 colors="darkgray",
595 levels=cntr_levels,
596 norm=cntr_norm,
597 alpha=0.5,
598 linestyles="--",
599 linewidths=1,
600 )
601 plt.clabel(contour)
602 # Overplot valid elements of point-based field, if required.
603 # Check for non-masked points only to avoid plotting missing data.
604 if point_cube: 604 ↛ 605line 604 didn't jump to line 605 because the condition on line 604 was never true
605 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
606 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
607 lon_coord = point_cube.coord(lon_axis)
608 lat_coord = point_cube.coord(lat_axis)
609 valid = ~point_cube.data.mask
610 valid_lon = iris.coords.AuxCoord(
611 lon_coord.points[valid],
612 standard_name=lon_coord.standard_name,
613 units=lon_coord.units,
614 coord_system=lon_coord.coord_system,
615 )
616 valid_lat = iris.coords.AuxCoord(
617 lat_coord.points[valid],
618 standard_name=lat_coord.standard_name,
619 units=lat_coord.units,
620 coord_system=lat_coord.coord_system,
621 )
622 iplt.scatter(
623 valid_lon,
624 valid_lat,
625 c=point_cube.data[valid],
626 s=mrk_size,
627 cmap=cmap,
628 edgecolors="k",
629 norm=norm,
630 vmin=vmin,
631 vmax=vmax,
632 )
634 # Check to see if transect, and if so, adjust y axis.
635 if is_transect(cube):
636 if "pressure" in [coord.name() for coord in cube.coords()]:
637 axes.invert_yaxis()
638 axes.set_yscale("log")
639 axes.set_ylim(1100, 100)
640 # If both model_level_number and level_height exists, iplt can construct
641 # plot as a function of height above orography (NOT sea level).
642 elif {"model_level_number", "level_height"}.issubset( 642 ↛ 647line 642 didn't jump to line 647 because the condition on line 642 was always true
643 {coord.name() for coord in cube.coords()}
644 ):
645 axes.set_yscale("log")
647 axes.set_title(
648 f"{title}\n"
649 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
650 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
651 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
652 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
653 fontsize=16,
654 )
656 # Inset code
657 axins = inset_axes(
658 axes,
659 width="20%",
660 height="20%",
661 loc="upper right",
662 axes_class=GeoAxes,
663 axes_kwargs={"map_projection": ccrs.PlateCarree()},
664 )
666 # Slightly transparent to reduce plot blocking.
667 axins.patch.set_alpha(0.4)
669 axins.coastlines(resolution="50m")
670 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
672 SLat, SLon, ELat, ELon = (
673 float(coord) for coord in cube.attributes["transect_coords"].split("_")
674 )
676 # Draw line between them
677 axins.plot(
678 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
679 )
681 # Plot points (note: lon, lat order for Cartopy)
682 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
683 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
685 lon_min, lon_max = sorted([SLon, ELon])
686 lat_min, lat_max = sorted([SLat, ELat])
688 # Midpoints
689 lon_mid = (lon_min + lon_max) / 2
690 lat_mid = (lat_min + lat_max) / 2
692 # Maximum half-range
693 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
694 if half_range == 0: # points identical → provide small default 694 ↛ 698line 694 didn't jump to line 698 because the condition on line 694 was always true
695 half_range = 1
697 # Set square extent
698 axins.set_extent(
699 [
700 lon_mid - half_range,
701 lon_mid + half_range,
702 lat_mid - half_range,
703 lat_mid + half_range,
704 ],
705 crs=ccrs.PlateCarree(),
706 )
708 # Ensure square aspect
709 axins.set_aspect("equal")
711 else:
712 # Add title.
713 axes.set_title(title, fontsize=16)
715 # Adjust padding if spatial plot or transect
716 if is_transect(cube):
717 yinfopad = -0.1
718 ycbarpad = 0.1
719 else:
720 yinfopad = 0.01
721 ycbarpad = 0.042
723 # Add watermark with min/max/mean. Currently not user togglable.
724 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
725 axes.annotate(
726 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
727 xy=(0.025, yinfopad),
728 xycoords="axes fraction",
729 xytext=(-5, 5),
730 textcoords="offset points",
731 ha="left",
732 va="bottom",
733 size=11,
734 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
735 )
737 # Add secondary colour bar for overlay_cube field if required.
738 if overlay_cube:
739 cbarB = fig.colorbar(
740 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
741 )
742 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
743 # add ticks and tick_labels for every levels if less than 20 levels exist
744 if over_levels is not None and len(over_levels) < 20: 744 ↛ 745line 744 didn't jump to line 745 because the condition on line 744 was never true
745 cbarB.set_ticks(over_levels)
746 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
747 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
748 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
749 logging.debug("Set secondary colorbar ticks and labels.")
751 # Add main colour bar.
752 cbar = fig.colorbar(
753 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
754 )
756 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
757 # add ticks and tick_labels for every levels if less than 20 levels exist
758 if levels is not None and len(levels) < 20:
759 cbar.set_ticks(levels)
760 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
761 if "rainfall" or "snowfall" or "visibility" in cube.name(): 761 ↛ 764line 761 didn't jump to line 764 because the condition on line 761 was always true
762 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
763 # Tick labels for rainfall rates from Nimrod radar data.
764 if "rainfall rate composite" in cube.name(): 764 ↛ 765line 764 didn't jump to line 765 because the condition on line 764 was never true
765 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
766 # Tick labels for rain accumulations from Nimrod radar data.
767 if "rain accumulation" in cube.name(): 767 ↛ 768line 767 didn't jump to line 768 because the condition on line 767 was never true
768 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
769 if "wts accumulation" in cube.name(): 769 ↛ 770line 769 didn't jump to line 770 because the condition on line 769 was never true
770 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
771 cbar.minorticks_off()
772 cbar.set_ticks(tick_levels)
773 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
774 cbar.set_label(label=f"{cube.name()}", size=14)
775 # Tick labels for model rainfall data.
776 if "surface_microphysical" in cube.name(): 776 ↛ 779line 776 didn't jump to line 779 because the condition on line 776 was always true
777 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
778 # Tick labels for Nimrod weights data.
779 logging.debug("Set colorbar ticks and labels.")
781 # Save plot.
782 if not in_sphinx_gallery(): 782 ↛ exitline 782 didn't return from function '_plot_and_save_spatial_plot' because the condition on line 782 was always true
783 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
784 logging.info("Saved spatial plot to %s", filename)
785 plt.close(fig)
788def _plot_and_save_postage_stamp_spatial_plot(
789 cube: iris.cube.Cube,
790 filename: str,
791 stamp_coordinate: str,
792 title: str,
793 method: Literal["contourf", "pcolormesh"],
794 overlay_cube: iris.cube.Cube | None = None,
795 contour_cube: iris.cube.Cube | None = None,
796 **kwargs,
797):
798 """Plot postage stamp spatial plots from an ensemble.
800 Parameters
801 ----------
802 cube: Cube
803 Iris cube of data to be plotted. It must have the stamp coordinate.
804 filename: str
805 Filename of the plot to write.
806 stamp_coordinate: str
807 Coordinate that becomes different plots.
808 method: "contourf" | "pcolormesh"
809 The plotting method to use.
810 overlay_cube: Cube, optional
811 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
812 contour_cube: Cube, optional
813 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
815 Raises
816 ------
817 ValueError
818 If the cube doesn't have the right dimensions.
819 """
820 # Use the smallest square grid that will fit the members.
821 nmember = len(cube.coord(stamp_coordinate).points)
822 grid_rows = int(math.sqrt(nmember))
823 grid_size = math.ceil(nmember / grid_rows)
825 fig = plt.figure(
826 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
827 )
829 # Specify the color bar
830 cmap, levels, norm = colorbar_map_levels(cube)
831 # If overplotting, set required colorbars
832 if overlay_cube: 832 ↛ 833line 832 didn't jump to line 833 because the condition on line 832 was never true
833 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
834 if contour_cube: 834 ↛ 835line 834 didn't jump to line 835 because the condition on line 834 was never true
835 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
837 # Make a subplot for each member.
838 for member, subplot in zip(
839 cube.slices_over(stamp_coordinate),
840 range(1, grid_size * grid_rows + 1),
841 strict=False,
842 ):
843 # Setup subplot map projection, extent and coastlines and borderlines.
844 axes = _setup_spatial_map(
845 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
846 )
847 if method == "contourf":
848 # Filled contour plot of the field.
849 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
850 elif method == "pcolormesh":
851 if levels is not None:
852 vmin = min(levels)
853 vmax = max(levels)
854 else:
855 raise TypeError("Unknown vmin and vmax range.")
856 vmin, vmax = None, None
857 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
858 # if levels are defined.
859 if norm is not None: 859 ↛ 860line 859 didn't jump to line 860 because the condition on line 859 was never true
860 vmin = None
861 vmax = None
862 # pcolormesh plot of the field.
863 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
864 else:
865 raise ValueError(f"Unknown plotting method: {method}")
867 # Overplot overlay field, if required
868 if overlay_cube: 868 ↛ 869line 868 didn't jump to line 869 because the condition on line 868 was never true
869 try:
870 over_vmin = min(over_levels)
871 over_vmax = max(over_levels)
872 except TypeError:
873 over_vmin, over_vmax = None, None
874 if over_norm is not None:
875 over_vmin = None
876 over_vmax = None
877 iplt.pcolormesh(
878 overlay_cube[member.coord(stamp_coordinate).points[0]],
879 cmap=over_cmap,
880 norm=over_norm,
881 alpha=0.6,
882 vmin=over_vmin,
883 vmax=over_vmax,
884 )
885 # Overplot contour field, if required
886 if contour_cube: 886 ↛ 887line 886 didn't jump to line 887 because the condition on line 886 was never true
887 iplt.contour(
888 contour_cube[member.coord(stamp_coordinate).points[0]],
889 colors="darkgray",
890 levels=cntr_levels,
891 norm=cntr_norm,
892 alpha=0.6,
893 linestyles="--",
894 linewidths=1,
895 )
896 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
897 axes.set_title(f"{mtitle}")
899 # Put the shared colorbar in its own axes.
900 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
901 colorbar = fig.colorbar(
902 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
903 )
904 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
906 # Overall figure title.
907 fig.suptitle(title, fontsize=16)
908 if not in_sphinx_gallery(): 908 ↛ exitline 908 didn't return from function '_plot_and_save_postage_stamp_spatial_plot' because the condition on line 908 was always true
909 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
910 logging.info("Saved contour postage stamp plot to %s", filename)
911 plt.close(fig)
914def _plot_and_save_line_series(
915 cubes: iris.cube.CubeList,
916 coords: list[iris.coords.Coord],
917 ensemble_coord: str,
918 filename: str,
919 title: str,
920 **kwargs,
921):
922 """Plot and save a 1D line series.
924 Parameters
925 ----------
926 cubes: Cube or CubeList
927 Cube or CubeList containing the cubes to plot on the y-axis.
928 coords: list[Coord]
929 Coordinates to plot on the x-axis, one per cube.
930 ensemble_coord: str
931 Ensemble coordinate in the cube.
932 filename: str
933 Filename of the plot to write.
934 title: str
935 Plot title.
936 """
937 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
939 model_colors_map = get_model_colors_map(cubes)
941 # Store min/max ranges.
942 y_levels = []
944 # Check match-up across sequence coords gives consistent sizes
945 validate_cubes_coords(cubes, coords)
947 for cube, coord in zip(cubes, coords, strict=True):
948 label = None
949 color = "black"
950 if model_colors_map:
951 label = cube.attributes.get("model_name")
952 color = model_colors_map.get(label)
953 for cube_slice in cube.slices_over(ensemble_coord):
954 # Label with (control) if part of an ensemble or not otherwise.
955 if cube_slice.coord(ensemble_coord).points == [0]:
956 iplt.plot(
957 coord,
958 cube_slice,
959 color=color,
960 marker="o",
961 ls="-",
962 lw=3,
963 label=f"{label} (control)"
964 if len(cube.coord(ensemble_coord).points) > 1
965 else label,
966 )
967 # Label with (perturbed) if part of an ensemble and not the control.
968 else:
969 iplt.plot(
970 coord,
971 cube_slice,
972 color=color,
973 ls="-",
974 lw=1.5,
975 alpha=0.75,
976 label=f"{label} (member)",
977 )
979 # Calculate the global min/max if multiple cubes are given.
980 _, levels, _ = colorbar_map_levels(cube, axis="y")
981 if levels is not None: 981 ↛ 982line 981 didn't jump to line 982 because the condition on line 981 was never true
982 y_levels.append(min(levels))
983 y_levels.append(max(levels))
985 # Get the current axes.
986 ax = plt.gca()
988 # Add some labels and tweak the style.
989 # check if cubes[0] works for single cube if not CubeList
990 if coords[0].name() == "time":
991 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
992 else:
993 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
994 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
995 ax.set_title(title, fontsize=16)
997 ax.ticklabel_format(axis="y", useOffset=False)
998 ax.tick_params(axis="x", labelrotation=15)
999 ax.tick_params(axis="both", labelsize=12)
1001 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1002 if y_levels: 1002 ↛ 1003line 1002 didn't jump to line 1003 because the condition on line 1002 was never true
1003 ax.set_ylim(min(y_levels), max(y_levels))
1004 # Add zero line.
1005 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1006 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1007 logging.debug(
1008 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1009 )
1010 else:
1011 ax.autoscale()
1013 # Add gridlines
1014 ax.grid(linestyle="--", color="grey", linewidth=1)
1015 # Ientify unique labels for legend
1016 handles = list(
1017 {
1018 label: handle
1019 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1020 }.values()
1021 )
1022 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1024 # Save plot.
1025 if not in_sphinx_gallery(): 1025 ↛ exitline 1025 didn't return from function '_plot_and_save_line_series' because the condition on line 1025 was always true
1026 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1027 logging.info("Saved line plot to %s", filename)
1028 plt.close(fig)
1031def _plot_and_save_line_power_spectrum_series(
1032 cubes: iris.cube.Cube | iris.cube.CubeList,
1033 coords: list[iris.coords.Coord],
1034 ensemble_coord: str,
1035 filename: str,
1036 title: str,
1037 series_coordinate: str,
1038 **kwargs,
1039):
1040 """Plot and save a 1D line series.
1042 Parameters
1043 ----------
1044 cubes: Cube or CubeList
1045 Cube or CubeList containing the cubes to plot on the y-axis.
1046 coords: list[Coord]
1047 Coordinates to plot on the x-axis, one per cube.
1048 ensemble_coord: str
1049 Ensemble coordinate in the cube.
1050 filename: str
1051 Filename of the plot to write.
1052 title: str
1053 Plot title.
1054 series_coordinate: str
1055 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1056 """
1057 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1058 model_colors_map = get_model_colors_map(cubes)
1059 ax = plt.gca()
1061 # Store min/max ranges.
1062 y_levels = []
1064 line_marker = None
1065 line_width = 1
1067 for cube in iter_maybe(cubes):
1068 # next 2 lines replace chunk of code.
1069 xcoord = _select_series_coord(cube, series_coordinate)
1070 xname = xcoord.points
1072 yfield = cube.data # power spectrum
1073 label = None
1074 color = "black"
1075 if model_colors_map: 1075 ↛ 1078line 1075 didn't jump to line 1078 because the condition on line 1075 was always true
1076 label = cube.attributes.get("model_name")
1077 color = model_colors_map.get(label)
1078 for cube_slice in cube.slices_over(ensemble_coord):
1079 # Label with (control) if part of an ensemble or not otherwise.
1080 if cube_slice.coord(ensemble_coord).points == [0]: 1080 ↛ 1094line 1080 didn't jump to line 1094 because the condition on line 1080 was always true
1081 ax.plot(
1082 xname,
1083 yfield,
1084 color=color,
1085 marker=line_marker,
1086 ls="-",
1087 lw=line_width,
1088 label=f"{label} (control)"
1089 if len(cube.coord(ensemble_coord).points) > 1
1090 else label,
1091 )
1092 # Label with (perturbed) if part of an ensemble and not the control.
1093 else:
1094 ax.plot(
1095 xname,
1096 yfield,
1097 color=color,
1098 ls="-",
1099 lw=1.5,
1100 alpha=0.75,
1101 label=f"{label} (member)",
1102 )
1104 # Calculate the global min/max if multiple cubes are given.
1105 _, levels, _ = colorbar_map_levels(cube, axis="y")
1106 if levels is not None: 1106 ↛ 1107line 1106 didn't jump to line 1107 because the condition on line 1106 was never true
1107 y_levels.append(min(levels))
1108 y_levels.append(max(levels))
1110 # Add some labels and tweak the style.
1112 title = f"{title}"
1113 ax.set_title(title, fontsize=16)
1115 # Set appropriate x-axis label based on coordinate
1116 if series_coordinate == "wavelength" or ( 1116 ↛ 1119line 1116 didn't jump to line 1119 because the condition on line 1116 was never true
1117 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1118 ):
1119 ax.set_xlabel("Wavelength (km)", fontsize=14)
1120 elif series_coordinate == "physical_wavenumber" or ( 1120 ↛ 1123line 1120 didn't jump to line 1123 because the condition on line 1120 was never true
1121 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1122 ):
1123 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1124 else: # frequency or check units
1125 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1125 ↛ 1126line 1125 didn't jump to line 1126 because the condition on line 1125 was never true
1126 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1127 else:
1128 ax.set_xlabel("Wavenumber", fontsize=14)
1130 ax.set_ylabel("Power Spectral Density", fontsize=14)
1131 ax.tick_params(axis="both", labelsize=12)
1133 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1135 # Set log-log scale
1136 ax.set_xscale("log")
1137 ax.set_yscale("log")
1139 # Add gridlines
1140 ax.grid(linestyle="--", color="grey", linewidth=1)
1141 # Ientify unique labels for legend
1142 handles = list(
1143 {
1144 label: handle
1145 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1146 }.values()
1147 )
1148 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1150 # Save plot.
1151 if not in_sphinx_gallery(): 1151 ↛ exitline 1151 didn't return from function '_plot_and_save_line_power_spectrum_series' because the condition on line 1151 was always true
1152 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1153 logging.info("Saved line plot to %s", filename)
1154 plt.close(fig)
1157def _plot_and_save_vertical_line_series(
1158 cubes: iris.cube.CubeList,
1159 coords: list[iris.coords.Coord],
1160 ensemble_coord: str,
1161 filename: str,
1162 series_coordinate: str,
1163 title: str,
1164 vmin: float,
1165 vmax: float,
1166 **kwargs,
1167):
1168 """Plot and save a 1D line series in vertical.
1170 Parameters
1171 ----------
1172 cubes: CubeList
1173 1 dimensional Cube or CubeList of the data to plot on x-axis.
1174 coord: list[Coord]
1175 Coordinates to plot on the y-axis, one per cube.
1176 ensemble_coord: str
1177 Ensemble coordinate in the cube.
1178 filename: str
1179 Filename of the plot to write.
1180 series_coordinate: str
1181 Coordinate to use as vertical axis.
1182 title: str
1183 Plot title.
1184 vmin: float
1185 Minimum value for the x-axis.
1186 vmax: float
1187 Maximum value for the x-axis.
1188 """
1189 # plot the vertical pressure axis using log scale
1190 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1192 model_colors_map = get_model_colors_map(cubes)
1194 # Check match-up across sequence coords gives consistent sizes
1195 validate_cubes_coords(cubes, coords)
1197 for cube, coord in zip(cubes, coords, strict=True):
1198 label = None
1199 color = "black"
1200 if model_colors_map: 1200 ↛ 1201line 1200 didn't jump to line 1201 because the condition on line 1200 was never true
1201 label = cube.attributes.get("model_name")
1202 color = model_colors_map.get(label)
1204 for cube_slice in cube.slices_over(ensemble_coord):
1205 # If ensemble data given plot control member with (control)
1206 # unless single forecast.
1207 if cube_slice.coord(ensemble_coord).points == [0]:
1208 iplt.plot(
1209 cube_slice,
1210 coord,
1211 color=color,
1212 marker="o",
1213 ls="-",
1214 lw=3,
1215 label=f"{label} (control)"
1216 if len(cube.coord(ensemble_coord).points) > 1
1217 else label,
1218 )
1219 # If ensemble data given plot perturbed members with (perturbed).
1220 else:
1221 iplt.plot(
1222 cube_slice,
1223 coord,
1224 color=color,
1225 ls="-",
1226 lw=1.5,
1227 alpha=0.75,
1228 label=f"{label} (member)",
1229 )
1231 # Get the current axis
1232 ax = plt.gca()
1234 # Special handling for pressure level data.
1235 if series_coordinate == "pressure": 1235 ↛ 1257line 1235 didn't jump to line 1257 because the condition on line 1235 was always true
1236 # Invert y-axis and set to log scale.
1237 ax.invert_yaxis()
1238 ax.set_yscale("log")
1240 # Define y-ticks and labels for pressure log axis.
1241 y_tick_labels = [
1242 "1000",
1243 "850",
1244 "700",
1245 "500",
1246 "300",
1247 "200",
1248 "100",
1249 ]
1250 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1252 # Set y-axis limits and ticks.
1253 ax.set_ylim(1100, 100)
1255 # Test if series_coordinate is model level data. The UM data uses
1256 # model_level_number and lfric uses full_levels as coordinate.
1257 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1258 # Define y-ticks and labels for vertical axis.
1259 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1260 y_tick_labels = [str(int(i)) for i in y_ticks]
1261 ax.set_ylim(min(y_ticks), max(y_ticks))
1263 ax.set_yticks(y_ticks)
1264 ax.set_yticklabels(y_tick_labels)
1266 # Set x-axis limits.
1267 ax.set_xlim(vmin, vmax)
1268 # Mark y=0 if present in plot.
1269 if vmin < 0.0 and vmax > 0.0: 1269 ↛ 1270line 1269 didn't jump to line 1270 because the condition on line 1269 was never true
1270 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1272 # Add some labels and tweak the style.
1273 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1274 ax.set_xlabel(
1275 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1276 )
1277 ax.set_title(title, fontsize=16)
1278 ax.ticklabel_format(axis="x")
1279 ax.tick_params(axis="y")
1280 ax.tick_params(axis="both", labelsize=12)
1282 # Add gridlines
1283 ax.grid(linestyle="--", color="grey", linewidth=1)
1284 # Ientify unique labels for legend
1285 handles = list(
1286 {
1287 label: handle
1288 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1289 }.values()
1290 )
1291 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1293 # Save plot.
1294 if not in_sphinx_gallery(): 1294 ↛ exitline 1294 didn't return from function '_plot_and_save_vertical_line_series' because the condition on line 1294 was always true
1295 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1296 logging.info("Saved line plot to %s", filename)
1297 plt.close(fig)
1300def _plot_and_save_scatter_plot(
1301 cube_x: iris.cube.Cube | iris.cube.CubeList,
1302 cube_y: iris.cube.Cube | iris.cube.CubeList,
1303 filename: str,
1304 title: str,
1305 one_to_one: bool,
1306 model_names: list[str] = None,
1307 **kwargs,
1308):
1309 """Plot and save a 2D scatter plot.
1311 Parameters
1312 ----------
1313 cube_x: Cube | CubeList
1314 1 dimensional Cube or CubeList of the data to plot on x-axis.
1315 cube_y: Cube | CubeList
1316 1 dimensional Cube or CubeList of the data to plot on y-axis.
1317 filename: str
1318 Filename of the plot to write.
1319 title: str
1320 Plot title.
1321 one_to_one: bool
1322 Whether a 1:1 line is plotted.
1323 """
1324 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1325 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1326 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1327 # over the pairs simultaneously.
1329 # Ensure cube_x and cube_y are iterable
1330 cube_x_iterable = iter_maybe(cube_x)
1331 cube_y_iterable = iter_maybe(cube_y)
1333 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1334 iplt.scatter(cube_x_iter, cube_y_iter)
1335 if one_to_one is True:
1336 plt.plot(
1337 [
1338 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1339 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1340 ],
1341 [
1342 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1343 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1344 ],
1345 "k",
1346 linestyle="--",
1347 )
1348 ax = plt.gca()
1350 # Add some labels and tweak the style.
1351 if model_names is None:
1352 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1353 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1354 else:
1355 # Add the model names, these should be order of base (x) and other (y).
1356 ax.set_xlabel(
1357 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1358 )
1359 ax.set_ylabel(
1360 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1361 )
1362 ax.set_title(title, fontsize=16)
1363 ax.ticklabel_format(axis="y", useOffset=False)
1364 ax.tick_params(axis="x", labelrotation=15)
1365 ax.tick_params(axis="both", labelsize=12)
1366 ax.autoscale()
1368 # Save plot.
1369 if not in_sphinx_gallery(): 1369 ↛ exitline 1369 didn't return from function '_plot_and_save_scatter_plot' because the condition on line 1369 was always true
1370 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1371 logging.info("Saved scatter plot to %s", filename)
1372 plt.close(fig)
1375def _plot_and_save_vector_plot(
1376 cube_u: iris.cube.Cube,
1377 cube_v: iris.cube.Cube,
1378 filename: str,
1379 title: str,
1380 method: Literal["contourf", "pcolormesh"],
1381 **kwargs,
1382):
1383 """Plot and save a 2D vector plot.
1385 Parameters
1386 ----------
1387 cube_u: Cube
1388 2 dimensional Cube of u component of the data.
1389 cube_v: Cube
1390 2 dimensional Cube of v component of the data.
1391 filename: str
1392 Filename of the plot to write.
1393 title: str
1394 Plot title.
1395 """
1396 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1397 # Create a cube containing the magnitude of the vector field.
1398 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1399 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1400 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1401 cube_vec_mag.rename(
1402 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1403 )
1405 # Specify the color bar
1406 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1408 # Setup plot map projection, extent and coastlines and borderlines.
1409 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1411 if method == "contourf":
1412 # Filled contour plot of the field.
1413 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1414 elif method == "pcolormesh":
1415 try:
1416 vmin = min(levels)
1417 vmax = max(levels)
1418 except TypeError:
1419 vmin, vmax = None, None
1420 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1421 # if levels are defined.
1422 if norm is not None:
1423 vmin = None
1424 vmax = None
1425 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1426 else:
1427 raise ValueError(f"Unknown plotting method: {method}")
1429 # Check to see if transect, and if so, adjust y axis.
1430 if is_transect(cube_vec_mag):
1431 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1432 axes.invert_yaxis()
1433 axes.set_yscale("log")
1434 axes.set_ylim(1100, 100)
1435 # If both model_level_number and level_height exists, iplt can construct
1436 # plot as a function of height above orography (NOT sea level).
1437 elif {"model_level_number", "level_height"}.issubset(
1438 {coord.name() for coord in cube_vec_mag.coords()}
1439 ):
1440 axes.set_yscale("log")
1442 axes.set_title(
1443 f"{title}\n"
1444 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1445 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1446 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1447 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1448 fontsize=16,
1449 )
1451 else:
1452 # Add title.
1453 axes.set_title(title, fontsize=16)
1455 # Add watermark with min/max/mean. Currently not user togglable.
1456 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1457 axes.annotate(
1458 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}",
1459 xy=(0.05, -0.05),
1460 xycoords="axes fraction",
1461 xytext=(-5, 5),
1462 textcoords="offset points",
1463 ha="right",
1464 va="bottom",
1465 size=11,
1466 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1467 )
1469 # Add colour bar.
1470 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1471 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1472 # add ticks and tick_labels for every levels if less than 20 levels exist
1473 if levels is not None and len(levels) < 20:
1474 cbar.set_ticks(levels)
1475 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1477 # 30 barbs along the longest axis of the plot, or a barb per point for data
1478 # with less than 30 points.
1479 step = max(max(cube_u.shape) // 30, 1)
1480 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1482 # Save plot.
1483 if not in_sphinx_gallery():
1484 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1485 logging.info("Saved vector plot to %s", filename)
1486 plt.close(fig)
1489def _plot_and_save_histogram_series(
1490 cubes: iris.cube.Cube | iris.cube.CubeList,
1491 filename: str,
1492 title: str,
1493 vmin: float,
1494 vmax: float,
1495 **kwargs,
1496):
1497 """Plot and save a histogram series.
1499 Parameters
1500 ----------
1501 cubes: Cube or CubeList
1502 2 dimensional Cube or CubeList of the data to plot as histogram.
1503 filename: str
1504 Filename of the plot to write.
1505 title: str
1506 Plot title.
1507 vmin: float
1508 minimum for colorbar
1509 vmax: float
1510 maximum for colorbar
1511 """
1512 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1513 ax = plt.gca()
1515 model_colors_map = get_model_colors_map(cubes)
1517 # Set default that histograms will produce probability density function
1518 # at each bin (integral over range sums to 1).
1519 density = True
1521 for cube in iter_maybe(cubes):
1522 # Easier to check title (where var name originates)
1523 # than seeing if long names exist etc.
1524 # Exception case, where distribution better fits log scales/bins.
1525 if (
1526 ("surface_microphysical" in title)
1527 or ("rain accumulation" in title)
1528 or ("Rainfall rate Composite" in title)
1529 or ("Nimrod_5min" in title)
1530 ):
1531 if "amount" in title: 1531 ↛ 1533line 1531 didn't jump to line 1533 because the condition on line 1531 was never true
1532 # Compute histogram following Klingaman et al. (2017): ASoP
1533 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1534 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1535 density = False
1536 else:
1537 bins = 10.0 ** (
1538 np.arange(-10, 27, 1) / 10.0
1539 ) # Suggestion from RMED toolbox.
1540 bins = np.insert(bins, 0, 0)
1541 ax.set_yscale("log")
1542 vmin = bins[1]
1543 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1544 ax.set_xscale("log")
1545 elif "lightning" in title:
1546 bins = [0, 1, 2, 3, 4, 5]
1547 else:
1548 bins = np.linspace(vmin, vmax, 51)
1549 logging.debug(
1550 "Plotting histogram with %s bins %s - %s.",
1551 np.size(bins),
1552 np.min(bins),
1553 np.max(bins),
1554 )
1556 # Reshape cube data into a single array to allow for a single histogram.
1557 # Otherwise we plot xdim histograms stacked.
1558 cube_data_1d = (cube.data).flatten()
1560 label = None
1561 color = "black"
1562 if model_colors_map:
1563 label = cube.attributes.get("model_name")
1564 color = model_colors_map[label]
1565 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1567 # Compute area under curve.
1568 if ( 1568 ↛ 1574line 1568 didn't jump to line 1574 because the condition on line 1568 was never true
1569 ("surface_microphysical" in title and "amount" in title)
1570 or ("rain_accumulation" in title)
1571 or ("Rainfall rate Composite" in title)
1572 or ("Nimrod_5min" in title)
1573 ):
1574 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1575 x = x * bin_mean / x.sum()
1576 x = x[1:]
1577 y = y[1:]
1579 ax.plot(
1580 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1581 )
1583 # Add some labels and tweak the style.
1584 ax.set_title(title, fontsize=16)
1585 ax.set_xlabel(
1586 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1587 )
1588 ax.set_ylabel("Normalised probability density", fontsize=14)
1589 if ( 1589 ↛ 1594line 1589 didn't jump to line 1594 because the condition on line 1589 was never true
1590 ("surface_microphysical" in title and "amount" in title)
1591 or ("rain accumulation" in title)
1592 or ("Nimrod_5min" in title)
1593 ):
1594 ax.set_ylabel(
1595 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1596 )
1597 ax.set_xlim(vmin, vmax)
1598 ax.tick_params(axis="both", labelsize=12)
1600 # Overlay grid-lines onto histogram plot.
1601 ax.grid(linestyle="--", color="grey", linewidth=1)
1602 if model_colors_map:
1603 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1605 # Save plot.
1606 if not in_sphinx_gallery(): 1606 ↛ exitline 1606 didn't return from function '_plot_and_save_histogram_series' because the condition on line 1606 was always true
1607 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1608 logging.info("Saved histogram plot to %s", filename)
1609 plt.close(fig)
1612def _plot_and_save_postage_stamp_histogram_series(
1613 cube: iris.cube.Cube,
1614 filename: str,
1615 title: str,
1616 stamp_coordinate: str,
1617 vmin: float,
1618 vmax: float,
1619 **kwargs,
1620):
1621 """Plot and save postage (ensemble members) stamps for a histogram series.
1623 Parameters
1624 ----------
1625 cube: Cube
1626 2 dimensional Cube of the data to plot as histogram.
1627 filename: str
1628 Filename of the plot to write.
1629 title: str
1630 Plot title.
1631 stamp_coordinate: str
1632 Coordinate that becomes different plots.
1633 vmin: float
1634 minimum for pdf x-axis
1635 vmax: float
1636 maximum for pdf x-axis
1637 """
1638 # Use the smallest square grid that will fit the members.
1639 nmember = len(cube.coord(stamp_coordinate).points)
1640 grid_rows = int(math.sqrt(nmember))
1641 grid_size = math.ceil(nmember / grid_rows)
1643 fig = plt.figure(
1644 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1645 )
1646 # Make a subplot for each member.
1647 for member, subplot in zip(
1648 cube.slices_over(stamp_coordinate),
1649 range(1, grid_size * grid_rows + 1),
1650 strict=False,
1651 ):
1652 # Implicit interface is much easier here, due to needing to have the
1653 # cartopy GeoAxes generated.
1654 plt.subplot(grid_rows, grid_size, subplot)
1655 # Reshape cube data into a single array to allow for a single histogram.
1656 # Otherwise we plot xdim histograms stacked.
1657 member_data_1d = (member.data).flatten()
1658 plt.hist(member_data_1d, density=True, stacked=True)
1659 axes = plt.gca()
1660 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1661 axes.set_title(f"{mtitle}")
1662 axes.set_xlim(vmin, vmax)
1664 # Overall figure title.
1665 fig.suptitle(title, fontsize=16)
1666 if not in_sphinx_gallery(): 1666 ↛ exitline 1666 didn't return from function '_plot_and_save_postage_stamp_histogram_series' because the condition on line 1666 was always true
1667 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1668 logging.info("Saved histogram postage stamp plot to %s", filename)
1669 plt.close(fig)
1672def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1673 cube: iris.cube.Cube,
1674 filename: str,
1675 title: str,
1676 stamp_coordinate: str,
1677 vmin: float,
1678 vmax: float,
1679 **kwargs,
1680):
1681 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1682 ax.set_title(title, fontsize=16)
1683 ax.set_xlim(vmin, vmax)
1684 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1685 ax.set_ylabel("normalised probability density", fontsize=14)
1686 # Loop over all slices along the stamp_coordinate
1687 for member in cube.slices_over(stamp_coordinate):
1688 # Flatten the member data to 1D
1689 member_data_1d = member.data.flatten()
1690 # Plot the histogram using plt.hist
1691 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1692 plt.hist(
1693 member_data_1d,
1694 density=True,
1695 stacked=True,
1696 label=f"{mtitle}",
1697 )
1699 # Add a legend
1700 ax.legend(fontsize=16)
1702 # Save the figure to a file
1703 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1704 logging.info("Saved histogram postage stamp plot to %s", filename)
1706 # Close the figure
1707 plt.close(fig)
1710def _spatial_plot(
1711 method: Literal["contourf", "pcolormesh", "scatter"],
1712 cube: iris.cube.Cube,
1713 filename: str | None,
1714 sequence_coordinate: str,
1715 stamp_coordinate: str,
1716 overlay_cube: iris.cube.Cube | None = None,
1717 contour_cube: iris.cube.Cube | None = None,
1718 point_cube: iris.cube.Cube | None = None,
1719 **kwargs,
1720):
1721 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1723 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1724 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1725 is present then postage stamp plots will be produced.
1727 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1728 be overplotted on the same figure.
1730 Parameters
1731 ----------
1732 method: "contourf" | "pcolormesh" | "scatter"
1733 The plotting method to use.
1734 Select choice of "contourf" or "pcolormesh" for gridded data.
1735 Use "scatter" for point-based data.
1736 cube: Cube
1737 Iris cube of the data to plot. It should have two spatial dimensions,
1738 such as lat and lon, and may also have a another two dimension to be
1739 plotted sequentially and/or as postage stamp plots.
1740 filename: str | None
1741 Name of the plot to write, used as a prefix for plot sequences. If None
1742 uses the recipe name.
1743 sequence_coordinate: str
1744 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1745 This coordinate must exist in the cube.
1746 stamp_coordinate: str
1747 Coordinate about which to plot postage stamp plots. Defaults to
1748 ``"realization"``.
1749 overlay_cube: Cube | None, optional
1750 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1751 contour_cube: Cube | None, optional
1752 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1753 point_cube: Cube | None, optional
1754 Optional 1 dimensional (e.g. list of points) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
1756 Raises
1757 ------
1758 ValueError
1759 If the cube doesn't have the right dimensions.
1760 TypeError
1761 If the cube isn't a single cube.
1762 """
1763 # Ensure we've got a single cube.
1764 cube = check_single_cube(cube)
1766 # Set title based on recipe metadata or use cube name
1767 recipe_title = get_recipe_metadata().get("title", cube.name())
1769 # Check if there is a valid stamp coordinate in cube dimensions.
1770 if stamp_coordinate == "realization": 1770 ↛ 1775line 1770 didn't jump to line 1775 because the condition on line 1770 was always true
1771 stamp_coordinate = check_stamp_coordinate(cube)
1773 # Make postage stamp plots if stamp_coordinate exists and has more than a
1774 # single point.
1775 plotting_func = _plot_and_save_spatial_plot
1776 try:
1777 if cube.coord(stamp_coordinate).shape[0] > 1:
1778 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1779 except iris.exceptions.CoordinateNotFoundError:
1780 pass
1782 # Produce a geographical scatter plot if the data have a
1783 # dimension called observation or model_obs_error
1784 if any( 1784 ↛ 1788line 1784 didn't jump to line 1788 because the condition on line 1784 was never true
1785 crd.var_name == "station" or crd.var_name == "model_obs_error"
1786 for crd in cube.coords()
1787 ):
1788 plotting_func = _plot_and_save_spatial_plot
1789 method = "scatter"
1791 # Must have a sequence coordinate.
1792 try:
1793 cube.coord(sequence_coordinate)
1794 except iris.exceptions.CoordinateNotFoundError as err:
1795 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1797 # Create a plot for each value of the sequence coordinate.
1798 plot_index = []
1799 nplot = np.size(cube.coord(sequence_coordinate).points)
1801 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1802 # Set plot titles and filename
1803 seq_coord = cube_slice.coord(sequence_coordinate)
1804 plot_title, plot_filename = _set_title_and_filename(
1805 seq_coord, nplot, recipe_title, filename
1806 )
1808 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1809 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1810 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1811 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1813 # Do the actual plotting.
1814 plotting_func(
1815 cube_slice,
1816 filename=plot_filename,
1817 stamp_coordinate=stamp_coordinate,
1818 title=plot_title,
1819 method=method,
1820 overlay_cube=overlay_slice,
1821 contour_cube=contour_slice,
1822 point_cube=point_slice,
1823 **kwargs,
1824 )
1825 plot_index.append(plot_filename)
1827 # Add list of plots to plot metadata.
1828 complete_plot_index = _append_to_plot_index(plot_index)
1830 # Make a page to display the plots.
1831 _make_plot_html_page(complete_plot_index)
1834####################
1835# Public functions #
1836####################
1839def spatial_contour_plot(
1840 cube: iris.cube.Cube,
1841 filename: str = None,
1842 sequence_coordinate: str = "time",
1843 stamp_coordinate: str = "realization",
1844 **kwargs,
1845) -> iris.cube.Cube:
1846 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1848 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1849 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1850 is present then postage stamp plots will be produced.
1852 Parameters
1853 ----------
1854 cube: Cube
1855 Iris cube of the data to plot. It should have two spatial dimensions,
1856 such as lat and lon, and may also have a another two dimension to be
1857 plotted sequentially and/or as postage stamp plots.
1858 filename: str, optional
1859 Name of the plot to write, used as a prefix for plot sequences. Defaults
1860 to the recipe name.
1861 sequence_coordinate: str, optional
1862 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1863 This coordinate must exist in the cube.
1864 stamp_coordinate: str, optional
1865 Coordinate about which to plot postage stamp plots. Defaults to
1866 ``"realization"``.
1868 Returns
1869 -------
1870 Cube
1871 The original cube (so further operations can be applied).
1873 Raises
1874 ------
1875 ValueError
1876 If the cube doesn't have the right dimensions.
1877 TypeError
1878 If the cube isn't a single cube.
1879 """
1880 _spatial_plot(
1881 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1882 )
1883 return cube
1886def spatial_pcolormesh_plot(
1887 cube: iris.cube.Cube,
1888 filename: str = None,
1889 sequence_coordinate: str = "time",
1890 stamp_coordinate: str = "realization",
1891 **kwargs,
1892) -> iris.cube.Cube:
1893 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1895 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1896 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1897 is present then postage stamp plots will be produced.
1899 This function is significantly faster than ``spatial_contour_plot``,
1900 especially at high resolutions, and should be preferred unless contiguous
1901 contour areas are important.
1903 Parameters
1904 ----------
1905 cube: Cube
1906 Iris cube of the data to plot. It should have two spatial dimensions,
1907 such as lat and lon, and may also have a another two dimension to be
1908 plotted sequentially and/or as postage stamp plots.
1909 filename: str, optional
1910 Name of the plot to write, used as a prefix for plot sequences. Defaults
1911 to the recipe name.
1912 sequence_coordinate: str, optional
1913 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1914 This coordinate must exist in the cube.
1915 stamp_coordinate: str, optional
1916 Coordinate about which to plot postage stamp plots. Defaults to
1917 ``"realization"``.
1919 Returns
1920 -------
1921 Cube
1922 The original cube (so further operations can be applied).
1924 Raises
1925 ------
1926 ValueError
1927 If the cube doesn't have the right dimensions.
1928 TypeError
1929 If the cube isn't a single cube.
1930 """
1931 _spatial_plot(
1932 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1933 )
1934 return cube
1937def spatial_multi_pcolormesh_plot(
1938 cube: iris.cube.Cube,
1939 overlay_cube: iris.cube.Cube | None = None,
1940 contour_cube: iris.cube.Cube | None = None,
1941 point_cube: iris.cube.Cube | None = None,
1942 filename: str | None = None,
1943 sequence_coordinate: str = "time",
1944 stamp_coordinate: str = "realization",
1945 **kwargs,
1946) -> iris.cube.Cube:
1947 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1949 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1950 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1951 is present then postage stamp plots will be produced.
1953 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1955 If specified, contours of a contour_cube can be overplotted on top of those.
1957 If specified, a spatial scatter map of point_cube can be overplotted.
1959 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1961 This function is significantly faster than ``spatial_contour_plot``,
1962 especially at high resolutions, and should be preferred unless contiguous
1963 contour areas are important.
1965 Parameters
1966 ----------
1967 cube: Cube
1968 Iris cube of the data to plot. It should have two spatial dimensions,
1969 such as lat and lon, and may also have two additional dimensions to be
1970 plotted sequentially and/or as postage stamp plots.
1971 overlay_cube: Cube, optional
1972 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1973 such as lat and lon, and may also have two additional dimensions to be
1974 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.
1975 If not provided, output plot generated without overlay cube.
1976 contour_cube: Cube, optional
1977 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,
1978 such as lat and lon, and may also have two additional dimensions to be
1979 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
1980 point_cube: Cube, optional
1981 Iris cube of the data to plot as a scatter map overlay on top of basis cube (overlay_cube and/or contour_cube). It should have two
1982 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
1983 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
1984 postage stamp plots. If not provided, output plot generated without point-based layer.
1985 filename: str, optional
1986 Name of the plot to write, used as a prefix for plot sequences. Defaults
1987 to the recipe name.
1988 sequence_coordinate: str, optional
1989 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1990 This coordinate must exist in the cube.
1991 stamp_coordinate: str, optional
1992 Coordinate about which to plot postage stamp plots. Defaults to
1993 ``"realization"``.
1995 Returns
1996 -------
1997 Cube
1998 The original cube (so further operations can be applied).
2000 Raises
2001 ------
2002 ValueError
2003 If the cube doesn't have the right dimensions.
2004 TypeError
2005 If the cube isn't a single cube.
2006 """
2007 _spatial_plot(
2008 "pcolormesh",
2009 cube,
2010 filename,
2011 sequence_coordinate,
2012 stamp_coordinate,
2013 overlay_cube=overlay_cube,
2014 contour_cube=contour_cube,
2015 point_cube=point_cube,
2016 )
2017 return cube, overlay_cube, contour_cube, point_cube
2020# TODO: Expand function to handle ensemble data.
2021# line_coordinate: str, optional
2022# Coordinate about which to plot multiple lines. Defaults to
2023# ``"realization"``.
2024def plot_line_series(
2025 cube: iris.cube.Cube | iris.cube.CubeList,
2026 filename: str = None,
2027 series_coordinate: str = "time",
2028 sequence_coordinate: str = "time",
2029 # add the following for ensembles
2030 stamp_coordinate: str = "realization",
2031 single_plot: bool = False,
2032 **kwargs,
2033) -> iris.cube.Cube | iris.cube.CubeList:
2034 """Plot a line plot for the specified coordinate.
2036 The Cube or CubeList must be 1D.
2038 Parameters
2039 ----------
2040 iris.cube | iris.cube.CubeList
2041 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2042 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2043 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2044 filename: str, optional
2045 Name of the plot to write, used as a prefix for plot sequences. Defaults
2046 to the recipe name.
2047 series_coordinate: str, optional
2048 Coordinate about which to make a series. Defaults to ``"time"``. This
2049 coordinate must exist in the cube.
2051 Returns
2052 -------
2053 iris.cube.Cube | iris.cube.CubeList
2054 The original Cube or CubeList (so further operations can be applied).
2056 Raises
2057 ------
2058 ValueError
2059 If the cubes don't have the right dimensions.
2060 TypeError
2061 If the cube isn't a Cube or CubeList.
2062 """
2063 # Ensure we have a name for the plot file.
2064 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2066 num_models = get_num_models(cube)
2068 validate_cube_shape(cube, num_models)
2070 # Iterate over all cubes and extract coordinate to plot.
2071 cubes = iter_maybe(cube)
2072 coords = []
2073 for cube in cubes:
2074 try:
2075 coords.append(cube.coord(series_coordinate))
2076 except iris.exceptions.CoordinateNotFoundError as err:
2077 raise ValueError(
2078 f"Cube must have a {series_coordinate} coordinate."
2079 ) from err
2080 if cube.coords("realization"): 2080 ↛ 2084line 2080 didn't jump to line 2084 because the condition on line 2080 was always true
2081 if cube.ndim > 3: 2081 ↛ 2082line 2081 didn't jump to line 2082 because the condition on line 2081 was never true
2082 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2083 else:
2084 raise ValueError("Cube must have a realization coordinate.")
2086 plot_index = []
2088 # Check if this is a spectral plot by looking for spectral coordinates
2089 is_spectral_plot = series_coordinate in [
2090 "frequency",
2091 "physical_wavenumber",
2092 "wavelength",
2093 ]
2095 if is_spectral_plot:
2096 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2097 # coordinate frequency/wavenumber.
2098 # If several power spectra are plotted with time as sequence_coordinate for the
2099 # time slider option.
2101 # Internal plotting function.
2102 plotting_func = _plot_and_save_line_power_spectrum_series
2104 for cube in cubes:
2105 try:
2106 cube.coord(sequence_coordinate)
2107 except iris.exceptions.CoordinateNotFoundError as err:
2108 raise ValueError(
2109 f"Cube must have a {sequence_coordinate} coordinate."
2110 ) from err
2112 if num_models == 1: 2112 ↛ 2126line 2112 didn't jump to line 2126 because the condition on line 2112 was always true
2113 # check for ensembles
2114 if ( 2114 ↛ 2118line 2114 didn't jump to line 2118 because the condition on line 2114 was never true
2115 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2116 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2117 ):
2118 if single_plot:
2119 # Plot spectra, mean and ensemble spread on 1 plot
2120 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2121 else:
2122 # Plot postage stamps
2123 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2124 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2125 else:
2126 all_points = sorted(
2127 set(
2128 itertools.chain.from_iterable(
2129 cb.coord(sequence_coordinate).points for cb in cubes
2130 )
2131 )
2132 )
2133 all_slices = list(
2134 itertools.chain.from_iterable(
2135 cb.slices_over(sequence_coordinate) for cb in cubes
2136 )
2137 )
2138 # Matched slices (matched by seq coord point; it may happen that
2139 # evaluated models do not cover the same seq coord range, hence matching
2140 # necessary)
2141 cube_iterables = [
2142 iris.cube.CubeList(
2143 s
2144 for s in all_slices
2145 if s.coord(sequence_coordinate).points[0] == point
2146 )
2147 for point in all_points
2148 ]
2150 nplot = np.size(cube.coord(sequence_coordinate).points)
2152 # Create a plot for each value of the sequence coordinate. Allowing for
2153 # multiple cubes in a CubeList to be plotted in the same plot for similar
2154 # sequence values. Passing a CubeList into the internal plotting function
2155 # for similar values of the sequence coordinate. cube_slice can be an
2156 # iris.cube.Cube or an iris.cube.CubeList.
2158 for cube_slice in cube_iterables:
2159 # Normalize cube_slice to a list of cubes
2160 if isinstance(cube_slice, iris.cube.CubeList): 2160 ↛ 2161line 2160 didn't jump to line 2161 because the condition on line 2160 was never true
2161 cubes = list(cube_slice)
2162 elif isinstance(cube_slice, iris.cube.Cube): 2162 ↛ 2165line 2162 didn't jump to line 2165 because the condition on line 2162 was always true
2163 cubes = [cube_slice]
2164 else:
2165 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2167 # Use sequence value so multiple sequences can merge.
2168 seq_coord = cube_slice[0].coord(sequence_coordinate)
2169 plot_title, plot_filename = _set_title_and_filename(
2170 seq_coord, nplot, recipe_title, filename
2171 )
2173 # Format the coordinate value in a unit appropriate way.
2174 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2176 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2177 if nplot == 1 and seq_coord.has_bounds: 2177 ↛ 2182line 2177 didn't jump to line 2182 because the condition on line 2177 was always true
2178 if np.size(seq_coord.bounds) > 1: 2178 ↛ 2179line 2178 didn't jump to line 2179 because the condition on line 2178 was never true
2179 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2181 # Do the actual plotting.
2182 plotting_func(
2183 cube_slice,
2184 coords,
2185 stamp_coordinate,
2186 plot_filename,
2187 title,
2188 series_coordinate,
2189 )
2191 plot_index.append(plot_filename)
2192 else:
2193 # Format the title and filename using plotted series coordinate
2194 nplot = 1
2195 seq_coord = coords[0]
2196 plot_title, plot_filename = _set_title_and_filename(
2197 seq_coord, nplot, recipe_title, filename
2198 )
2199 # Do the actual plotting for all other series coordinate options.
2200 _plot_and_save_line_series(
2201 cubes, coords, stamp_coordinate, plot_filename, plot_title
2202 )
2204 plot_index.append(plot_filename)
2206 # append plot to list of plots
2207 complete_plot_index = _append_to_plot_index(plot_index)
2209 # Make a page to display the plots.
2210 _make_plot_html_page(complete_plot_index)
2212 return cube
2215def plot_vertical_line_series(
2216 cubes: iris.cube.Cube | iris.cube.CubeList,
2217 filename: str = None,
2218 series_coordinate: str = "model_level_number",
2219 sequence_coordinate: str = "time",
2220 # line_coordinate: str = "realization",
2221 **kwargs,
2222) -> iris.cube.Cube | iris.cube.CubeList:
2223 """Plot a line plot against a type of vertical coordinate.
2225 The Cube or CubeList must be 1D.
2227 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2228 then a sequence of plots will be produced.
2230 Parameters
2231 ----------
2232 iris.cube | iris.cube.CubeList
2233 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2234 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2235 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2236 filename: str, optional
2237 Name of the plot to write, used as a prefix for plot sequences. Defaults
2238 to the recipe name.
2239 series_coordinate: str, optional
2240 Coordinate to plot on the y-axis. Can be ``pressure`` or
2241 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2242 for LFRic. Defaults to ``model_level_number``.
2243 This coordinate must exist in the cube.
2244 sequence_coordinate: str, optional
2245 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2246 This coordinate must exist in the cube.
2248 Returns
2249 -------
2250 iris.cube.Cube | iris.cube.CubeList
2251 The original Cube or CubeList (so further operations can be applied).
2252 Plotted data.
2254 Raises
2255 ------
2256 ValueError
2257 If the cubes doesn't have the right dimensions.
2258 TypeError
2259 If the cube isn't a Cube or CubeList.
2260 """
2261 # Ensure we have a name for the plot file.
2262 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2264 cubes = iter_maybe(cubes)
2265 # Initialise empty list to hold all data from all cubes in a CubeList
2266 all_data = []
2268 # Store min/max ranges for x range.
2269 x_levels = []
2271 num_models = get_num_models(cubes)
2273 validate_cube_shape(cubes, num_models)
2275 # Iterate over all cubes in cube or CubeList and plot.
2276 coords = []
2277 for cube in cubes:
2278 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2279 try:
2280 coords.append(cube.coord(series_coordinate))
2281 except iris.exceptions.CoordinateNotFoundError as err:
2282 raise ValueError(
2283 f"Cube must have a {series_coordinate} coordinate."
2284 ) from err
2286 try:
2287 if cube.ndim > 1 or not cube.coords("realization"): 2287 ↛ 2295line 2287 didn't jump to line 2295 because the condition on line 2287 was always true
2288 cube.coord(sequence_coordinate)
2289 except iris.exceptions.CoordinateNotFoundError as err:
2290 raise ValueError(
2291 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2292 ) from err
2294 # Get minimum and maximum from levels information.
2295 _, levels, _ = colorbar_map_levels(cube, axis="x")
2296 if levels is not None: 2296 ↛ 2300line 2296 didn't jump to line 2300 because the condition on line 2296 was always true
2297 x_levels.append(min(levels))
2298 x_levels.append(max(levels))
2299 else:
2300 all_data.append(cube.data)
2302 if len(x_levels) == 0: 2302 ↛ 2304line 2302 didn't jump to line 2304 because the condition on line 2302 was never true
2303 # Combine all data into a single NumPy array
2304 combined_data = np.concatenate(all_data)
2306 # Set the lower and upper limit for the x-axis to ensure all plots have
2307 # same range. This needs to read the whole cube over the range of the
2308 # sequence and if applicable postage stamp coordinate.
2309 vmin = np.floor(combined_data.min())
2310 vmax = np.ceil(combined_data.max())
2311 else:
2312 vmin = min(x_levels)
2313 vmax = max(x_levels)
2315 # Matching the slices (matching by seq coord point; it may happen that
2316 # evaluated models do not cover the same seq coord range, hence matching
2317 # necessary)
2318 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2320 # Create a plot for each value of the sequence coordinate.
2321 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2322 # similar sequence values. Passing a CubeList into the internal plotting function
2323 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2324 # or an iris.cube.CubeList.
2325 plot_index = []
2326 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2327 for cubes_slice in cube_iterables:
2328 # Format the coordinate value in a unit appropriate way.
2329 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2330 plot_title, plot_filename = _set_title_and_filename(
2331 seq_coord, nplot, recipe_title, filename
2332 )
2334 # Do the actual plotting.
2335 _plot_and_save_vertical_line_series(
2336 cubes_slice,
2337 coords,
2338 "realization",
2339 plot_filename,
2340 series_coordinate,
2341 title=plot_title,
2342 vmin=vmin,
2343 vmax=vmax,
2344 )
2345 plot_index.append(plot_filename)
2347 # Add list of plots to plot metadata.
2348 complete_plot_index = _append_to_plot_index(plot_index)
2350 # Make a page to display the plots.
2351 _make_plot_html_page(complete_plot_index)
2353 return cubes
2356def qq_plot(
2357 cubes: iris.cube.CubeList,
2358 coordinates: list[str],
2359 percentiles: list[float],
2360 model_names: list[str],
2361 filename: str = None,
2362 one_to_one: bool = True,
2363 **kwargs,
2364) -> iris.cube.CubeList:
2365 """Plot a Quantile-Quantile plot between two models for common time points.
2367 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2368 collapsed within the operator over all specified coordinates such as
2369 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2370 ensemble members to ensure a 1D cube (array).
2372 Parameters
2373 ----------
2374 cubes: iris.cube.CubeList
2375 Two cubes of the same variable with different models.
2376 coordinate: list[str]
2377 The list of coordinates to collapse over. This list should be
2378 every coordinate within the cube to result in a 1D cube around
2379 the percentile coordinate.
2380 percent: list[float]
2381 A list of percentiles to appear in the plot.
2382 model_names: list[str]
2383 A list of model names to appear on the axis of the plot.
2384 filename: str, optional
2385 Filename of the plot to write.
2386 one_to_one: bool, optional
2387 If True a 1:1 line is plotted; if False it is not. Default is True.
2389 Raises
2390 ------
2391 ValueError
2392 When the cubes are not compatible.
2394 Notes
2395 -----
2396 The quantile-quantile plot is a variant on the scatter plot representing
2397 two datasets by their quantiles (percentiles) for common time points.
2398 This plot does not use a theoretical distribution to compare against, but
2399 compares percentiles of two datasets. This plot does
2400 not use all raw data points, but plots the selected percentiles (quantiles) of
2401 each variable instead for the two datasets, thereby normalising the data for a
2402 direct comparison between the selected percentiles of the two dataset distributions.
2404 Quantile-quantile plots are valuable for comparing against
2405 observations and other models. Identical percentiles between the variables
2406 will lie on the one-to-one line implying the values correspond well to each
2407 other. Where there is a deviation from the one-to-one line a range of
2408 possibilities exist depending on how and where the data is shifted (e.g.,
2409 Wilks 2011 [Wilks2011]_).
2411 For distributions above the one-to-one line the distribution is left-skewed;
2412 below is right-skewed. A distinct break implies a bimodal distribution, and
2413 closer values/values further apart at the tails imply poor representation of
2414 the extremes.
2416 References
2417 ----------
2418 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2419 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2420 """
2421 # Check cubes using same functionality as the difference operator.
2422 if len(cubes) != 2:
2423 raise ValueError("cubes should contain exactly 2 cubes.")
2424 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2425 other: Cube = cubes.extract_cube(
2426 iris.Constraint(
2427 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2428 )
2429 )
2431 # Get spatial coord names.
2432 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2433 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2435 # Ensure cubes to compare are on common differencing grid.
2436 # This is triggered if either
2437 # i) latitude and longitude shapes are not the same. Note grid points
2438 # are not compared directly as these can differ through rounding
2439 # errors.
2440 # ii) or variables are known to often sit on different grid staggering
2441 # in different models (e.g. cell center vs cell edge), as is the case
2442 # for UM and LFRic comparisons.
2443 # In future greater choice of regridding method might be applied depending
2444 # on variable type. Linear regridding can in general be appropriate for smooth
2445 # variables. Care should be taken with interpretation of differences
2446 # given this dependency on regridding.
2447 if (
2448 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2449 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2450 ) or (
2451 base.long_name
2452 in [
2453 "eastward_wind_at_10m",
2454 "northward_wind_at_10m",
2455 "northward_wind_at_cell_centres",
2456 "eastward_wind_at_cell_centres",
2457 "zonal_wind_at_pressure_levels",
2458 "meridional_wind_at_pressure_levels",
2459 "potential_vorticity_at_pressure_levels",
2460 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2461 ]
2462 ):
2463 logging.debug(
2464 "Linear regridding base cube to other grid to compute differences"
2465 )
2466 base = regrid_onto_cube(base, other, method="Linear")
2468 # Extract just common time points.
2469 base, other = _extract_common_time_points(base, other)
2471 # Equalise attributes so we can merge.
2472 fully_equalise_attributes([base, other])
2473 logging.debug("Base: %s\nOther: %s", base, other)
2475 # Collapse cubes.
2476 base = collapse(
2477 base,
2478 coordinate=coordinates,
2479 method="PERCENTILE",
2480 additional_percent=percentiles,
2481 )
2482 other = collapse(
2483 other,
2484 coordinate=coordinates,
2485 method="PERCENTILE",
2486 additional_percent=percentiles,
2487 )
2489 # Ensure we have a name for the plot file.
2490 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2491 title = f"{recipe_title}"
2493 if filename is None:
2494 filename = slugify(recipe_title)
2496 # Add file extension.
2497 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2499 # Do the actual plotting on a scatter plot
2500 _plot_and_save_scatter_plot(
2501 base, other, plot_filename, title, one_to_one, model_names
2502 )
2504 # Add list of plots to plot metadata.
2505 plot_index = _append_to_plot_index([plot_filename])
2507 # Make a page to display the plots.
2508 _make_plot_html_page(plot_index)
2510 return iris.cube.CubeList([base, other])
2513def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2514 """
2515 Plot a Hinton style triangle/scorecard plot.
2517 This plot type can be useful for summarising high level information, such as comparing
2518 how 'skillful' two models are when verified against observations for a variety of metrics,
2519 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2520 than customisable by the user as input arguments; many have been designed to automatically
2521 scale the plot depending on the number of x and y components.
2523 Parameters
2524 ----------
2525 change: np.ndarray
2526 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2527 size/direction.
2528 signif: np.ndarray
2529 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2530 xaxis_labels: list
2531 List of labels for the xaxis (must match the second dimension length of signif and change,
2532 along with magnitude if not None).
2533 yaxis_labels: list
2534 List of labels for the yaxis (must match the first dimension length of signif and change,
2535 along with magnitude if not None).
2536 magnitude: np.ndarray | None
2537 Optional 2D array, matching the shape of change, signif, which contains numerical values
2538 the user wishes to display under each respective triangle.
2540 Returns
2541 -------
2542 matplotlib axes object to either display or do further modifications to.
2543 """
2544 # Setup colors of triangles
2545 color_pos = "#7CAE00"
2546 color_neg = "#7B68EE"
2548 # Setup cell/text size ratios
2549 figsize = None
2550 cell_size_in = 0.35
2551 text_row_ratio = 0.25
2553 # Ensure arrays, and change to bool for sig.
2554 change = np.asarray(change)
2555 signif = np.asarray(signif).astype(bool)
2556 if magnitude is not None: 2556 ↛ 2557line 2556 didn't jump to line 2557 because the condition on line 2556 was never true
2557 magnitude = np.asarray(magnitude)
2559 # Get the number of x and y elements
2560 ny, nx = change.shape
2562 # Build non-uniform y coordinates
2563 tri_height = 1.0
2564 txt_height = text_row_ratio
2566 tri_y = []
2567 txt_y = []
2568 y_edges = [0.0]
2570 y = 0.0
2571 for _j in range(ny):
2572 tri_y.append(y + tri_height / 2)
2573 y += tri_height
2574 y_edges.append(y)
2576 if magnitude is not None: 2576 ↛ 2577line 2576 didn't jump to line 2577 because the condition on line 2576 was never true
2577 txt_y.append(y + txt_height / 2)
2578 y += txt_height
2579 y_edges.append(y)
2581 total_height = y
2583 # Dynamic figure size
2584 if figsize is None: 2584 ↛ 2589line 2584 didn't jump to line 2589 because the condition on line 2584 was always true
2585 width = nx * cell_size_in
2586 height = total_height * cell_size_in + 2
2587 figsize = (width, height)
2589 fig, ax = plt.subplots(figsize=figsize)
2591 # Setup axes and grid.
2592 ax.set_aspect("equal", adjustable="box")
2593 ax.set_xlim(-0.5, nx - 0.5)
2594 ax.set_ylim(0, total_height)
2596 ax.set_xticks(np.arange(nx))
2597 ax.set_xticklabels(xaxis_labels, rotation=90)
2599 ax.set_yticks(tri_y)
2600 ax.set_yticklabels(yaxis_labels)
2602 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2603 ax.set_yticks(y_edges, minor=True)
2605 ax.set_axisbelow(True)
2606 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2607 ax.grid(False, which="major")
2608 ax.tick_params(which="minor", length=0)
2610 ax.invert_yaxis()
2612 # Compute marker scaling (fixed overlap)
2613 fig.canvas.draw()
2615 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2616 width_in, height_in = bbox.width, bbox.height
2618 cell_w = (width_in * fig.dpi) / nx
2619 cell_h = (height_in * fig.dpi) / total_height
2620 cell_pixels = min(cell_w, cell_h)
2622 max_marker_size = (0.6 * cell_pixels) ** 2
2624 text_fontsize = cell_pixels * 0.15
2626 # Plot triangles + text
2627 for j in range(ny):
2628 for i in range(nx):
2629 val = change[j, i]
2630 if np.isnan(val): 2630 ↛ 2631line 2630 didn't jump to line 2631 because the condition on line 2630 was never true
2631 continue
2633 if abs(val) < 0.01: 2633 ↛ 2634line 2633 didn't jump to line 2634 because the condition on line 2633 was never true
2634 continue
2636 sig = signif[j, i]
2637 size = max_marker_size * abs(val)
2639 # Triangle style
2640 if val >= 0:
2641 marker = "^"
2642 color = color_pos
2643 else:
2644 marker = "v"
2645 color = color_neg
2647 if sig:
2648 edgecolor = "black"
2649 linewidth = 0.6
2650 else:
2651 edgecolor = "none"
2652 linewidth = 0.0
2654 # Triangle
2655 ax.scatter(
2656 i,
2657 tri_y[j],
2658 s=size,
2659 marker=marker,
2660 c=color,
2661 edgecolors=edgecolor,
2662 linewidths=linewidth,
2663 zorder=3,
2664 clip_on=True, # ensures no rendering bleed
2665 )
2667 # Text row
2668 if magnitude is not None: 2668 ↛ 2669line 2668 didn't jump to line 2669 because the condition on line 2668 was never true
2669 mag_val = magnitude[j, i]
2671 if not np.isnan(mag_val):
2672 ax.text(
2673 i,
2674 txt_y[j],
2675 f"{mag_val:.1f}",
2676 ha="center",
2677 va="center",
2678 fontsize=text_fontsize,
2679 color="black",
2680 zorder=4,
2681 )
2683 plt.tight_layout()
2684 return fig, ax
2687def scatter_plot(
2688 cube_x: iris.cube.Cube | iris.cube.CubeList,
2689 cube_y: iris.cube.Cube | iris.cube.CubeList,
2690 filename: str = None,
2691 one_to_one: bool = True,
2692 **kwargs,
2693) -> iris.cube.CubeList:
2694 """Plot a scatter plot between two variables.
2696 Both cubes must be 1D.
2698 Parameters
2699 ----------
2700 cube_x: Cube | CubeList
2701 1 dimensional Cube of the data to plot on y-axis.
2702 cube_y: Cube | CubeList
2703 1 dimensional Cube of the data to plot on x-axis.
2704 filename: str, optional
2705 Filename of the plot to write.
2706 one_to_one: bool, optional
2707 If True a 1:1 line is plotted; if False it is not. Default is True.
2709 Returns
2710 -------
2711 cubes: CubeList
2712 CubeList of the original x and y cubes for further processing.
2714 Raises
2715 ------
2716 ValueError
2717 If the cube doesn't have the right dimensions and cubes not the same
2718 size.
2719 TypeError
2720 If the cube isn't a single cube.
2722 Notes
2723 -----
2724 Scatter plots are used for determining if there is a relationship between
2725 two variables. Positive relations have a slope going from bottom left to top
2726 right; Negative relations have a slope going from top left to bottom right.
2727 """
2728 # Iterate over all cubes in cube or CubeList and plot.
2729 for cube_iter in iter_maybe(cube_x):
2730 # Check cubes are correct shape.
2731 cube_iter = check_single_cube(cube_iter)
2732 if cube_iter.ndim > 1:
2733 raise ValueError("cube_x must be 1D.")
2735 # Iterate over all cubes in cube or CubeList and plot.
2736 for cube_iter in iter_maybe(cube_y):
2737 # Check cubes are correct shape.
2738 cube_iter = check_single_cube(cube_iter)
2739 if cube_iter.ndim > 1:
2740 raise ValueError("cube_y must be 1D.")
2742 # Ensure we have a name for the plot file.
2743 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2744 title = f"{recipe_title}"
2746 if filename is None:
2747 filename = slugify(recipe_title)
2749 # Add file extension.
2750 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2752 # Do the actual plotting.
2753 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2755 # Add list of plots to plot metadata.
2756 plot_index = _append_to_plot_index([plot_filename])
2758 # Make a page to display the plots.
2759 _make_plot_html_page(plot_index)
2761 return iris.cube.CubeList([cube_x, cube_y])
2764def vector_plot(
2765 cube_u: iris.cube.Cube,
2766 cube_v: iris.cube.Cube,
2767 filename: str = None,
2768 sequence_coordinate: str = "time",
2769 **kwargs,
2770) -> iris.cube.CubeList:
2771 """Plot a vector plot based on the input u and v components."""
2772 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2774 # Cubes must have a matching sequence coordinate.
2775 try:
2776 # Check that the u and v cubes have the same sequence coordinate.
2777 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2777 ↛ anywhereline 2777 didn't jump anywhere: it always raised an exception.
2778 raise ValueError("Coordinates do not match.")
2779 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2780 raise ValueError(
2781 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2782 ) from err
2784 # Create a plot for each value of the sequence coordinate.
2785 plot_index = []
2786 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2787 for cube_u_slice, cube_v_slice in zip(
2788 cube_u.slices_over(sequence_coordinate),
2789 cube_v.slices_over(sequence_coordinate),
2790 strict=True,
2791 ):
2792 # Format the coordinate value in a unit appropriate way.
2793 seq_coord = cube_u_slice.coord(sequence_coordinate)
2794 plot_title, plot_filename = _set_title_and_filename(
2795 seq_coord, nplot, recipe_title, filename
2796 )
2798 # Do the actual plotting.
2799 _plot_and_save_vector_plot(
2800 cube_u_slice,
2801 cube_v_slice,
2802 filename=plot_filename,
2803 title=plot_title,
2804 method="pcolormesh",
2805 )
2806 plot_index.append(plot_filename)
2808 # Add list of plots to plot metadata.
2809 complete_plot_index = _append_to_plot_index(plot_index)
2811 # Make a page to display the plots.
2812 _make_plot_html_page(complete_plot_index)
2814 return iris.cube.CubeList([cube_u, cube_v])
2817def plot_histogram_series(
2818 cubes: iris.cube.Cube | iris.cube.CubeList,
2819 filename: str = None,
2820 sequence_coordinate: str = "time",
2821 stamp_coordinate: str = "realization",
2822 single_plot: bool = False,
2823 **kwargs,
2824) -> iris.cube.Cube | iris.cube.CubeList:
2825 """Plot a histogram plot for each vertical level provided.
2827 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2828 is present then a sequence of plots will be produced using the time slider
2829 functionality to scroll through histograms against time. If a
2830 stamp_coordinate is present then postage stamp plots will be produced. If
2831 stamp_coordinate and single_plot is True, all postage stamp plots will be
2832 plotted in a single plot instead of separate postage stamp plots.
2834 Parameters
2835 ----------
2836 cubes: Cube | iris.cube.CubeList
2837 Iris cube or CubeList of the data to plot. It should have a single dimension other
2838 than the stamp coordinate.
2839 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2840 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2841 filename: str, optional
2842 Name of the plot to write, used as a prefix for plot sequences. Defaults
2843 to the recipe name.
2844 sequence_coordinate: str, optional
2845 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2846 This coordinate must exist in the cube and will be used for the time
2847 slider.
2848 stamp_coordinate: str, optional
2849 Coordinate about which to plot postage stamp plots. Defaults to
2850 ``"realization"``.
2851 single_plot: bool, optional
2852 If True, all postage stamp plots will be plotted in a single plot. If
2853 False, each postage stamp plot will be plotted separately. Is only valid
2854 if stamp_coordinate exists and has more than a single point.
2856 Returns
2857 -------
2858 iris.cube.Cube | iris.cube.CubeList
2859 The original Cube or CubeList (so further operations can be applied).
2860 Plotted data.
2862 Raises
2863 ------
2864 ValueError
2865 If the cube doesn't have the right dimensions.
2866 TypeError
2867 If the cube isn't a Cube or CubeList.
2868 """
2869 recipe_title = get_recipe_metadata().get("title", "Histogram")
2871 cubes = iter_maybe(cubes)
2872 # Ensure we have a name for the plot file.
2873 if filename is None:
2874 filename = slugify(recipe_title)
2876 # Internal plotting function.
2877 plotting_func = _plot_and_save_histogram_series
2879 num_models = get_num_models(cubes)
2881 validate_cube_shape(cubes, num_models)
2883 # If several histograms are plotted, check sequence_coordinate
2884 check_sequence_coordinate(cubes, sequence_coordinate)
2886 # Get axis minimum and maximum from levels information.
2887 # If no levels set, derive minima and maxima from data in CubeList.
2888 vmin, vmax = _set_axis_range(cubes)
2890 # Make postage stamp plots if stamp_coordinate exists and has more than a
2891 # single point. If single_plot is True:
2892 # -- all postage stamp plots will be plotted in a single plot instead of
2893 # separate postage stamp plots.
2894 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2895 # produced per single model only
2896 if num_models == 1:
2897 if ( 2897 ↛ 2901line 2897 didn't jump to line 2901 because the condition on line 2897 was never true
2898 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2899 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2900 ):
2901 if single_plot:
2902 plotting_func = (
2903 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2904 )
2905 else:
2906 plotting_func = _plot_and_save_postage_stamp_histogram_series
2907 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2908 else:
2909 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2911 plot_index = []
2912 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2913 # Create a plot for each value of the sequence coordinate. Allowing for
2914 # multiple cubes in a CubeList to be plotted in the same plot for similar
2915 # sequence values. Passing a CubeList into the internal plotting function
2916 # for similar values of the sequence coordinate. cube_slice can be an
2917 # iris.cube.Cube or an iris.cube.CubeList.
2918 for cube_slice in cube_iterables:
2919 single_cube = cube_slice
2920 if isinstance(cube_slice, iris.cube.CubeList):
2921 single_cube = cube_slice[0]
2923 # Ensure valid stamp coordinate in cube dimensions
2924 if stamp_coordinate == "realization": 2924 ↛ 2927line 2924 didn't jump to line 2927 because the condition on line 2924 was always true
2925 stamp_coordinate = check_stamp_coordinate(single_cube)
2926 # Set plot titles and filename, based on sequence coordinate
2927 seq_coord = single_cube.coord(sequence_coordinate)
2928 # Use time coordinate in title and filename if single histogram output.
2929 if sequence_coordinate == "realization" and nplot == 1: 2929 ↛ 2930line 2929 didn't jump to line 2930 because the condition on line 2929 was never true
2930 seq_coord = single_cube.coord("time")
2931 plot_title, plot_filename = _set_title_and_filename(
2932 seq_coord, nplot, recipe_title, filename
2933 )
2935 # Do the actual plotting.
2936 plotting_func(
2937 cube_slice,
2938 filename=plot_filename,
2939 stamp_coordinate=stamp_coordinate,
2940 title=plot_title,
2941 vmin=vmin,
2942 vmax=vmax,
2943 )
2944 plot_index.append(plot_filename)
2946 # Add list of plots to plot metadata.
2947 complete_plot_index = _append_to_plot_index(plot_index)
2949 # Make a page to display the plots.
2950 _make_plot_html_page(complete_plot_index)
2952 return cubes
2955def _plot_and_save_postage_stamp_power_spectrum_series(
2956 cubes: iris.cube.Cube,
2957 coords: list[iris.coords.Coord],
2958 stamp_coordinate: str,
2959 filename: str,
2960 title: str,
2961 series_coordinate: str | None = None,
2962 **kwargs,
2963):
2964 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2966 Parameters
2967 ----------
2968 cubes: Cube or CubeList
2969 Cube or Cubelist of the power spectrum data.
2970 coords: list[Coord]
2971 Coordinates to plot on the x-axis, one per cube.
2972 stamp_coordinate: str
2973 Coordinate that becomes different plots.
2974 filename: str
2975 Filename of the plot to write.
2976 title: str
2977 Plot title.
2978 series_coordinate: str, optional
2979 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2981 """
2982 # Use the smallest square grid that will fit the members.
2983 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
2985 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
2986 model_colors_map = get_model_colors_map(cubes)
2987 # ax = plt.gca()
2988 # Make a subplot for each member.
2989 for member, subplot in zip(
2990 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
2991 ):
2992 ax = plt.subplot(grid_size, grid_size, subplot)
2994 # Store min/max ranges.
2995 y_levels = []
2997 line_marker = None
2998 line_width = 1
3000 for cube in iter_maybe(member):
3001 xcoord = _select_series_coord(cube, series_coordinate)
3002 xname = xcoord.points
3004 yfield = cube.data # power spectrum
3005 label = None
3006 color = "black"
3007 if model_colors_map: 3007 ↛ 3008line 3007 didn't jump to line 3008 because the condition on line 3007 was never true
3008 label = cube.attributes.get("model_name")
3009 color = model_colors_map.get(label)
3011 if member.coord(stamp_coordinate).points == [0]:
3012 ax.plot(
3013 xname,
3014 yfield,
3015 color=color,
3016 marker=line_marker,
3017 ls="-",
3018 lw=line_width,
3019 label=f"{label} (control)"
3020 if len(cube.coord(stamp_coordinate).points) > 1
3021 else label,
3022 )
3023 # Label with member if part of an ensemble and not the control.
3024 else:
3025 ax.plot(
3026 xname,
3027 yfield,
3028 color=color,
3029 ls="-",
3030 lw=1.5,
3031 alpha=0.75,
3032 label=f"{label} (member)",
3033 )
3035 # Calculate the global min/max if multiple cubes are given.
3036 _, levels, _ = colorbar_map_levels(cube, axis="y")
3037 if levels is not None: 3037 ↛ 3038line 3037 didn't jump to line 3038 because the condition on line 3037 was never true
3038 y_levels.append(min(levels))
3039 y_levels.append(max(levels))
3041 # Add some labels and tweak the style.
3042 title = f"{title}"
3043 ax.set_title(title, fontsize=16)
3045 # Set appropriate x-axis label based on coordinate
3046 if series_coordinate == "wavelength" or ( 3046 ↛ 3049line 3046 didn't jump to line 3049 because the condition on line 3046 was never true
3047 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3048 ):
3049 ax.set_xlabel("Wavelength (km)", fontsize=14)
3050 elif series_coordinate == "physical_wavenumber" or ( 3050 ↛ 3055line 3050 didn't jump to line 3055 because the condition on line 3050 was always true
3051 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3052 ):
3053 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3054 else: # frequency or check units
3055 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3056 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3057 else:
3058 ax.set_xlabel("Wavenumber", fontsize=14)
3060 ax.set_ylabel("Power Spectral Density", fontsize=14)
3061 ax.tick_params(axis="both", labelsize=12)
3063 # Set log-log scale
3064 ax.set_xscale("log")
3065 ax.set_yscale("log")
3067 # Add gridlines
3068 ax.grid(linestyle="--", color="grey", linewidth=1)
3069 # Ientify unique labels for legend
3070 handles = list(
3071 {
3072 label: handle
3073 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3074 }.values()
3075 )
3076 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3078 ax = plt.gca()
3079 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3081 if not in_sphinx_gallery(): 3081 ↛ exitline 3081 didn't return from function '_plot_and_save_postage_stamp_power_spectrum_series' because the condition on line 3081 was always true
3082 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3083 logging.info("Saved power spectrum histogram plot to %s", filename)
3084 plt.close(fig)
3087def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3088 cubes: iris.cube.Cube,
3089 coords: list[iris.coords.Coord],
3090 stamp_coordinate: str,
3091 filename: str,
3092 title: str,
3093 series_coordinate: str = None,
3094 **kwargs,
3095):
3096 """Plot and save power spectra for ensemble members in single plot.
3098 Parameters
3099 ----------
3100 cubes: Cube or CubeList
3101 Cube or Cubelist of the power spectrum data.
3102 coords: list[Coord]
3103 Coordinates to plot on the x-axis, one per cube.
3104 stamp_coordinate: str
3105 Coordinate that becomes different plots.
3106 filename: str
3107 Filename of the plot to write.
3108 title: str
3109 Plot title.
3110 series_coordinate: str, optional
3111 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3113 """
3114 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3115 model_colors_map = get_model_colors_map(cubes)
3117 line_marker = None
3118 line_width = 1
3120 # Compute ensemble statistics to show spread
3121 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3122 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3123 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3125 xcoord_global = mean_cube.coord(series_coordinate)
3126 x_global = xcoord_global.points
3128 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3129 xcoord = _select_series_coord(member, series_coordinate)
3130 xname = xcoord.points
3132 yfield = member.data # power spectrum
3133 color = "black"
3134 if model_colors_map: 3134 ↛ 3138line 3134 didn't jump to line 3138 because the condition on line 3134 was always true
3135 label = member.attributes.get("model_name") if i == 0 else None
3136 color = model_colors_map.get(label)
3138 if member.coord(stamp_coordinate).points == [0]:
3139 ax.plot(
3140 xname,
3141 yfield,
3142 color=color,
3143 marker=line_marker,
3144 ls="-",
3145 lw=line_width,
3146 label=f"{label} (control)"
3147 if len(member.coord(stamp_coordinate).points) > 1
3148 else label,
3149 )
3150 # Label with member number if part of an ensemble and not the control.
3151 else:
3152 ax.plot(
3153 xname,
3154 yfield,
3155 color=color,
3156 ls="-",
3157 lw=1.5,
3158 alpha=0.75,
3159 label=label,
3160 )
3162 # Set appropriate x-axis label based on coordinate
3163 if series_coordinate == "wavelength" or ( 3163 ↛ 3166line 3163 didn't jump to line 3166 because the condition on line 3163 was never true
3164 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3165 ):
3166 ax.set_xlabel("Wavelength (km)", fontsize=14)
3167 elif series_coordinate == "physical_wavenumber" or ( 3167 ↛ 3172line 3167 didn't jump to line 3172 because the condition on line 3167 was always true
3168 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3169 ):
3170 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3171 else: # frequency or check units
3172 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3173 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3174 else:
3175 ax.set_xlabel("Wavenumber", fontsize=14)
3177 # Add ensemble spread shading
3178 ax.fill_between(
3179 x_global,
3180 min_cube.data,
3181 max_cube.data,
3182 color="grey",
3183 alpha=0.3,
3184 label="Ensemble spread",
3185 )
3187 # Add ensemble mean line
3188 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3190 ax.set_ylabel("Power Spectral Density", fontsize=14)
3191 ax.tick_params(axis="both", labelsize=12)
3193 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3194 # Set log-log scale
3195 ax.set_xscale("log")
3196 ax.set_yscale("log")
3198 # Add gridlines
3199 ax.grid(linestyle="--", color="grey", linewidth=1)
3200 # Identify unique labels for legend
3201 handles = list(
3202 {
3203 label: handle
3204 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3205 }.values()
3206 )
3207 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3209 # Figure title.
3210 ax.set_title(title, fontsize=16)
3212 # Save the figure to a file
3213 if not in_sphinx_gallery(): 3213 ↛ exitline 3213 didn't return from function '_plot_and_save_postage_stamps_in_single_plot_power_spectrum_series' because the condition on line 3213 was always true
3214 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3215 logging.info("Saved power spectrum plot to %s", filename)
3216 plt.close(fig)