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