Coverage for src/CSET/operators/plot.py: 84%
1172 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 12:33 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24import sys
25from typing import Literal
27import cartopy.crs as ccrs
28import cartopy.feature as cfeature
29import iris
30import iris.coords
31import iris.cube
32import iris.exceptions
33import iris.plot as iplt
34import matplotlib as mpl
35import matplotlib.pyplot as plt
36import numpy as np
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._colormaps import (
50 colorbar_map_levels,
51 get_model_colors_map,
52)
53from CSET.operators._utils import (
54 calc_array_stats,
55 check_sequence_coordinate,
56 check_single_cube,
57 check_stamp_coordinate,
58 fully_equalise_attributes,
59 get_cube_yxcoordname,
60 get_num_models,
61 is_transect,
62 slice_over_maybe,
63 validate_cube_shape,
64 validate_cubes_coords,
65)
66from CSET.operators.collapse import collapse
67from CSET.operators.misc import _extract_common_time_points
68from CSET.operators.regrid import regrid_onto_cube
70logger = logging.getLogger(__name__)
72# Use a non-interactive plotting backend.
73mpl.use("agg")
76############################
77# Private helper functions #
78############################
81def in_sphinx_gallery():
82 """Test if running plot code in sphinx-gallery context."""
83 return "sphinx_gallery" in sys.modules
86def _append_to_plot_index(plot_index: list) -> list:
87 """Add plots into the plot index, returning the complete plot index."""
88 with open("meta.json", "r+t", encoding="UTF-8") as fp:
89 fcntl.flock(fp, fcntl.LOCK_EX)
90 fp.seek(0)
91 meta = json.load(fp)
92 complete_plot_index = meta.get("plots", [])
93 complete_plot_index = complete_plot_index + plot_index
94 meta["plots"] = complete_plot_index
95 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
96 os.getenv("DO_CASE_AGGREGATION")
97 ):
98 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
99 fp.seek(0)
100 fp.truncate()
101 json.dump(meta, fp, indent=2)
102 return complete_plot_index
105def _make_plot_html_page(plots: list):
106 """Create a HTML page to display a plot image."""
107 # Debug check that plots actually contains some strings.
108 assert isinstance(plots[0], str)
110 # Load HTML template file.
111 operator_files = importlib.resources.files()
112 template_file = operator_files.joinpath("_plot_page_template.html")
114 # Get some metadata.
115 meta = get_recipe_metadata()
116 title = meta.get("title", "Untitled")
117 description = MarkdownIt().render(meta.get("description", "*No description.*"))
119 # Prepare template variables.
120 variables = {
121 "title": title,
122 "description": description,
123 "initial_plot": plots[0],
124 "plots": plots,
125 "title_slug": slugify(title),
126 }
128 # Render template.
129 html = render_file(template_file, **variables)
131 # Save completed HTML.
132 with open("index.html", "wt", encoding="UTF-8") as fp:
133 fp.write(html)
136def _save_close_figure(figure, plot_type: str, filename: str):
137 """Save generated plot figure file and close figure.
139 If running documentation gallery generation, avoid saving to file.
141 Parameters
142 ----------
143 figure:
144 Matplotlib Figure object holding all plot elements.
145 plot_type: str
146 String identifier for plot type for logging information.
147 filename: str
148 Filename for saved figure.
149 """
150 if not in_sphinx_gallery():
151 figure.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
152 logger.info("Saved %s plot to %s", plot_type, filename)
153 plt.close(figure)
156def _setup_spatial_map(
157 cube: iris.cube.Cube,
158 figure,
159 cmap,
160 grid_size: tuple[int, int] | None = None,
161 subplot: int | None = None,
162):
163 """Define map projections, extent and add coastlines and borderlines for spatial plots.
165 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
166 is specified, and map extent defined based on the input data.
168 Parameters
169 ----------
170 cube: Cube
171 2 dimensional (lat and lon) Cube of the data to plot.
172 figure:
173 Matplotlib Figure object holding all plot elements.
174 cmap:
175 Matplotlib colormap.
176 grid_size: (int, int), optional
177 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
178 subplot: int, optional
179 Subplot index if multiple spatial subplots in figure.
181 Returns
182 -------
183 axes:
184 Matplotlib GeoAxes definition.
185 """
186 # Identify min/max plot bounds.
187 try:
188 lat_axis, lon_axis = get_cube_yxcoordname(cube)
189 xmin = np.nanmin(cube.coord(lon_axis).points)
190 xmax = np.nanmax(cube.coord(lon_axis).points)
191 ymin = np.nanmin(cube.coord(lat_axis).points)
192 ymax = np.nanmax(cube.coord(lat_axis).points)
194 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
195 if np.abs(xmax - xmin) > 180.0:
196 xmin = xmin - 180.0
197 xmax = xmax - 180.0
198 logger.debug("Adjusting plot bounds to fit global extent.")
200 # Consider map projection orientation.
201 # Adapting orientation enables plotting across international dateline.
202 # Users can adapt the default central_longitude if alternative projections views.
203 if xmax > 180.0 or xmin < -180.0:
204 central_longitude = 180.0
205 else:
206 central_longitude = 0.0
208 # Define spatial map projection.
209 coord_system = cube.coord(lat_axis).coord_system
210 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
211 # Define rotated pole map projection for rotated pole inputs.
212 projection = ccrs.RotatedPole(
213 pole_longitude=coord_system.grid_north_pole_longitude,
214 pole_latitude=coord_system.grid_north_pole_latitude,
215 central_rotated_longitude=central_longitude,
216 )
217 crs = projection
218 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 218 ↛ 220line 218 didn't jump to line 220 because the condition on line 218 was never true
219 # Define Transverse Mercator projection for TM inputs.
220 projection = ccrs.TransverseMercator(
221 central_longitude=coord_system.longitude_of_central_meridian,
222 central_latitude=coord_system.latitude_of_projection_origin,
223 false_easting=coord_system.false_easting,
224 false_northing=coord_system.false_northing,
225 scale_factor=coord_system.scale_factor_at_central_meridian,
226 )
227 crs = projection
228 else:
229 # Assume polar projection for regional grids encompassing N. Pole
230 if ymin > 20.0 and ymax > 80.0:
231 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
232 elif ymin < -80.0 and ymax < -20.0:
233 projection = ccrs.SouthPolarStereo(central_longitude=central_longitude)
234 # Define regular map projection for non-rotated pole inputs.
235 # Alternatives might include e.g. for global model outputs:
236 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
237 # projection = ccrs.NearsidePerspective(
238 # central_longitude=180.0,
239 # central_latitude=0,
240 # satellite_height=35785831,
241 # )
242 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
243 else:
244 projection = ccrs.PlateCarree(central_longitude=central_longitude)
245 crs = ccrs.PlateCarree()
247 # Define axes for plot (or subplot) with required map projection.
248 if subplot is not None:
249 axes = figure.add_subplot(
250 grid_size[0], grid_size[1], subplot, projection=projection
251 )
252 else:
253 axes = figure.add_subplot(projection=projection)
255 # Add coastlines and borderlines if cube contains x and y map coordinates.
256 # Avoid adding lines for specific fixed ancillary spatial plots
257 if any(name in cube.name() for name in ("land_", "orography", "altitude")):
258 pass
259 else:
260 if cmap.name in ["viridis", "Greys"]:
261 coastcol = "magenta"
262 else:
263 coastcol = "black"
264 logger.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
265 axes.coastlines(resolution="10m", color=coastcol, alpha=0.8)
266 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol, alpha=0.3)
268 # Add gridlines.
269 gl = axes.gridlines(
270 alpha=0.3,
271 draw_labels=True,
272 dms=False,
273 x_inline=False,
274 y_inline=False,
275 )
276 gl.top_labels = False
277 gl.right_labels = False
278 if subplot:
279 gl.bottom_labels = False
280 gl.left_labels = False
281 if subplot % grid_size[1] == 1:
282 gl.left_labels = True
283 if subplot > ((grid_size[0] - 1) * grid_size[1]): 283 ↛ 288line 283 didn't jump to line 288 because the condition on line 283 was always true
284 gl.bottom_labels = True
286 # If is lat/lon spatial map, fix extent to keep plot tight.
287 # Specifying crs within set_extent helps ensure only data region is shown.
288 if isinstance(
289 coord_system, (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
290 ):
291 axes.set_extent([xmin, xmax, ymin, ymax], crs=crs)
293 except ValueError:
294 # Skip if not both x and y map coordinates.
295 axes = figure.gca()
297 return axes
300def _get_plot_resolution() -> int:
301 """Get resolution of rasterised plots in pixels per inch."""
302 return get_recipe_metadata().get("plot_resolution", 100)
305def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
306 """Return title and filename based on start and end points or bounds."""
307 if use_bounds and seq_coord.has_bounds():
308 vals = seq_coord.bounds.flatten()
309 else:
310 vals = seq_coord.points
311 start = seq_coord.units.title(vals[0])
312 end = seq_coord.units.title(vals[-1])
314 if start == end:
315 sequence_title = f"\n [{start}]"
316 sequence_fname = f"_{filename_slugify(start)}"
317 else:
318 sequence_title = f"\n [{start} to {end}]"
319 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
321 # Do not include time if coord set to zero.
322 if (
323 seq_coord.units == "hours since 0001-01-01 00:00:00"
324 and vals[0] == 0
325 and vals[-1] == 0
326 ):
327 sequence_title = ""
328 sequence_fname = ""
330 return sequence_title, sequence_fname
333def _set_title_and_filename(
334 seq_coord: iris.coords.Coord,
335 nplot: int,
336 recipe_title: str,
337 filename: str,
338 model_name: str | None = None,
339):
340 """Set plot title and filename based on cube coordinate.
342 Parameters
343 ----------
344 sequence_coordinate: iris.coords.Coord
345 Coordinate about which to make a plot sequence.
346 nplot: int
347 Number of output plots to generate - controls title/naming.
348 recipe_title: str
349 Default plot title, potentially to update.
350 filename: str
351 Input plot filename, potentially to update.
353 Returns
354 -------
355 plot_title: str
356 Output formatted plot title string, based on plotted data.
357 plot_filename: str
358 Output formatted plot filename string.
359 """
360 ndim = seq_coord.ndim
361 npoints = np.size(seq_coord.points)
362 sequence_title = ""
363 sequence_fname = ""
365 # Case 1: Multiple dimension sequence input - list number of aggregated cases
366 # (e.g. aggregation histogram plots)
367 if ndim > 1:
368 ncase = np.shape(seq_coord)[0]
369 sequence_title = f"\n [{ncase} cases]"
370 sequence_fname = f"_{ncase}cases"
372 # Case 2: Single dimension input
373 else:
374 # Single sequence point
375 if npoints == 1:
376 if nplot > 1:
377 # Default labels for sequence inputs
378 sequence_value = seq_coord.units.title(seq_coord.points[0])
379 sequence_value = sequence_value.replace(" unknown", "")
380 sequence_title = f"\n [{sequence_value}]"
381 sequence_fname = f"_{filename_slugify(sequence_value)}"
382 else:
383 # Aggregated attribute available where input collapsed over aggregation
384 try:
385 ncase = seq_coord.attributes["number_reference_times"]
386 sequence_title = f"\n [{ncase} cases]"
387 sequence_fname = f"_{ncase}cases"
388 except KeyError:
389 sequence_title, sequence_fname = _get_start_end_strings(
390 seq_coord, use_bounds=seq_coord.has_bounds()
391 )
392 # Multiple sequence (e.g. time) points
393 else:
394 sequence_title, sequence_fname = _get_start_end_strings(
395 seq_coord, use_bounds=False
396 )
398 # Set plot title and filename
399 plot_title = f"{recipe_title}{sequence_title}"
401 # Set plot filename, defaulting to user input if provided.
402 if filename is None:
403 filename = slugify(recipe_title)
404 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
405 else:
406 if nplot > 1:
407 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
408 else:
409 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
411 if model_name: 411 ↛ 412line 411 didn't jump to line 412 because the condition on line 411 was never true
412 if "the comparison model" in plot_title:
413 plot_title = plot_title.replace("the comparison model", model_name)
414 plot_filename = slugify(
415 plot_filename.replace("the comparison model", model_name)
416 )
417 else:
418 plot_title = f"{model_name}_{plot_title}"
419 plot_filename = slugify(f"{model_name}_{plot_filename}")
421 return plot_title, plot_filename
424def _select_series_coord(cube, series_coordinate):
425 """Determine the grid coordinates to use to calculate grid spacing."""
426 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
427 if series_coordinate in spacing_coordinates: 427 ↛ 433line 427 didn't jump to line 433 because the condition on line 427 was always true
428 # Try the requested coordinate first then the fallbacks in order.
429 fallbacks = [series_coordinate] + [
430 c for c in spacing_coordinates if c != series_coordinate
431 ]
432 else:
433 fallbacks = {series_coordinate}
435 # Try each possible coordinate.
436 for coord in fallbacks:
437 try:
438 return cube.coord(coord)
439 except iris.exceptions.CoordinateNotFoundError:
440 logger.debug("Coordinate %s not found.", coord)
442 # If we get here, none of the fallback options were found.
443 raise iris.exceptions.CoordinateNotFoundError(
444 f"No valid coordinate found for '{series_coordinate}' "
445 f"or fallback options {fallbacks}"
446 )
449def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
450 """Control postage stamp plot output titles based on stamp coordinate."""
451 if stamp_coord.name() == "realization":
452 mtitle = "Member"
453 else:
454 mtitle = stamp_coord.name().capitalize()
456 if stamp_coord.name() == "time":
457 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
458 else:
459 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
461 return mtitle
464def _set_axis_range(cubes):
465 """Get minimum and maximum from levels information."""
466 levels = None
467 for cube in cubes: 467 ↛ 483line 467 didn't jump to line 483 because the loop on line 467 didn't complete
468 # First check if user-specified "auto" range variable.
469 # This maintains the value of levels as None, so proceed.
470 _, levels, _ = colorbar_map_levels(cube, axis="y")
471 if levels is None:
472 break
473 # If levels is changed, recheck to use the vmin,vmax or
474 # levels-based ranges for histogram plots.
475 _, levels, _ = colorbar_map_levels(cube)
476 logger.debug("levels: %s", levels)
477 if levels is not None: 477 ↛ 467line 477 didn't jump to line 467 because the condition on line 477 was always true
478 vmin = min(levels)
479 vmax = max(levels)
480 logger.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
481 break
483 if levels is None:
484 vmin = min(cb.data.min() for cb in cubes)
485 vmax = max(cb.data.max() for cb in cubes)
487 return vmin, vmax
490def _find_matched_slices(cubes, sequence_coordinate):
491 """Identify matched cubes in CubeList by sequence_coordinate values.
493 Ensures common points are compared for multiple cube inputs.
494 """
495 all_points = sorted(
496 set(
497 itertools.chain.from_iterable(
498 cb.coord(sequence_coordinate).points for cb in cubes
499 )
500 )
501 )
502 all_slices = list(
503 itertools.chain.from_iterable(
504 cb.slices_over(sequence_coordinate) for cb in cubes
505 )
506 )
507 # Matched slices (matched by seq coord point; it may happen that
508 # evaluated models do not cover the same seq coord range, hence matching
509 # necessary)
510 cube_iterables = [
511 iris.cube.CubeList(
512 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
513 )
514 for point in all_points
515 ]
517 return cube_iterables
520def _plot_and_save_spatial_plot(
521 cube: iris.cube.Cube,
522 filename: str,
523 title: str,
524 method: Literal["contourf", "pcolormesh", "scatter"],
525 overlay_cube: iris.cube.Cube | None = None,
526 contour_cube: iris.cube.Cube | None = None,
527 point_cube: iris.cube.Cube | None = None,
528 **kwargs,
529):
530 """Plot and save a spatial plot.
532 Parameters
533 ----------
534 cube: Cube
535 2 dimensional (lat and lon) Cube of the data to plot.
536 filename: str
537 Filename of the plot to write.
538 title: str
539 Plot title.
540 method: "contourf" | "pcolormesh" | "scatter"
541 The plotting method to use
542 Select choice of "contourf" or "pcolormesh" for gridded data. Use "scatter" for point-based data.
543 overlay_cube: Cube, optional
544 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
545 contour_cube: Cube, optional
546 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
547 point_cube: Cube, optional
548 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
549 """
550 # Setup plot details, size, resolution, etc.
551 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
553 # Specify the color bar
554 cmap, levels, norm = colorbar_map_levels(cube)
556 # If overplotting, set required colorbars
557 if overlay_cube:
558 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
559 if contour_cube:
560 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
562 # Setup plot map projection, extent and coastlines and borderlines.
563 axes = _setup_spatial_map(cube, fig, cmap)
565 # Set colorscale bounds
566 try:
567 vmin = min(levels)
568 vmax = max(levels)
569 except TypeError:
570 vmin, vmax = None, None
571 # Ensure to use norm and not vmin/vmax if levels are defined.
572 if norm is not None:
573 vmin = None
574 vmax = None
575 logger.debug("Plotting using defined levels.")
577 # Plot the field.
578 if method == "contourf":
579 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
580 elif method == "pcolormesh":
581 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
582 elif method == "scatter":
583 # Scatter plot of the field. The marker size is chosen to give
584 # symbols that decrease in size as the number of data points
585 # increases, although the fraction of the figure covered by
586 # symbols increases roughly as N^(1/2), disregarding overlaps,
587 # and has been selected for the default figure size of (10, 10).
588 # Should this be changed, the marker size should be adjusted in
589 # proportion to the area of the figure.
590 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
591 lat_axis, lon_axis = get_cube_yxcoordname(cube)
592 plot = iplt.scatter(
593 cube.coord(lon_axis),
594 cube.coord(lat_axis),
595 c=cube.data[:],
596 s=mrk_size,
597 cmap=cmap,
598 edgecolors="k",
599 norm=norm,
600 vmin=vmin,
601 vmax=vmax,
602 )
603 else:
604 raise ValueError(f"Unknown plotting method: {method}")
606 # Overplot overlay field, if required
607 if overlay_cube:
608 try:
609 over_vmin = min(over_levels)
610 over_vmax = max(over_levels)
611 except TypeError:
612 over_vmin, over_vmax = None, None
613 if over_norm is not None: 613 ↛ 614line 613 didn't jump to line 614 because the condition on line 613 was never true
614 over_vmin = None
615 over_vmax = None
616 overlay = iplt.pcolormesh(
617 overlay_cube,
618 cmap=over_cmap,
619 norm=over_norm,
620 alpha=0.8,
621 vmin=over_vmin,
622 vmax=over_vmax,
623 )
624 # Overplot contour field, if required, with contour labelling.
625 if contour_cube:
626 contour = iplt.contour(
627 contour_cube,
628 colors="darkgray",
629 levels=cntr_levels,
630 norm=cntr_norm,
631 alpha=0.5,
632 linestyles="--",
633 linewidths=1,
634 )
635 plt.clabel(contour)
636 # Overplot valid elements of point-based field, if required.
637 # Check for non-masked points only to avoid plotting missing data.
638 if point_cube:
639 mrk_size = int(np.sqrt(2500000.0 / len(point_cube.data)))
640 lat_axis, lon_axis = get_cube_yxcoordname(point_cube)
641 lon_coord = point_cube.coord(lon_axis)
642 lat_coord = point_cube.coord(lat_axis)
643 valid = ~point_cube.data.mask
644 valid_lon = iris.coords.AuxCoord(
645 lon_coord.points[valid],
646 standard_name=lon_coord.standard_name,
647 units=lon_coord.units,
648 coord_system=lon_coord.coord_system,
649 )
650 valid_lat = iris.coords.AuxCoord(
651 lat_coord.points[valid],
652 standard_name=lat_coord.standard_name,
653 units=lat_coord.units,
654 coord_system=lat_coord.coord_system,
655 )
656 iplt.scatter(
657 valid_lon,
658 valid_lat,
659 c=point_cube.data[valid],
660 s=mrk_size,
661 cmap=cmap,
662 edgecolors="k",
663 norm=norm,
664 vmin=vmin,
665 vmax=vmax,
666 )
668 # Check to see if transect, and if so, adjust y axis.
669 if is_transect(cube):
670 if "pressure" in [coord.name() for coord in cube.coords()]:
671 axes.invert_yaxis()
672 axes.set_yscale("log")
673 axes.set_ylim(1100, 100)
674 # If both model_level_number and level_height exists, iplt can construct
675 # plot as a function of height above orography (NOT sea level).
676 elif {"model_level_number", "level_height"}.issubset( 676 ↛ 681line 676 didn't jump to line 681 because the condition on line 676 was always true
677 {coord.name() for coord in cube.coords()}
678 ):
679 axes.set_yscale("log")
681 axes.set_title(
682 f"{title}\n"
683 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
684 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
685 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
686 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
687 fontsize=16,
688 )
690 # Inset code
691 axins = inset_axes(
692 axes,
693 width="20%",
694 height="20%",
695 loc="upper right",
696 axes_class=GeoAxes,
697 axes_kwargs={"map_projection": ccrs.PlateCarree()},
698 )
700 # Slightly transparent to reduce plot blocking.
701 axins.patch.set_alpha(0.4)
703 axins.coastlines(resolution="50m")
704 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
706 SLat, SLon, ELat, ELon = (
707 float(coord) for coord in cube.attributes["transect_coords"].split("_")
708 )
710 # Draw line between them
711 axins.plot(
712 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
713 )
715 # Plot points (note: lon, lat order for Cartopy)
716 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
717 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
719 lon_min, lon_max = sorted([SLon, ELon])
720 lat_min, lat_max = sorted([SLat, ELat])
722 # Midpoints
723 lon_mid = (lon_min + lon_max) / 2
724 lat_mid = (lat_min + lat_max) / 2
726 # Maximum half-range
727 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
728 if half_range == 0: # points identical → provide small default 728 ↛ 732line 728 didn't jump to line 732 because the condition on line 728 was always true
729 half_range = 1
731 # Set square extent
732 axins.set_extent(
733 [
734 lon_mid - half_range,
735 lon_mid + half_range,
736 lat_mid - half_range,
737 lat_mid + half_range,
738 ],
739 crs=ccrs.PlateCarree(),
740 )
742 # Ensure square aspect
743 axins.set_aspect("equal")
745 else:
746 # Add title.
747 axes.set_title(title, fontsize=16)
749 # Adjust padding if spatial plot or transect
750 if is_transect(cube):
751 yinfopad = -0.1
752 ycbarpad = 0.1
753 else:
754 yinfopad = 0.01
755 ycbarpad = 0.042
757 # Add watermark with min/max/mean. Currently not user togglable.
758 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
759 cube_min, cube_max, cube_mean = calc_array_stats(cube.data)
760 axes.annotate(
761 f"Min: {cube_min:.3g} Max: {cube_max:.3g} Mean: {cube_mean:.3g}",
762 xy=(0.025, yinfopad),
763 xycoords="axes fraction",
764 xytext=(-5, 5),
765 textcoords="offset points",
766 ha="left",
767 va="bottom",
768 size=11,
769 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
770 )
772 # Add secondary colour bar for overlay_cube field if required.
773 if overlay_cube:
774 cbarB = fig.colorbar(
775 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
776 )
777 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
778 # add ticks and tick_labels for every levels if less than 20 levels exist
779 if over_levels is not None and len(over_levels) < 20: 779 ↛ 780line 779 didn't jump to line 780 because the condition on line 779 was never true
780 cbarB.set_ticks(over_levels)
781 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
782 if any(
783 var in overlay_cube.name()
784 for var in ("rainfall", "snowfall", "visibility")
785 ):
786 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
787 logger.debug("Set secondary colorbar ticks and labels.")
789 # Add main colour bar.
790 cbar = fig.colorbar(
791 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
792 )
794 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
795 # add ticks and tick_labels for every levels if less than 20 levels exist
796 if levels is not None and len(levels) < 20:
797 cbar.set_ticks(levels)
798 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
799 if any(var in cube.name() for var in ("rainfall", "snowfall", "visibility")): 799 ↛ 802line 799 didn't jump to line 802 because the condition on line 799 was always true
800 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
801 # Tick labels for rainfall rates from Nimrod radar data.
802 if "rainfall rate composite" in cube.name(): 802 ↛ 803line 802 didn't jump to line 803 because the condition on line 802 was never true
803 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
804 # Tick labels for rain accumulations from Nimrod radar data.
805 if "rain accumulation" in cube.name(): 805 ↛ 806line 805 didn't jump to line 806 because the condition on line 805 was never true
806 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
807 if "wts accumulation" in cube.name(): 807 ↛ 808line 807 didn't jump to line 808 because the condition on line 807 was never true
808 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
809 cbar.minorticks_off()
810 cbar.set_ticks(tick_levels)
811 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
812 cbar.set_label(label=f"{cube.name()}", size=14)
813 # Tick labels for model rainfall data.
814 if "surface_microphysical" in cube.name(): 814 ↛ 817line 814 didn't jump to line 817 because the condition on line 814 was always true
815 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
816 # Tick labels for Nimrod weights data.
817 logger.debug("Set colorbar ticks and labels.")
819 # Save plot.
820 _save_close_figure(fig, "spatial", filename)
823def _plot_and_save_postage_stamp_spatial_plot(
824 cube: iris.cube.Cube,
825 filename: str,
826 stamp_coordinate: str,
827 title: str,
828 method: Literal["contourf", "pcolormesh"],
829 overlay_cube: iris.cube.Cube | None = None,
830 contour_cube: iris.cube.Cube | None = None,
831 **kwargs,
832):
833 """Plot postage stamp spatial plots from an ensemble.
835 Parameters
836 ----------
837 cube: Cube
838 Iris cube of data to be plotted. It must have the stamp coordinate.
839 filename: str
840 Filename of the plot to write.
841 stamp_coordinate: str
842 Coordinate that becomes different plots.
843 method: "contourf" | "pcolormesh"
844 The plotting method to use.
845 overlay_cube: Cube, optional
846 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
847 contour_cube: Cube, optional
848 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
850 Raises
851 ------
852 ValueError
853 If the cube doesn't have the right dimensions.
854 """
855 # Use the smallest square grid that will fit the members.
856 nmember = len(cube.coord(stamp_coordinate).points)
857 grid_rows = int(math.sqrt(nmember))
858 grid_size = math.ceil(nmember / grid_rows)
860 fig = plt.figure(
861 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
862 )
864 # Specify the color bar
865 cmap, levels, norm = colorbar_map_levels(cube)
866 # If overplotting, set required colorbars
867 if overlay_cube: 867 ↛ 868line 867 didn't jump to line 868 because the condition on line 867 was never true
868 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
869 if contour_cube: 869 ↛ 870line 869 didn't jump to line 870 because the condition on line 869 was never true
870 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
872 # Make a subplot for each member.
873 for member, subplot in zip(
874 cube.slices_over(stamp_coordinate),
875 range(1, grid_size * grid_rows + 1),
876 strict=False,
877 ):
878 # Setup subplot map projection, extent and coastlines and borderlines.
879 axes = _setup_spatial_map(
880 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
881 )
882 if method == "contourf":
883 # Filled contour plot of the field.
884 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
885 elif method == "pcolormesh":
886 if levels is not None:
887 vmin = min(levels)
888 vmax = max(levels)
889 else:
890 raise TypeError("Unknown vmin and vmax range.")
891 vmin, vmax = None, None
892 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
893 # if levels are defined.
894 if norm is not None: 894 ↛ 895line 894 didn't jump to line 895 because the condition on line 894 was never true
895 vmin = None
896 vmax = None
897 # pcolormesh plot of the field.
898 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
899 else:
900 raise ValueError(f"Unknown plotting method: {method}")
902 # Overplot overlay field, if required
903 if overlay_cube: 903 ↛ 904line 903 didn't jump to line 904 because the condition on line 903 was never true
904 try:
905 over_vmin = min(over_levels)
906 over_vmax = max(over_levels)
907 except TypeError:
908 over_vmin, over_vmax = None, None
909 if over_norm is not None:
910 over_vmin = None
911 over_vmax = None
912 iplt.pcolormesh(
913 overlay_cube[member.coord(stamp_coordinate).points[0]],
914 cmap=over_cmap,
915 norm=over_norm,
916 alpha=0.6,
917 vmin=over_vmin,
918 vmax=over_vmax,
919 )
920 # Overplot contour field, if required
921 if contour_cube: 921 ↛ 922line 921 didn't jump to line 922 because the condition on line 921 was never true
922 iplt.contour(
923 contour_cube[member.coord(stamp_coordinate).points[0]],
924 colors="darkgray",
925 levels=cntr_levels,
926 norm=cntr_norm,
927 alpha=0.6,
928 linestyles="--",
929 linewidths=1,
930 )
931 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
932 axes.set_title(f"{mtitle}")
934 # Put the shared colorbar in its own axes.
935 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
936 colorbar = fig.colorbar(
937 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
938 )
939 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
941 # Overall figure title.
942 fig.suptitle(title, fontsize=16)
944 # Save plot.
945 _save_close_figure(fig, "contour postate stamp", filename)
948def _plot_and_save_line_series(
949 cubes: iris.cube.CubeList,
950 coords: list[iris.coords.Coord],
951 ensemble_coord: str,
952 filename: str,
953 title: str,
954 **kwargs,
955):
956 """Plot and save a 1D line series.
958 Parameters
959 ----------
960 cubes: Cube or CubeList
961 Cube or CubeList containing the cubes to plot on the y-axis.
962 coords: list[Coord]
963 Coordinates to plot on the x-axis, one per cube.
964 ensemble_coord: str
965 Ensemble coordinate in the cube.
966 filename: str
967 Filename of the plot to write.
968 title: str
969 Plot title.
970 """
971 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
973 model_colors_map = get_model_colors_map(cubes)
975 # Store min/max ranges.
976 y_levels = []
978 # Check match-up across sequence coords gives consistent sizes
979 validate_cubes_coords(cubes, coords)
981 for cube, coord in zip(cubes, coords, strict=True):
982 label = None
983 color = "black"
984 if model_colors_map:
985 label = cube.attributes.get("model_name")
986 color = model_colors_map.get(label)
987 if not cube.coords(ensemble_coord): 987 ↛ 989line 987 didn't jump to line 989 because the condition on line 987 was never true
988 # No ensemble coordinate — plot the cube directly as a single line.
989 iplt.plot(coord, cube, color=color, marker="o", ls="-", lw=3, label=label)
990 else:
991 for cube_slice in cube.slices_over(ensemble_coord):
992 # Label with (control) if part of an ensemble or not otherwise.
993 if cube_slice.coord(ensemble_coord).points == [0]:
994 iplt.plot(
995 coord,
996 cube_slice,
997 color=color,
998 marker="o",
999 ls="-",
1000 lw=3,
1001 label=f"{label} (control)"
1002 if len(cube.coord(ensemble_coord).points) > 1
1003 else label,
1004 )
1005 # Label with (perturbed) if part of an ensemble and not the control.
1006 else:
1007 iplt.plot(
1008 coord,
1009 cube_slice,
1010 color=color,
1011 ls="-",
1012 lw=1.5,
1013 alpha=0.75,
1014 label=f"{label} (member)",
1015 )
1017 # Calculate the global min/max if multiple cubes are given.
1018 _, levels, _ = colorbar_map_levels(cube, axis="y")
1019 if levels is not None: 1019 ↛ 1020line 1019 didn't jump to line 1020 because the condition on line 1019 was never true
1020 y_levels.append(min(levels))
1021 y_levels.append(max(levels))
1023 # Get the current axes.
1024 ax = plt.gca()
1026 # Add some labels and tweak the style.
1027 # check if cubes[0] works for single cube if not CubeList
1028 if coords[0].name() == "time":
1029 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1030 else:
1031 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1032 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1033 ax.set_title(title, fontsize=16)
1035 ax.ticklabel_format(axis="y", useOffset=False)
1036 ax.tick_params(axis="x", labelrotation=15)
1037 ax.tick_params(axis="both", labelsize=12)
1039 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1040 if y_levels: 1040 ↛ 1041line 1040 didn't jump to line 1041 because the condition on line 1040 was never true
1041 ax.set_ylim(min(y_levels), max(y_levels))
1042 logger.debug("Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels))
1043 else:
1044 ax.autoscale()
1046 # Add gridlines
1047 ax.grid(linestyle="--", color="grey", linewidth=1)
1048 # Add zero line
1049 ymin, ymax = ax.get_ylim()
1050 if ymin < 0.0 and ymax > 0.0:
1051 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1052 # Identify unique labels for legend
1053 handles = list(
1054 {
1055 label: handle
1056 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1057 }.values()
1058 )
1059 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1061 # Save plot.
1062 _save_close_figure(fig, "line", filename)
1065def _plot_and_save_line_power_spectrum_series(
1066 cubes: iris.cube.Cube | iris.cube.CubeList,
1067 coords: list[iris.coords.Coord],
1068 ensemble_coord: str,
1069 filename: str,
1070 title: str,
1071 series_coordinate: str,
1072 **kwargs,
1073):
1074 """Plot and save a 1D line series.
1076 Parameters
1077 ----------
1078 cubes: Cube or CubeList
1079 Cube or CubeList containing the cubes to plot on the y-axis.
1080 coords: list[Coord]
1081 Coordinates to plot on the x-axis, one per cube.
1082 ensemble_coord: str
1083 Ensemble coordinate in the cube.
1084 filename: str
1085 Filename of the plot to write.
1086 title: str
1087 Plot title.
1088 series_coordinate: str
1089 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1090 """
1091 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1092 model_colors_map = get_model_colors_map(cubes)
1093 ax = plt.gca()
1095 # Store min/max ranges.
1096 y_levels = []
1098 line_marker = None
1099 line_width = 1
1101 for cube in iter_maybe(cubes):
1102 # next 2 lines replace chunk of code.
1103 xcoord = _select_series_coord(cube, series_coordinate)
1104 xname = xcoord.points
1106 yfield = cube.data # power spectrum
1108 # If data from power spectra is all np.nans (like T+0h rainfall field which
1109 # might be full of zeros), then set yfield to zeros so it doesn't crash the
1110 # plotting.
1111 if np.all(np.isnan(yfield)):
1112 yfield = np.zeros_like(yfield)
1114 label = None
1115 color = "black"
1116 if model_colors_map: 1116 ↛ 1119line 1116 didn't jump to line 1119 because the condition on line 1116 was always true
1117 label = cube.attributes.get("model_name")
1118 color = model_colors_map.get(label)
1119 for cube_slice in cube.slices_over(ensemble_coord):
1120 # Label with (control) if part of an ensemble or not otherwise.
1121 if cube_slice.coord(ensemble_coord).points == [0]: 1121 ↛ 1135line 1121 didn't jump to line 1135 because the condition on line 1121 was always true
1122 ax.plot(
1123 xname,
1124 yfield,
1125 color=color,
1126 marker=line_marker,
1127 ls="-",
1128 lw=line_width,
1129 label=f"{label} (control)"
1130 if len(cube.coord(ensemble_coord).points) > 1
1131 else label,
1132 )
1133 # Label with (perturbed) if part of an ensemble and not the control.
1134 else:
1135 ax.plot(
1136 xname,
1137 yfield,
1138 color=color,
1139 ls="-",
1140 lw=1.5,
1141 alpha=0.75,
1142 label=f"{label} (member)",
1143 )
1145 # Calculate the global min/max if multiple cubes are given.
1146 _, levels, _ = colorbar_map_levels(cube, axis="y")
1147 if levels is not None: 1147 ↛ 1148line 1147 didn't jump to line 1148 because the condition on line 1147 was never true
1148 y_levels.append(min(levels))
1149 y_levels.append(max(levels))
1151 # Add some labels and tweak the style.
1153 title = f"{title}"
1154 ax.set_title(title, fontsize=16)
1156 # Set appropriate x-axis label based on coordinate
1157 if series_coordinate == "wavelength" or ( 1157 ↛ 1160line 1157 didn't jump to line 1160 because the condition on line 1157 was never true
1158 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1159 ):
1160 ax.set_xlabel("Wavelength (km)", fontsize=14)
1161 elif series_coordinate == "physical_wavenumber" or ( 1161 ↛ 1164line 1161 didn't jump to line 1164 because the condition on line 1161 was never true
1162 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1163 ):
1164 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1165 else: # frequency or check units
1166 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1166 ↛ 1167line 1166 didn't jump to line 1167 because the condition on line 1166 was never true
1167 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1168 else:
1169 ax.set_xlabel("Wavenumber", fontsize=14)
1171 ax.set_ylabel("Power Spectral Density", fontsize=14)
1172 ax.tick_params(axis="both", labelsize=12)
1174 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1176 # Set log-log scale
1177 ax.set_xscale("log")
1178 ax.set_yscale("log")
1180 # Add gridlines
1181 ax.grid(linestyle="--", color="grey", linewidth=1)
1182 # Ientify unique labels for legend
1183 handles = list(
1184 {
1185 label: handle
1186 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1187 }.values()
1188 )
1189 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1191 # Save plot.
1192 _save_close_figure(fig, "line power spectrum", filename)
1195def _plot_and_save_vertical_line_series(
1196 cubes: iris.cube.CubeList,
1197 coords: list[iris.coords.Coord],
1198 ensemble_coord: str,
1199 filename: str,
1200 series_coordinate: str,
1201 title: str,
1202 vmin: float,
1203 vmax: float,
1204 **kwargs,
1205):
1206 """Plot and save a 1D line series in vertical.
1208 Parameters
1209 ----------
1210 cubes: CubeList
1211 1 dimensional Cube or CubeList of the data to plot on x-axis.
1212 coord: list[Coord]
1213 Coordinates to plot on the y-axis, one per cube.
1214 ensemble_coord: str
1215 Ensemble coordinate in the cube.
1216 filename: str
1217 Filename of the plot to write.
1218 series_coordinate: str
1219 Coordinate to use as vertical axis.
1220 title: str
1221 Plot title.
1222 vmin: float
1223 Minimum value for the x-axis.
1224 vmax: float
1225 Maximum value for the x-axis.
1226 """
1227 # plot the vertical pressure axis using log scale
1228 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1230 model_colors_map = get_model_colors_map(cubes)
1232 # Check match-up across sequence coords gives consistent sizes
1233 validate_cubes_coords(cubes, coords)
1235 for cube, coord in zip(cubes, coords, strict=True):
1236 label = None
1237 color = "black"
1238 if model_colors_map: 1238 ↛ 1239line 1238 didn't jump to line 1239 because the condition on line 1238 was never true
1239 label = cube.attributes.get("model_name")
1240 color = model_colors_map.get(label)
1242 for cube_slice in cube.slices_over(ensemble_coord):
1243 # If ensemble data given plot control member with (control)
1244 # unless single forecast.
1245 if cube_slice.coord(ensemble_coord).points == [0]:
1246 iplt.plot(
1247 cube_slice,
1248 coord,
1249 color=color,
1250 marker="o",
1251 ls="-",
1252 lw=3,
1253 label=f"{label} (control)"
1254 if len(cube.coord(ensemble_coord).points) > 1
1255 else label,
1256 )
1257 # If ensemble data given plot perturbed members with (perturbed).
1258 else:
1259 iplt.plot(
1260 cube_slice,
1261 coord,
1262 color=color,
1263 ls="-",
1264 lw=1.5,
1265 alpha=0.75,
1266 label=f"{label} (member)",
1267 )
1269 # Get the current axis
1270 ax = plt.gca()
1272 # Special handling for pressure level data.
1273 if series_coordinate == "pressure": 1273 ↛ 1295line 1273 didn't jump to line 1295 because the condition on line 1273 was always true
1274 # Invert y-axis and set to log scale.
1275 ax.invert_yaxis()
1276 ax.set_yscale("log")
1278 # Define y-ticks and labels for pressure log axis.
1279 y_tick_labels = [
1280 "1000",
1281 "850",
1282 "700",
1283 "500",
1284 "300",
1285 "200",
1286 "100",
1287 ]
1288 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1290 # Set y-axis limits and ticks.
1291 ax.set_ylim(1100, 100)
1293 # Test if series_coordinate is model level data. The UM data uses
1294 # model_level_number and lfric uses full_levels as coordinate.
1295 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1296 # Define y-ticks and labels for vertical axis.
1297 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1298 y_tick_labels = [str(int(i)) for i in y_ticks]
1299 ax.set_ylim(min(y_ticks), max(y_ticks))
1301 ax.set_yticks(y_ticks)
1302 ax.set_yticklabels(y_tick_labels)
1304 # Set x-axis limits.
1305 ax.set_xlim(vmin, vmax)
1306 # Mark y=0 if present in plot.
1307 if vmin < 0.0 and vmax > 0.0: 1307 ↛ 1308line 1307 didn't jump to line 1308 because the condition on line 1307 was never true
1308 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1310 # Add some labels and tweak the style.
1311 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1312 ax.set_xlabel(
1313 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1314 )
1315 ax.set_title(title, fontsize=16)
1316 ax.ticklabel_format(axis="x")
1317 ax.tick_params(axis="y")
1318 ax.tick_params(axis="both", labelsize=12)
1320 # Add gridlines
1321 ax.grid(linestyle="--", color="grey", linewidth=1)
1322 # Ientify unique labels for legend
1323 handles = list(
1324 {
1325 label: handle
1326 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1327 }.values()
1328 )
1329 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1331 # Save plot.
1332 _save_close_figure(fig, "vertical line", filename)
1335def _plot_and_save_scatter_plot(
1336 cube_x: iris.cube.Cube | iris.cube.CubeList,
1337 cube_y: iris.cube.Cube | iris.cube.CubeList,
1338 filename: str,
1339 title: str,
1340 one_to_one: bool,
1341 model_names: list[str] | None = None,
1342 **kwargs,
1343):
1344 """Plot and save a 2D scatter plot.
1346 Parameters
1347 ----------
1348 cube_x: Cube | CubeList
1349 1 dimensional Cube or CubeList of the data to plot on x-axis.
1350 cube_y: Cube | CubeList
1351 1 dimensional Cube or CubeList of the data to plot on y-axis.
1352 filename: str
1353 Filename of the plot to write.
1354 title: str
1355 Plot title.
1356 one_to_one: bool
1357 Whether a 1:1 line is plotted.
1358 """
1359 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1360 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1361 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1362 # over the pairs simultaneously.
1364 # Ensure cube_x and cube_y are iterable
1365 cube_x_iterable = iter_maybe(cube_x)
1366 cube_y_iterable = iter_maybe(cube_y)
1368 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1369 iplt.scatter(cube_x_iter, cube_y_iter)
1370 if one_to_one is True:
1371 plt.plot(
1372 [
1373 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1374 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1375 ],
1376 [
1377 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1378 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1379 ],
1380 "k",
1381 linestyle="--",
1382 )
1383 ax = plt.gca()
1385 # Add some labels and tweak the style.
1386 if model_names is None:
1387 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1388 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1389 else:
1390 # Add the model names, these should be order of base (x) and other (y).
1391 ax.set_xlabel(
1392 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1393 )
1394 ax.set_ylabel(
1395 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1396 )
1397 ax.set_title(title, fontsize=16)
1398 ax.ticklabel_format(axis="y", useOffset=False)
1399 ax.tick_params(axis="x", labelrotation=15)
1400 ax.tick_params(axis="both", labelsize=12)
1401 ax.autoscale()
1403 # Save plot.
1404 _save_close_figure(fig, "scatter", filename)
1407def _plot_and_save_vector_plot(
1408 cube_u: iris.cube.Cube,
1409 cube_v: iris.cube.Cube,
1410 filename: str,
1411 title: str,
1412 method: Literal["contourf", "pcolormesh"],
1413 **kwargs,
1414):
1415 """Plot and save a 2D vector plot.
1417 Parameters
1418 ----------
1419 cube_u: Cube
1420 2 dimensional Cube of u component of the data.
1421 cube_v: Cube
1422 2 dimensional Cube of v component of the data.
1423 filename: str
1424 Filename of the plot to write.
1425 title: str
1426 Plot title.
1427 """
1428 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1429 # Create a cube containing the magnitude of the vector field.
1430 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1431 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1432 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1433 cube_vec_mag.rename(
1434 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1435 )
1437 # Specify the color bar
1438 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1440 # Setup plot map projection, extent and coastlines and borderlines.
1441 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1443 if method == "contourf":
1444 # Filled contour plot of the field.
1445 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1446 elif method == "pcolormesh":
1447 try:
1448 vmin = min(levels)
1449 vmax = max(levels)
1450 except TypeError:
1451 vmin, vmax = None, None
1452 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1453 # if levels are defined.
1454 if norm is not None:
1455 vmin = None
1456 vmax = None
1457 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1458 else:
1459 raise ValueError(f"Unknown plotting method: {method}")
1461 # Check to see if transect, and if so, adjust y axis.
1462 if is_transect(cube_vec_mag):
1463 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1464 axes.invert_yaxis()
1465 axes.set_yscale("log")
1466 axes.set_ylim(1100, 100)
1467 # If both model_level_number and level_height exists, iplt can construct
1468 # plot as a function of height above orography (NOT sea level).
1469 elif {"model_level_number", "level_height"}.issubset(
1470 {coord.name() for coord in cube_vec_mag.coords()}
1471 ):
1472 axes.set_yscale("log")
1474 axes.set_title(
1475 f"{title}\n"
1476 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1477 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1478 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1479 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1480 fontsize=16,
1481 )
1483 else:
1484 # Add title.
1485 axes.set_title(title, fontsize=16)
1487 # Add watermark with min/max/mean. Currently not user togglable.
1488 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1489 cube_min, cube_max, cube_mean = calc_array_stats(cube_vec_mag.data)
1490 axes.annotate(
1491 f"Min: {cube_min:.3g} Max: {cube_max:.3g} Mean: {cube_mean:.3g}",
1492 xy=(0.05, -0.05),
1493 xycoords="axes fraction",
1494 xytext=(-5, 5),
1495 textcoords="offset points",
1496 ha="right",
1497 va="bottom",
1498 size=11,
1499 bbox={"boxstyle": "round", "fc": "#cccccc", "ec": "#808080", "alpha": 0.9},
1500 )
1502 # Add colour bar.
1503 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1504 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1505 # add ticks and tick_labels for every levels if less than 20 levels exist
1506 if levels is not None and len(levels) < 20:
1507 cbar.set_ticks(levels)
1508 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1510 # 30 barbs along the longest axis of the plot, or a barb per point for data
1511 # with less than 30 points.
1512 step = max(max(cube_u.shape) // 30, 1)
1513 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1515 # Save plot.
1516 _save_close_figure(fig, "vector", filename)
1519def _plot_and_save_histogram_series(
1520 cubes: iris.cube.Cube | iris.cube.CubeList,
1521 filename: str,
1522 title: str,
1523 vmin: float,
1524 vmax: float,
1525 **kwargs,
1526):
1527 """Plot and save a histogram series.
1529 Parameters
1530 ----------
1531 cubes: Cube or CubeList
1532 2 dimensional Cube or CubeList of the data to plot as histogram.
1533 filename: str
1534 Filename of the plot to write.
1535 title: str
1536 Plot title.
1537 vmin: float
1538 minimum for colorbar
1539 vmax: float
1540 maximum for colorbar
1541 """
1542 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1543 ax = plt.gca()
1545 model_colors_map = get_model_colors_map(cubes)
1547 # Set default that histograms will produce probability density function
1548 # at each bin (integral over range sums to 1).
1549 density = True
1551 for cube in iter_maybe(cubes):
1552 # Easier to check title (where var name originates)
1553 # than seeing if long names exist etc.
1554 # Exception case, where distribution better fits log scales/bins.
1555 if (
1556 ("surface_microphysical" in title)
1557 or ("rain accumulation" in title)
1558 or ("Rainfall rate Composite" in title)
1559 or ("Nimrod_5min" in title)
1560 ):
1561 if "amount" in title:
1562 # Compute histogram following Klingaman et al. (2017): ASoP
1563 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1564 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1565 density = False
1566 else:
1567 bins = 10.0 ** (
1568 np.arange(-10, 27, 1) / 10.0
1569 ) # Suggestion from RMED toolbox.
1570 bins = np.insert(bins, 0, 0)
1571 ax.set_yscale("log")
1572 vmin = bins[1]
1573 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1574 ax.set_xscale("log")
1575 elif "lightning" in title:
1576 bins = [0, 1, 2, 3, 4, 5]
1577 else:
1578 bins = np.linspace(vmin, vmax, 51)
1579 logger.debug(
1580 "Plotting histogram with %s bins %s - %s.",
1581 np.size(bins),
1582 np.min(bins),
1583 np.max(bins),
1584 )
1586 # Reshape cube data into a single array to allow for a single histogram.
1587 # Otherwise we plot xdim histograms stacked.
1588 cube_data_1d = (cube.data).flatten()
1590 label = None
1591 color = "black"
1592 if model_colors_map:
1593 label = cube.attributes.get("model_name")
1594 color = model_colors_map[label]
1595 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1597 # Compute area under curve.
1598 if (
1599 ("surface_microphysical" in title and "amount" in title)
1600 or ("rain_accumulation" in title)
1601 or ("Rainfall rate Composite" in title)
1602 or ("Nimrod_5min" in title)
1603 ):
1604 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1605 x = x * bin_mean / x.sum()
1606 x = x[1:]
1607 y = y[1:]
1609 ax.plot(
1610 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1611 )
1613 # Add some labels and tweak the style.
1614 ax.set_title(title, fontsize=16)
1615 ax.set_xlabel(
1616 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1617 )
1618 ax.set_ylabel("Normalised probability density", fontsize=14)
1619 if (
1620 ("surface_microphysical" in title and "amount" in title)
1621 or ("rain accumulation" in title)
1622 or ("Nimrod_5min" in title)
1623 ):
1624 ax.set_ylabel(
1625 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1626 )
1627 try:
1628 ax.set_xlim(vmin, vmax)
1629 except ValueError:
1630 pass
1631 ax.tick_params(axis="both", labelsize=12)
1633 # Overlay grid-lines onto histogram plot.
1634 ax.grid(linestyle="--", color="grey", linewidth=1)
1635 if model_colors_map:
1636 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1638 # Save plot.
1639 _save_close_figure(fig, "histogram", filename)
1642def _plot_and_save_postage_stamp_histogram_series(
1643 cube: iris.cube.Cube,
1644 filename: str,
1645 title: str,
1646 stamp_coordinate: str,
1647 vmin: float,
1648 vmax: float,
1649 **kwargs,
1650):
1651 """Plot and save postage (ensemble members) stamps for a histogram series.
1653 Parameters
1654 ----------
1655 cube: Cube
1656 2 dimensional Cube of the data to plot as histogram.
1657 filename: str
1658 Filename of the plot to write.
1659 title: str
1660 Plot title.
1661 stamp_coordinate: str
1662 Coordinate that becomes different plots.
1663 vmin: float
1664 minimum for pdf x-axis
1665 vmax: float
1666 maximum for pdf x-axis
1667 """
1668 # Use the smallest square grid that will fit the members.
1669 nmember = len(cube.coord(stamp_coordinate).points)
1670 grid_rows = int(math.sqrt(nmember))
1671 grid_size = math.ceil(nmember / grid_rows)
1673 fig = plt.figure(
1674 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1675 )
1676 # Make a subplot for each member.
1677 for member, subplot in zip(
1678 cube.slices_over(stamp_coordinate),
1679 range(1, grid_size * grid_rows + 1),
1680 strict=False,
1681 ):
1682 # Implicit interface is much easier here, due to needing to have the
1683 # cartopy GeoAxes generated.
1684 plt.subplot(grid_rows, grid_size, subplot)
1685 # Reshape cube data into a single array to allow for a single histogram.
1686 # Otherwise we plot xdim histograms stacked.
1687 member_data_1d = (member.data).flatten()
1688 plt.hist(member_data_1d, density=True, stacked=True)
1689 axes = plt.gca()
1690 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1691 axes.set_title(f"{mtitle}")
1692 axes.set_xlim(vmin, vmax)
1694 # Overall figure title.
1695 fig.suptitle(title, fontsize=16)
1697 # Save plot.
1698 _save_close_figure(fig, "histogram postage stamp", filename)
1701def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1702 cube: iris.cube.Cube,
1703 filename: str,
1704 title: str,
1705 stamp_coordinate: str,
1706 vmin: float,
1707 vmax: float,
1708 **kwargs,
1709):
1710 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1711 ax.set_title(title, fontsize=16)
1712 ax.set_xlim(vmin, vmax)
1713 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1714 ax.set_ylabel("normalised probability density", fontsize=14)
1715 # Loop over all slices along the stamp_coordinate
1716 for member in cube.slices_over(stamp_coordinate):
1717 # Flatten the member data to 1D
1718 member_data_1d = member.data.flatten()
1719 # Plot the histogram using plt.hist
1720 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1721 plt.hist(
1722 member_data_1d,
1723 density=True,
1724 stacked=True,
1725 label=f"{mtitle}",
1726 )
1728 # Add a legend
1729 ax.legend(fontsize=16)
1731 # Save plot.
1732 _save_close_figure(fig, "histogram postage stamp", filename)
1735def _plot_and_save_scatter_series(
1736 cubes: iris.cube.Cube | iris.cube.CubeList,
1737 filename: str,
1738 title: str,
1739 vmin: float,
1740 vmax: float,
1741 hexbin: bool,
1742 **kwargs,
1743):
1744 """Plot and save a scatter plot series.
1746 Parameters
1747 ----------
1748 cubes: Cube or CubeList
1749 2 dimensional Cube or CubeList of the data to plot as scatter.
1750 filename: str
1751 Filename of the plot to write.
1752 title: str
1753 Plot title.
1754 vmin: float
1755 minimum for colorbar
1756 vmax: float
1757 maximum for colorbar
1758 hexbin: bool
1759 Flag to set output scatter generated as a hexbin frequency distribution plot of 2 cubes on single plot.
1760 Else scatter of all points, with potential to overplot many comparisons on same plot.
1761 """
1762 if hexbin:
1763 # Check cubes using same functionality as the difference operator.
1764 if len(cubes) != 2:
1765 raise ValueError(
1766 "Cubes should contain exactly 2 cubes for hexbin plotting."
1767 )
1768 title = title.replace("scatter", "hexbin")
1769 filename = filename.replace("scatter", "hexbin")
1771 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1772 ax = plt.gca()
1774 model_colors_map = get_model_colors_map(cubes)
1776 percentiles = np.arange(0, 100, 5)
1777 percentiles[0] = 1
1778 percentiles[-1] = 99
1779 quantiles = iris.cube.CubeList()
1781 # Loop through all output cubes for both data points and overplotting quantiles.
1782 # Set indexing of nplot to avoid plotting 1:1 scatter of cubes[0] vs cubes[0]
1783 for plottype in ["points", "quantiles"]:
1784 nplot = 0
1785 for cube in iter_maybe(cubes):
1786 label = None
1787 color = "black"
1788 if model_colors_map: 1788 ↛ 1793line 1788 didn't jump to line 1793 because the condition on line 1788 was always true
1789 label = cube.attributes.get("model_name")
1790 color = model_colors_map[label]
1792 # Plot all data points
1793 if plottype == "points":
1794 if nplot > 0:
1795 if hexbin:
1796 hb = plt.hexbin(
1797 cubes[0].data.flatten(),
1798 cube.data.flatten(),
1799 alpha=0.3,
1800 gridsize=100,
1801 mincnt=1,
1802 )
1803 else:
1804 plt.scatter(
1805 cubes[0].data.flatten(),
1806 cube.data.flatten(),
1807 color=color,
1808 marker="+",
1809 label=None,
1810 alpha=0.3,
1811 )
1813 elif plottype == "quantiles": 1813 ↛ 1832line 1813 didn't jump to line 1832 because the condition on line 1813 was always true
1814 # Construct Q-Q plot
1815 quantiles.append(
1816 cube.collapsed(
1817 cube.coords(dim_coords=True),
1818 iris.analysis.PERCENTILE,
1819 percent=percentiles,
1820 )
1821 )
1822 if nplot > 0:
1823 iplt.scatter(
1824 quantiles[0],
1825 quantiles[-1],
1826 color=color,
1827 marker="o",
1828 label=label,
1829 edgecolors="black",
1830 )
1832 nplot = nplot + 1
1834 # Add some labels and tweak the style.
1835 ax.set_title(title, fontsize=16)
1836 ax.set_xlabel(
1837 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1838 )
1839 ax.set_ylabel(
1840 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1841 )
1842 ax.tick_params(axis="both", labelsize=12)
1843 ax.autoscale()
1845 # Set 1:1 line and equal axes if scatter plot of common cube names
1846 nameA = iter_maybe(cubes)[0].name()
1847 nameB = iter_maybe(cubes)[1].name()
1848 if any(part in nameB.split("_") for part in nameA.split("_")): 1848 ↛ 1859line 1848 didn't jump to line 1859 because the condition on line 1848 was always true
1849 lims = [
1850 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1851 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1852 ]
1853 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1854 ax.set_aspect("equal")
1855 ax.set_xlim(lims)
1856 ax.set_ylim(lims)
1858 # Overlay grid-lines onto scatter plot.
1859 ax.grid(linestyle="--", color="grey", linewidth=1)
1860 if model_colors_map: 1860 ↛ 1864line 1860 didn't jump to line 1864 because the condition on line 1860 was always true
1861 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1863 # Add colorbar if hexbin output
1864 if hexbin:
1865 cb = plt.colorbar(
1866 hb, orientation="horizontal", location="bottom", pad=0.08, shrink=0.7
1867 )
1868 cb.set_label("Number of data points", size=12)
1870 # Save plot.
1871 _save_close_figure(fig, "scatter", filename)
1874def _spatial_plot(
1875 method: Literal["contourf", "pcolormesh", "scatter"],
1876 cube: iris.cube.Cube,
1877 filename: str | None,
1878 sequence_coordinate: str,
1879 stamp_coordinate: str,
1880 overlay_cube: iris.cube.Cube | None = None,
1881 contour_cube: iris.cube.Cube | None = None,
1882 point_cube: iris.cube.Cube | None = None,
1883 **kwargs,
1884):
1885 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1887 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1888 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1889 is present then postage stamp plots will be produced.
1891 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1892 be overplotted on the same figure.
1894 Parameters
1895 ----------
1896 method: "contourf" | "pcolormesh" | "scatter"
1897 The plotting method to use.
1898 Select choice of "contourf" or "pcolormesh" for gridded data.
1899 Use "scatter" for point-based data.
1900 cube: Cube
1901 Iris cube of the data to plot. It should have two spatial dimensions,
1902 such as lat and lon, and may also have a another two dimension to be
1903 plotted sequentially and/or as postage stamp plots.
1904 filename: str | None
1905 Name of the plot to write, used as a prefix for plot sequences. If None
1906 uses the recipe name.
1907 sequence_coordinate: str
1908 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1909 This coordinate must exist in the cube.
1910 stamp_coordinate: str
1911 Coordinate about which to plot postage stamp plots. Defaults to
1912 ``"realization"``.
1913 overlay_cube: Cube | None, optional
1914 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1915 contour_cube: Cube | None, optional
1916 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1917 point_cube: Cube | None, optional
1918 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
1920 Raises
1921 ------
1922 ValueError
1923 If the cube doesn't have the right dimensions.
1924 TypeError
1925 If the cube isn't a single cube.
1926 """
1927 # Ensure we've got a single cube.
1928 cube = check_single_cube(cube)
1930 # Set title based on recipe metadata or use cube name
1931 recipe_title = get_recipe_metadata().get("title", cube.name())
1933 # Check if there is a valid stamp coordinate in cube dimensions.
1934 if stamp_coordinate == "realization": 1934 ↛ 1939line 1934 didn't jump to line 1939 because the condition on line 1934 was always true
1935 stamp_coordinate = check_stamp_coordinate(cube)
1937 # Make postage stamp plots if stamp_coordinate exists and has more than a
1938 # single point.
1939 plotting_func = _plot_and_save_spatial_plot
1940 try:
1941 if cube.coord(stamp_coordinate).shape[0] > 1:
1942 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1943 except iris.exceptions.CoordinateNotFoundError:
1944 pass
1946 # Produce a geographical scatter plot if the data have a
1947 # dimension called observation or model_obs_error
1948 if any(
1949 crd.var_name == "station"
1950 or crd.var_name == "Station_Name"
1951 or crd.var_name == "model_obs_error"
1952 for crd in cube.coords()
1953 ):
1954 plotting_func = _plot_and_save_spatial_plot
1955 method = "scatter"
1957 # Must have a sequence coordinate.
1958 try:
1959 cube.coord(sequence_coordinate)
1960 except iris.exceptions.CoordinateNotFoundError as err:
1961 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1963 # Create a plot for each value of the sequence coordinate.
1964 plot_index = []
1965 nplot = np.size(cube.coord(sequence_coordinate).points)
1967 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1968 # Set plot titles and filename
1969 seq_coord = cube_slice.coord(sequence_coordinate)
1971 if "model_name" in cube.attributes: 1971 ↛ 1972line 1971 didn't jump to line 1972 because the condition on line 1971 was never true
1972 model_name = cube.attributes["model_name"]
1973 else:
1974 model_name = None
1976 plot_title, plot_filename = _set_title_and_filename(
1977 seq_coord, nplot, recipe_title, filename, model_name=model_name
1978 )
1980 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1981 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1982 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1983 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1985 # Do the actual plotting.
1986 plotting_func(
1987 cube_slice,
1988 filename=plot_filename,
1989 stamp_coordinate=stamp_coordinate,
1990 title=plot_title,
1991 method=method,
1992 overlay_cube=overlay_slice,
1993 contour_cube=contour_slice,
1994 point_cube=point_slice,
1995 **kwargs,
1996 )
1997 plot_index.append(plot_filename)
1999 # Add list of plots to plot metadata.
2000 complete_plot_index = _append_to_plot_index(plot_index)
2002 # Make a page to display the plots.
2003 _make_plot_html_page(complete_plot_index)
2006####################
2007# Public functions #
2008####################
2011def spatial_contour_plot(
2012 cube: iris.cube.Cube,
2013 filename: str | None = None,
2014 sequence_coordinate: str = "time",
2015 stamp_coordinate: str = "realization",
2016 **kwargs,
2017) -> iris.cube.Cube:
2018 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2020 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2021 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2022 is present then postage stamp plots will be produced.
2024 Parameters
2025 ----------
2026 cube: Cube
2027 Iris cube of the data to plot. It should have two spatial dimensions,
2028 such as lat and lon, and may also have a another two dimension to be
2029 plotted sequentially and/or as postage stamp plots.
2030 filename: str, optional
2031 Name of the plot to write, used as a prefix for plot sequences. Defaults
2032 to the recipe name.
2033 sequence_coordinate: str, optional
2034 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2035 This coordinate must exist in the cube.
2036 stamp_coordinate: str, optional
2037 Coordinate about which to plot postage stamp plots. Defaults to
2038 ``"realization"``.
2040 Returns
2041 -------
2042 Cube
2043 The original cube (so further operations can be applied).
2045 Raises
2046 ------
2047 ValueError
2048 If the cube doesn't have the right dimensions.
2049 TypeError
2050 If the cube isn't a single cube.
2051 """
2052 _spatial_plot(
2053 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2054 )
2055 return cube
2058def spatial_pcolormesh_plot(
2059 cubes: iris.cube.Cube | iris.cube.CubeList,
2060 filename: str | None = None,
2061 sequence_coordinate: str = "time",
2062 stamp_coordinate: str = "realization",
2063 **kwargs,
2064) -> iris.cube.Cube:
2065 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2067 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2068 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2069 is present then postage stamp plots will be produced.
2071 This function is significantly faster than ``spatial_contour_plot``,
2072 especially at high resolutions, and should be preferred unless contiguous
2073 contour areas are important.
2075 Parameters
2076 ----------
2077 cube: Cubes
2078 Iris cube or cubelist of the data to plot. Each cube should have two spatial dimensions,
2079 such as lat and lon, and may also have a another two dimension to be
2080 plotted sequentially and/or as postage stamp plots.
2081 filename: str, optional
2082 Name of the plot to write, used as a prefix for plot sequences. Defaults
2083 to the recipe name.
2084 sequence_coordinate: str, optional
2085 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2086 This coordinate must exist in the cube.
2087 stamp_coordinate: str, optional
2088 Coordinate about which to plot postage stamp plots. Defaults to
2089 ``"realization"``.
2091 Returns
2092 -------
2093 Cubes
2094 The original cube/cubelist (so further operations can be applied).
2096 Raises
2097 ------
2098 ValueError
2099 If the cube doesn't have the right dimensions.
2100 """
2101 if isinstance(cubes, iris.cube.CubeList): 2101 ↛ 2102line 2101 didn't jump to line 2102 because the condition on line 2101 was never true
2102 for model_cube in cubes:
2103 _spatial_plot(
2104 "pcolormesh",
2105 model_cube,
2106 filename,
2107 sequence_coordinate,
2108 stamp_coordinate,
2109 **kwargs,
2110 )
2111 elif isinstance(cubes, iris.cube.Cube): 2111 ↛ 2120line 2111 didn't jump to line 2120 because the condition on line 2111 was always true
2112 _spatial_plot(
2113 "pcolormesh",
2114 cubes,
2115 filename,
2116 sequence_coordinate,
2117 stamp_coordinate,
2118 **kwargs,
2119 )
2120 return cubes
2123def spatial_multi_pcolormesh_plot(
2124 cube: iris.cube.Cube,
2125 overlay_cube: iris.cube.Cube | None = None,
2126 contour_cube: iris.cube.Cube | None = None,
2127 point_cube: iris.cube.Cube | None = None,
2128 filename: str | None = None,
2129 sequence_coordinate: str = "time",
2130 stamp_coordinate: str = "realization",
2131 **kwargs,
2132) -> iris.cube.Cube:
2133 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2135 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2136 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2137 is present then postage stamp plots will be produced.
2139 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2141 If specified, contours of a contour_cube can be overplotted on top of those.
2143 If specified, a spatial scatter map of point_cube can be overplotted.
2145 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2147 This function is significantly faster than ``spatial_contour_plot``,
2148 especially at high resolutions, and should be preferred unless contiguous
2149 contour areas are important.
2151 Parameters
2152 ----------
2153 cube: Cube
2154 Iris cube of the data to plot. It should have two spatial dimensions,
2155 such as lat and lon, and may also have two additional dimensions to be
2156 plotted sequentially and/or as postage stamp plots.
2157 overlay_cube: Cube, optional
2158 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2159 such as lat and lon, and may also have two additional dimensions to be
2160 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.
2161 If not provided, output plot generated without overlay cube.
2162 contour_cube: Cube, optional
2163 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,
2164 such as lat and lon, and may also have two additional dimensions to be
2165 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2166 point_cube: Cube, optional
2167 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
2168 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2169 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2170 postage stamp plots. If not provided, output plot generated without point-based layer.
2171 filename: str, optional
2172 Name of the plot to write, used as a prefix for plot sequences. Defaults
2173 to the recipe name.
2174 sequence_coordinate: str, optional
2175 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2176 This coordinate must exist in the cube.
2177 stamp_coordinate: str, optional
2178 Coordinate about which to plot postage stamp plots. Defaults to
2179 ``"realization"``.
2181 Returns
2182 -------
2183 Cube
2184 The original cube (so further operations can be applied).
2186 Raises
2187 ------
2188 ValueError
2189 If the cube doesn't have the right dimensions.
2190 TypeError
2191 If the cube isn't a single cube.
2192 """
2193 _spatial_plot(
2194 "pcolormesh",
2195 cube,
2196 filename,
2197 sequence_coordinate,
2198 stamp_coordinate,
2199 overlay_cube=overlay_cube,
2200 contour_cube=contour_cube,
2201 point_cube=point_cube,
2202 )
2203 return cube, overlay_cube, contour_cube, point_cube
2206# TODO: Expand function to handle ensemble data.
2207# line_coordinate: str, optional
2208# Coordinate about which to plot multiple lines. Defaults to
2209# ``"realization"``.
2210def plot_line_series(
2211 cube: iris.cube.Cube | iris.cube.CubeList,
2212 filename: str | None = None,
2213 series_coordinate: str = "time",
2214 sequence_coordinate: str = "time",
2215 # add the following for ensembles
2216 stamp_coordinate: str = "realization",
2217 single_plot: bool = False,
2218 **kwargs,
2219) -> iris.cube.Cube | iris.cube.CubeList:
2220 """Plot a line plot for the specified coordinate.
2222 The Cube or CubeList must be 1D.
2224 Parameters
2225 ----------
2226 iris.cube | iris.cube.CubeList
2227 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2228 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2229 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2230 filename: str, optional
2231 Name of the plot to write, used as a prefix for plot sequences. Defaults
2232 to the recipe name.
2233 series_coordinate: str, optional
2234 Coordinate about which to make a series. Defaults to ``"time"``. This
2235 coordinate must exist in the cube.
2237 Returns
2238 -------
2239 iris.cube.Cube | iris.cube.CubeList
2240 The original Cube or CubeList (so further operations can be applied).
2242 Raises
2243 ------
2244 ValueError
2245 If the cubes don't have the right dimensions.
2246 TypeError
2247 If the cube isn't a Cube or CubeList.
2248 """
2249 # Ensure we have a name for the plot file.
2250 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2252 num_models = get_num_models(cube)
2254 validate_cube_shape(cube, num_models)
2256 # Iterate over all cubes and extract coordinate to plot.
2257 cubes = iris.cube.CubeList(iter_maybe(cube))
2259 coords = []
2260 for model_cube in cubes:
2261 try:
2262 coords.append(model_cube.coord(series_coordinate))
2263 except iris.exceptions.CoordinateNotFoundError as err:
2264 raise ValueError(
2265 f"Cube must have a {series_coordinate} coordinate."
2266 ) from err
2267 # Count cube dimensions and exclude realization and
2268 # forecast_reference_time if they exist.
2269 ndim = model_cube.ndim
2271 if model_cube.coords("realization"): 2271 ↛ 2279line 2271 didn't jump to line 2279 because the condition on line 2271 was always true
2272 # returns coord dimension
2273 realization_dims = model_cube.coord_dims("realization")
2275 # Only subtract if realization is a dimension coordinate
2276 if realization_dims:
2277 ndim -= len(realization_dims)
2279 if model_cube.coords("forecast_reference_time"):
2280 frt_dims = model_cube.coord_dims("forecast_reference_time")
2282 # Only subtract if frt is a dimension coordinate
2283 if frt_dims: 2283 ↛ 2284line 2283 didn't jump to line 2284 because the condition on line 2283 was never true
2284 ndim -= len(frt_dims)
2286 if ndim > 2:
2287 raise ValueError(
2288 "Cube must be 1D or 2D (excluding any realization or forecast_reference_time dimensions)."
2289 )
2291 plot_index = []
2293 # Check if this is a spectral plot by looking for spectral coordinates
2294 is_spectral_plot = series_coordinate in [
2295 "frequency",
2296 "physical_wavenumber",
2297 "wavelength",
2298 ]
2300 if is_spectral_plot:
2301 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2302 # coordinate frequency/wavenumber.
2303 # If several power spectra are plotted with time as sequence_coordinate for the
2304 # time slider option.
2306 # Internal plotting function.
2307 plotting_func = _plot_and_save_line_power_spectrum_series
2309 for model_cube in cubes:
2310 try:
2311 model_cube.coord(sequence_coordinate)
2312 except iris.exceptions.CoordinateNotFoundError as err:
2313 raise ValueError(
2314 f"Cube must have a {sequence_coordinate} coordinate."
2315 ) from err
2317 if num_models == 1: 2317 ↛ 2332line 2317 didn't jump to line 2332 because the condition on line 2317 was always true
2318 # check for ensembles
2319 if ( 2319 ↛ 2323line 2319 didn't jump to line 2323 because the condition on line 2319 was never true
2320 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2321 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2322 ):
2323 if single_plot:
2324 # Plot spectra, mean and ensemble spread on 1 plot
2325 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2326 else:
2327 # Plot postage stamps
2328 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2329 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2330 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2331 else:
2332 all_points = sorted(
2333 set(
2334 itertools.chain.from_iterable(
2335 cb.coord(sequence_coordinate).points for cb in cubes
2336 )
2337 )
2338 )
2339 all_slices = list(
2340 itertools.chain.from_iterable(
2341 cb.slices_over(sequence_coordinate) for cb in cubes
2342 )
2343 )
2344 # Matched slices (matched by seq coord point; it may happen that
2345 # evaluated models do not cover the same seq coord range, hence matching
2346 # necessary)
2347 cube_iterables = [
2348 iris.cube.CubeList(
2349 s
2350 for s in all_slices
2351 if s.coord(sequence_coordinate).points[0] == point
2352 )
2353 for point in all_points
2354 ]
2355 nplot = len(all_points)
2357 # Create a plot for each value of the sequence coordinate. Allowing for
2358 # multiple cubes in a CubeList to be plotted in the same plot for similar
2359 # sequence values. Passing a CubeList into the internal plotting function
2360 # for similar values of the sequence coordinate. cube_slice can be an
2361 # iris.cube.Cube or an iris.cube.CubeList.
2363 for cube_slice in cube_iterables:
2364 # Normalize cube_slice to a list of cubes
2365 if isinstance(cube_slice, iris.cube.CubeList): 2365 ↛ 2366line 2365 didn't jump to line 2366 because the condition on line 2365 was never true
2366 cubes = list(cube_slice)
2367 elif isinstance(cube_slice, iris.cube.Cube): 2367 ↛ 2370line 2367 didn't jump to line 2370 because the condition on line 2367 was always true
2368 cubes = [cube_slice]
2369 else:
2370 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2372 # Use sequence value so multiple sequences can merge.
2373 seq_coord = cube_slice[0].coord(sequence_coordinate)
2374 plot_title, plot_filename = _set_title_and_filename(
2375 seq_coord, nplot, recipe_title, filename
2376 )
2378 # Format the coordinate value in a unit appropriate way.
2379 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2381 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2382 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2382 ↛ 2383line 2382 didn't jump to line 2383 because the condition on line 2382 was never true
2383 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2385 # Do the actual plotting.
2386 plotting_func(
2387 cube_slice,
2388 coords,
2389 stamp_coordinate,
2390 plot_filename,
2391 title,
2392 series_coordinate,
2393 )
2395 plot_index.append(plot_filename)
2396 else:
2397 # Format the title and filename using plotted series coordinate
2398 nplot = 1
2399 seq_coord = coords[0]
2400 plot_title, plot_filename = _set_title_and_filename(
2401 seq_coord, nplot, recipe_title, filename
2402 )
2404 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2405 if (
2406 "station" in [c.name() for c in cubes[0].coords()]
2407 and len(cubes[0].coord("station").points) > 1
2408 ):
2409 for station in cubes[0].coord("station").points:
2410 station_cubes = cubes.extract(iris.Constraint(station=station))
2411 station_name = station_cubes[0].coord("Station_Name").points[0]
2412 station_plotname = plot_filename.replace(
2413 ".png", "_" + station_name + ".png"
2414 )
2415 _plot_and_save_line_series(
2416 station_cubes,
2417 coords,
2418 "realization",
2419 station_plotname,
2420 f"{plot_title} {station_name}",
2421 )
2422 plot_index.append(station_plotname)
2424 else:
2425 # Do the actual plotting for all other series coordinate options.
2426 _plot_and_save_line_series(
2427 cubes, coords, stamp_coordinate, plot_filename, plot_title
2428 )
2430 plot_index.append(plot_filename)
2432 # append plot to list of plots
2433 complete_plot_index = _append_to_plot_index(plot_index)
2435 # Make a page to display the plots.
2436 _make_plot_html_page(complete_plot_index)
2438 return cube
2441def plot_vertical_line_series(
2442 cubes: iris.cube.Cube | iris.cube.CubeList,
2443 filename: str | None = None,
2444 series_coordinate: str = "model_level_number",
2445 sequence_coordinate: str = "time",
2446 # line_coordinate: str = "realization",
2447 **kwargs,
2448) -> iris.cube.Cube | iris.cube.CubeList:
2449 """Plot a line plot against a type of vertical coordinate.
2451 The Cube or CubeList must be 1D.
2453 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2454 then a sequence of plots will be produced.
2456 Parameters
2457 ----------
2458 iris.cube | iris.cube.CubeList
2459 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2460 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2461 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2462 filename: str, optional
2463 Name of the plot to write, used as a prefix for plot sequences. Defaults
2464 to the recipe name.
2465 series_coordinate: str, optional
2466 Coordinate to plot on the y-axis. Can be ``pressure`` or
2467 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2468 for LFRic. Defaults to ``model_level_number``.
2469 This coordinate must exist in the cube.
2470 sequence_coordinate: str, optional
2471 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2472 This coordinate must exist in the cube.
2474 Returns
2475 -------
2476 iris.cube.Cube | iris.cube.CubeList
2477 The original Cube or CubeList (so further operations can be applied).
2478 Plotted data.
2480 Raises
2481 ------
2482 ValueError
2483 If the cubes doesn't have the right dimensions.
2484 TypeError
2485 If the cube isn't a Cube or CubeList.
2486 """
2487 # Ensure we have a name for the plot file.
2488 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2490 cubes = iter_maybe(cubes)
2491 # Initialise empty list to hold all data from all cubes in a CubeList
2492 all_data = []
2494 # Store min/max ranges for x range.
2495 x_levels = []
2497 num_models = get_num_models(cubes)
2499 validate_cube_shape(cubes, num_models)
2501 # Iterate over all cubes in cube or CubeList and plot.
2502 coords = []
2503 for cube in cubes:
2504 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2505 try:
2506 coords.append(cube.coord(series_coordinate))
2507 except iris.exceptions.CoordinateNotFoundError as err:
2508 raise ValueError(
2509 f"Cube must have a {series_coordinate} coordinate."
2510 ) from err
2512 try:
2513 if cube.ndim > 1 or not cube.coords("realization"): 2513 ↛ 2521line 2513 didn't jump to line 2521 because the condition on line 2513 was always true
2514 cube.coord(sequence_coordinate)
2515 except iris.exceptions.CoordinateNotFoundError as err:
2516 raise ValueError(
2517 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2518 ) from err
2520 # Get minimum and maximum from levels information.
2521 _, levels, _ = colorbar_map_levels(cube, axis="x")
2522 if levels is not None: 2522 ↛ 2526line 2522 didn't jump to line 2526 because the condition on line 2522 was always true
2523 x_levels.append(min(levels))
2524 x_levels.append(max(levels))
2525 else:
2526 all_data.append(cube.data)
2528 if len(x_levels) == 0: 2528 ↛ 2530line 2528 didn't jump to line 2530 because the condition on line 2528 was never true
2529 # Combine all data into a single NumPy array
2530 combined_data = np.concatenate(all_data)
2532 # Set the lower and upper limit for the x-axis to ensure all plots have
2533 # same range. This needs to read the whole cube over the range of the
2534 # sequence and if applicable postage stamp coordinate.
2535 vmin = np.floor(combined_data.min())
2536 vmax = np.ceil(combined_data.max())
2537 else:
2538 vmin = min(x_levels)
2539 vmax = max(x_levels)
2541 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2542 # a single profile directly without iterating over a sequence.
2543 sequence_coords = [
2544 cube.coord(sequence_coordinate)
2545 for cube in cubes
2546 if cube.coords(sequence_coordinate)
2547 ]
2548 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2549 np.size(coord.points) > 1 for coord in sequence_coords
2550 )
2551 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2552 np.size(coord.points) == 1 for coord in sequence_coords
2553 )
2555 plot_index = []
2556 if has_sequence_coord: 2556 ↛ 2581line 2556 didn't jump to line 2581 because the condition on line 2556 was always true
2557 # Matching the slices (matching by seq coord point; it may happen that
2558 # evaluated models do not cover the same seq coord range, hence matching
2559 # necessary)
2560 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2561 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2562 for cubes_slice in cube_iterables:
2563 # Format the coordinate value in a unit appropriate way.
2564 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2565 plot_title, plot_filename = _set_title_and_filename(
2566 seq_coord, nplot, recipe_title, filename
2567 )
2569 # Do the actual plotting.
2570 _plot_and_save_vertical_line_series(
2571 cubes_slice,
2572 coords,
2573 "realization",
2574 plot_filename,
2575 series_coordinate,
2576 title=plot_title,
2577 vmin=vmin,
2578 vmax=vmax,
2579 )
2580 plot_index.append(plot_filename)
2581 elif has_scalar_sequence_coord:
2582 # Scalar sequence coordinate (typically aggregated time bounds):
2583 # make one plot and include sequence period in title/filename.
2584 plot_title, plot_filename = _set_title_and_filename(
2585 sequence_coords[0], 1, recipe_title, filename
2586 )
2588 _plot_and_save_vertical_line_series(
2589 cubes,
2590 coords,
2591 "realization",
2592 plot_filename,
2593 series_coordinate,
2594 title=plot_title,
2595 vmin=vmin,
2596 vmax=vmax,
2597 )
2598 plot_index.append(plot_filename)
2599 else:
2600 # 1D case: no sequence coordinate, plot a single profile.
2601 plot_title = recipe_title
2602 if filename:
2603 plot_filename = filename
2604 else:
2605 plot_filename = f"{slugify(plot_title)}.png"
2607 _plot_and_save_vertical_line_series(
2608 cubes,
2609 coords,
2610 "realization",
2611 plot_filename,
2612 series_coordinate,
2613 title=plot_title,
2614 vmin=vmin,
2615 vmax=vmax,
2616 )
2617 plot_index.append(plot_filename)
2619 # Add list of plots to plot metadata.
2620 complete_plot_index = _append_to_plot_index(plot_index)
2622 # Make a page to display the plots.
2623 _make_plot_html_page(complete_plot_index)
2625 return cubes
2628def qq_plot(
2629 cubes: iris.cube.CubeList,
2630 coordinates: list[str],
2631 percentiles: list[float],
2632 model_names: list[str],
2633 filename: str | None = None,
2634 one_to_one: bool = True,
2635 **kwargs,
2636) -> iris.cube.CubeList:
2637 """Plot a Quantile-Quantile plot between two models for common time points.
2639 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2640 collapsed within the operator over all specified coordinates such as
2641 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2642 ensemble members to ensure a 1D cube (array).
2644 Parameters
2645 ----------
2646 cubes: iris.cube.CubeList
2647 Two cubes of the same variable with different models.
2648 coordinate: list[str]
2649 The list of coordinates to collapse over. This list should be
2650 every coordinate within the cube to result in a 1D cube around
2651 the percentile coordinate.
2652 percent: list[float]
2653 A list of percentiles to appear in the plot.
2654 model_names: list[str]
2655 A list of model names to appear on the axis of the plot.
2656 filename: str, optional
2657 Filename of the plot to write.
2658 one_to_one: bool, optional
2659 If True a 1:1 line is plotted; if False it is not. Default is True.
2661 Raises
2662 ------
2663 ValueError
2664 When the cubes are not compatible.
2666 Notes
2667 -----
2668 The quantile-quantile plot is a variant on the scatter plot representing
2669 two datasets by their quantiles (percentiles) for common time points.
2670 This plot does not use a theoretical distribution to compare against, but
2671 compares percentiles of two datasets. This plot does
2672 not use all raw data points, but plots the selected percentiles (quantiles) of
2673 each variable instead for the two datasets, thereby normalising the data for a
2674 direct comparison between the selected percentiles of the two dataset distributions.
2676 Quantile-quantile plots are valuable for comparing against
2677 observations and other models. Identical percentiles between the variables
2678 will lie on the one-to-one line implying the values correspond well to each
2679 other. Where there is a deviation from the one-to-one line a range of
2680 possibilities exist depending on how and where the data is shifted (e.g.,
2681 Wilks 2011 [Wilks2011]_).
2683 For distributions above the one-to-one line the distribution is left-skewed;
2684 below is right-skewed. A distinct break implies a bimodal distribution, and
2685 closer values/values further apart at the tails imply poor representation of
2686 the extremes.
2688 """
2689 # Check cubes using same functionality as the difference operator.
2690 if len(cubes) != 2:
2691 raise ValueError("cubes should contain exactly 2 cubes.")
2692 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2693 other: Cube = cubes.extract_cube(
2694 iris.Constraint(
2695 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2696 )
2697 )
2699 # Get spatial coord names.
2700 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2701 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2703 # Ensure cubes to compare are on common differencing grid.
2704 # This is triggered if either
2705 # i) latitude and longitude shapes are not the same. Note grid points
2706 # are not compared directly as these can differ through rounding
2707 # errors.
2708 # ii) or variables are known to often sit on different grid staggering
2709 # in different models (e.g. cell center vs cell edge), as is the case
2710 # for UM and LFRic comparisons.
2711 # In future greater choice of regridding method might be applied depending
2712 # on variable type. Linear regridding can in general be appropriate for smooth
2713 # variables. Care should be taken with interpretation of differences
2714 # given this dependency on regridding.
2715 if (
2716 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2717 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2718 ) or (
2719 base.long_name
2720 in [
2721 "eastward_wind_at_10m",
2722 "northward_wind_at_10m",
2723 "northward_wind_at_cell_centres",
2724 "eastward_wind_at_cell_centres",
2725 "zonal_wind_at_pressure_levels",
2726 "meridional_wind_at_pressure_levels",
2727 "potential_vorticity_at_pressure_levels",
2728 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2729 ]
2730 ):
2731 logger.debug("Linear regridding base cube to other grid to compute differences")
2732 base = regrid_onto_cube(base, other, method="Linear")
2734 # Extract just common time points.
2735 base, other = _extract_common_time_points(base, other)
2737 # Equalise attributes so we can merge.
2738 fully_equalise_attributes([base, other])
2739 logger.debug("Base: %s\nOther: %s", base, other)
2741 # Collapse cubes.
2742 base = collapse(
2743 base,
2744 coordinate=coordinates,
2745 method="PERCENTILE",
2746 additional_percent=percentiles,
2747 )
2748 other = collapse(
2749 other,
2750 coordinate=coordinates,
2751 method="PERCENTILE",
2752 additional_percent=percentiles,
2753 )
2755 # Ensure we have a name for the plot file.
2756 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2757 title = f"{recipe_title}"
2759 if filename is None:
2760 filename = slugify(recipe_title)
2762 # Add file extension.
2763 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2765 # Do the actual plotting on a scatter plot
2766 _plot_and_save_scatter_plot(
2767 base, other, plot_filename, title, one_to_one, model_names
2768 )
2770 # Add list of plots to plot metadata.
2771 plot_index = _append_to_plot_index([plot_filename])
2773 # Make a page to display the plots.
2774 _make_plot_html_page(plot_index)
2776 return iris.cube.CubeList([base, other])
2779def hinton_plot(
2780 cubes: iris.cube.CubeList, base_name: str, other_name: str, magnitude: bool = False
2781) -> None:
2782 """
2783 Plot a Hinton style triangle/scorecard plot.
2785 This plot type can be useful for summarising high level information, such as comparing
2786 how 'skillful' two models are when verified against observations for a variety of metrics,
2787 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2788 than customisable by the user as input arguments; many have been designed to automatically
2789 scale the plot depending on the number of x and y components.
2791 Parameters
2792 ----------
2793 cubes: iris.cube.CubeList
2794 A iris cubelist, containing at least two cubes of a skill metric to plot (model vs obs). Can
2795 include multiple variables, in which the plot will automatically scale for. If cubes containing
2796 a name significance_<var> exist, containing a bool array, then it will also plot whether each
2797 triangle is significant by using a thick black outline. Each cube should be 1D, with
2798 forecast_period as the only dimension.
2799 base_name: str
2800 The name of the base model to use as the control in the Hinton plot, as a string.
2801 other_name: str
2802 The name of the other model to use as the test in the Hinton plot, as a string.
2803 magnitude: bool
2804 Option bool when if True, then plot the numerical value difference in two values under each
2805 triangle.
2806 """
2807 # Ensure we have a name for the plot file.
2808 recipe_title = get_recipe_metadata().get("title", "Hinton")
2809 title = f"{recipe_title}"
2810 filename = slugify(recipe_title)
2812 # Check that all cubes only have one dimension called forecast_period
2813 for cube in cubes:
2814 if len(cube.dim_coords) > 1:
2815 raise ValueError(f"Should only have one dimension coord, {cube}")
2816 if cube.dim_coords[0].name() != "forecast_period":
2817 raise ValueError(
2818 f"Single coord should be forecast_period, not {cube.dim_coords[0].name()}"
2819 )
2821 # Separate out base cubes and other cubes.
2822 base_cubes = iris.cube.CubeList()
2823 other_cubes = iris.cube.CubeList()
2824 for c in cubes:
2825 if c.attributes["model_name"] == base_name:
2826 base_cubes.append(c)
2827 elif c.attributes["model_name"] == other_name:
2828 other_cubes.append(c)
2830 # base cubes should be the same length as other cubes, otherwise one is missing a variable.
2831 if len(base_cubes) != len(other_cubes):
2832 raise ValueError(
2833 f"base cubes {base_cubes} are not same number as {other_cubes}"
2834 )
2836 # Find common variable names in the two groups.
2837 base_vars = {cube.long_name for cube in base_cubes if cube.long_name is not None}
2838 other_vars = {cube.long_name for cube in other_cubes if cube.long_name is not None}
2839 common_vars = sorted(base_vars & other_vars)
2841 # Iterate over each variable (row)
2842 rows = []
2843 for var in common_vars:
2844 # Extract cube with matching variable name
2845 base_cube = next(
2846 (c for c in base_cubes if c.long_name == var),
2847 None,
2848 )
2850 other_cube = next(
2851 (c for c in other_cubes if c.long_name == var),
2852 None,
2853 )
2855 # If we can't find a variable in both cubes, then skip
2856 if base_cube is None or other_cube is None: 2856 ↛ 2857line 2856 didn't jump to line 2857 because the condition on line 2856 was never true
2857 continue
2859 # Compute difference (1D array)
2860 # We can already make assumption both on same forecast_periods as checked
2861 # prior to computing metric.
2862 diff = other_cube.data - base_cube.data
2864 # See if there is a significance cube present, if not, set as None.
2865 sig_cube = next(
2866 (cube for cube in cubes if cube.long_name == f"significance_{var}"),
2867 None,
2868 )
2870 # Append row information.
2871 rows.append(
2872 {
2873 "name": var,
2874 "forecast_periods": base_cube.coord("forecast_period").points,
2875 "change": diff,
2876 "significance": sig_cube.data.astype(bool)
2877 if sig_cube is not None
2878 else None,
2879 }
2880 )
2882 # For each row, compute standardised anomalies
2883 for row in rows:
2884 change = np.asarray(row["change"])
2886 anoms = change - np.mean(change)
2888 scale = np.max(np.abs(anoms))
2890 if scale > 0:
2891 scaled = anoms / scale
2892 else:
2893 scaled = np.zeros_like(anoms)
2895 row["anoms"] = anoms
2896 row["scaled"] = scaled
2898 # Setup colors of triangles
2899 color_pos = "#7CAE00"
2900 color_neg = "#7B68EE"
2902 # Setup cell/text size ratios
2903 figsize = None
2904 cell_size_in = 1.5
2905 text_row_ratio = 0.25
2907 # Get the number of x and y elements
2908 ny = len(rows)
2909 nx = max(len(row["forecast_periods"]) for row in rows)
2911 # Build non-uniform y coordinates
2912 tri_height = 1.0
2913 txt_height = text_row_ratio
2915 tri_y = []
2916 txt_y = []
2917 y_edges = [0.0]
2919 y = 0.0
2920 for _j in range(ny):
2921 tri_y.append(y + tri_height / 2)
2922 y += tri_height
2923 y_edges.append(y)
2925 if magnitude:
2926 txt_y.append(y + txt_height / 2)
2927 y += txt_height
2928 y_edges.append(y)
2930 total_height = y
2932 # Dynamic figure size
2933 if figsize is None: 2933 ↛ 2938line 2933 didn't jump to line 2938 because the condition on line 2933 was always true
2934 width = nx * cell_size_in
2935 height = total_height * cell_size_in + 2
2936 figsize = (width, height)
2938 fig, ax = plt.subplots(figsize=figsize)
2940 # Setup axes and grid.
2941 ax.set_aspect("equal", adjustable="box")
2942 ax.set_xlim(-0.5, nx - 0.5)
2943 ax.set_ylim(0, total_height)
2945 longest_row = max(rows, key=lambda row: len(row["forecast_periods"]))
2947 ax.set_xticks(np.arange(nx))
2948 ax.set_xticklabels(
2949 longest_row["forecast_periods"],
2950 rotation=90,
2951 )
2953 ax.set_yticks(tri_y)
2954 ax.set_yticklabels([row["name"] for row in rows])
2956 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2957 ax.set_yticks(y_edges, minor=True)
2959 ax.set_axisbelow(True)
2960 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2961 ax.grid(False, which="major")
2962 ax.tick_params(which="minor", length=0)
2964 ax.invert_yaxis()
2966 # Compute marker scaling (fixed overlap)
2967 fig.canvas.draw()
2969 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2970 width_in, height_in = bbox.width, bbox.height
2972 cell_w = (width_in * fig.dpi) / nx
2973 cell_h = (height_in * fig.dpi) / total_height
2974 cell_pixels = min(cell_w, cell_h)
2976 max_marker_size = (0.6 * cell_pixels) ** 2
2978 text_fontsize = cell_pixels * 0.15
2980 # Plot triangles + text
2981 for j, row in enumerate(rows):
2982 scaled = row["scaled"]
2983 anoms = row["anoms"]
2984 signif = row["significance"]
2986 for i in range(len(scaled)):
2987 val = scaled[i]
2989 if np.isnan(val): 2989 ↛ 2990line 2989 didn't jump to line 2990 because the condition on line 2989 was never true
2990 continue
2992 if abs(val) < 0.01:
2993 continue
2995 size = max_marker_size * abs(val)
2997 if val >= 0:
2998 marker = "^"
2999 color = color_pos
3000 else:
3001 marker = "v"
3002 color = color_neg
3004 if signif is not None:
3005 sig = bool(signif[i])
3006 edgecolor = "black" if sig else "none"
3007 linewidth = 0.6 if sig else 0.0
3008 else:
3009 edgecolor = "none"
3010 linewidth = 0.0
3012 ax.scatter(
3013 i,
3014 tri_y[j],
3015 s=size,
3016 marker=marker,
3017 c=color,
3018 edgecolors=edgecolor,
3019 linewidths=linewidth,
3020 zorder=3,
3021 clip_on=True,
3022 )
3024 # Text row
3025 if magnitude:
3026 mag_val = anoms[i]
3028 if not np.isnan(mag_val): 3028 ↛ 2986line 3028 didn't jump to line 2986 because the condition on line 3028 was always true
3029 ax.text(
3030 i,
3031 txt_y[j],
3032 f"{mag_val:.1f}",
3033 ha="center",
3034 va="center",
3035 fontsize=text_fontsize,
3036 color="black",
3037 zorder=4,
3038 )
3040 ax.set_title(title)
3041 plt.tight_layout()
3043 # Save plot.
3044 _save_close_figure(fig, "hinton", filename)
3046 # Add file extension.
3047 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3049 # Add list of plots to plot metadata.
3050 plot_index = _append_to_plot_index([plot_filename])
3052 # Make a page to display the plots.
3053 _make_plot_html_page(plot_index)
3056def scatter_plot(
3057 cube_x: iris.cube.Cube | iris.cube.CubeList,
3058 cube_y: iris.cube.Cube | iris.cube.CubeList,
3059 filename: str | None = None,
3060 one_to_one: bool = True,
3061 **kwargs,
3062) -> iris.cube.CubeList:
3063 """Plot a scatter plot between two variables.
3065 Both cubes must be 1D.
3067 Parameters
3068 ----------
3069 cube_x: Cube | CubeList
3070 1 dimensional Cube of the data to plot on y-axis.
3071 cube_y: Cube | CubeList
3072 1 dimensional Cube of the data to plot on x-axis.
3073 filename: str, optional
3074 Filename of the plot to write.
3075 one_to_one: bool, optional
3076 If True a 1:1 line is plotted; if False it is not. Default is True.
3078 Returns
3079 -------
3080 cubes: CubeList
3081 CubeList of the original x and y cubes for further processing.
3083 Raises
3084 ------
3085 ValueError
3086 If the cube doesn't have the right dimensions and cubes not the same
3087 size.
3088 TypeError
3089 If the cube isn't a single cube.
3091 Notes
3092 -----
3093 Scatter plots are used for determining if there is a relationship between
3094 two variables. Positive relations have a slope going from bottom left to top
3095 right; Negative relations have a slope going from top left to bottom right.
3096 """
3097 # Iterate over all cubes in cube or CubeList and plot.
3098 for cube_iter in iter_maybe(cube_x):
3099 # Check cubes are correct shape.
3100 cube_iter = check_single_cube(cube_iter)
3101 if cube_iter.ndim > 1:
3102 raise ValueError("cube_x must be 1D.")
3104 # Iterate over all cubes in cube or CubeList and plot.
3105 for cube_iter in iter_maybe(cube_y):
3106 # Check cubes are correct shape.
3107 cube_iter = check_single_cube(cube_iter)
3108 if cube_iter.ndim > 1:
3109 raise ValueError("cube_y must be 1D.")
3111 # Ensure we have a name for the plot file.
3112 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
3113 title = f"{recipe_title}"
3115 if filename is None:
3116 filename = slugify(recipe_title)
3118 # Add file extension.
3119 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3121 # Do the actual plotting.
3122 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3124 # Add list of plots to plot metadata.
3125 plot_index = _append_to_plot_index([plot_filename])
3127 # Make a page to display the plots.
3128 _make_plot_html_page(plot_index)
3130 return iris.cube.CubeList([cube_x, cube_y])
3133def vector_plot(
3134 cube_u: iris.cube.Cube,
3135 cube_v: iris.cube.Cube,
3136 filename: str | None = None,
3137 sequence_coordinate: str = "time",
3138 **kwargs,
3139) -> iris.cube.CubeList:
3140 """Plot a vector plot based on the input u and v components."""
3141 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
3143 # Cubes must have a matching sequence coordinate.
3144 try:
3145 # Check that the u and v cubes have the same sequence coordinate.
3146 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3146 ↛ anywhereline 3146 didn't jump anywhere: it always raised an exception.
3147 raise ValueError("Coordinates do not match.")
3148 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3149 raise ValueError(
3150 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3151 ) from err
3153 # Create a plot for each value of the sequence coordinate.
3154 plot_index = []
3155 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3156 for cube_u_slice, cube_v_slice in zip(
3157 cube_u.slices_over(sequence_coordinate),
3158 cube_v.slices_over(sequence_coordinate),
3159 strict=True,
3160 ):
3161 # Format the coordinate value in a unit appropriate way.
3162 seq_coord = cube_u_slice.coord(sequence_coordinate)
3163 plot_title, plot_filename = _set_title_and_filename(
3164 seq_coord, nplot, recipe_title, filename
3165 )
3167 # Do the actual plotting.
3168 _plot_and_save_vector_plot(
3169 cube_u_slice,
3170 cube_v_slice,
3171 filename=plot_filename,
3172 title=plot_title,
3173 method="pcolormesh",
3174 )
3175 plot_index.append(plot_filename)
3177 # Add list of plots to plot metadata.
3178 complete_plot_index = _append_to_plot_index(plot_index)
3180 # Make a page to display the plots.
3181 _make_plot_html_page(complete_plot_index)
3183 return iris.cube.CubeList([cube_u, cube_v])
3186def plot_histogram_series(
3187 cubes: iris.cube.Cube | iris.cube.CubeList,
3188 filename: str | None = None,
3189 sequence_coordinate: str = "time",
3190 stamp_coordinate: str = "realization",
3191 single_plot: bool = False,
3192 **kwargs,
3193) -> iris.cube.Cube | iris.cube.CubeList:
3194 """Plot a histogram plot for each vertical level provided.
3196 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3197 is present then a sequence of plots will be produced using the time slider
3198 functionality to scroll through histograms against time. If a
3199 stamp_coordinate is present then postage stamp plots will be produced. If
3200 stamp_coordinate and single_plot is True, all postage stamp plots will be
3201 plotted in a single plot instead of separate postage stamp plots.
3203 Parameters
3204 ----------
3205 cubes: Cube | iris.cube.CubeList
3206 Iris cube or CubeList of the data to plot. It should have a single dimension other
3207 than the stamp coordinate.
3208 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3209 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3210 filename: str, optional
3211 Name of the plot to write, used as a prefix for plot sequences. Defaults
3212 to the recipe name.
3213 sequence_coordinate: str, optional
3214 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3215 This coordinate must exist in the cube and will be used for the time
3216 slider.
3217 stamp_coordinate: str, optional
3218 Coordinate about which to plot postage stamp plots. Defaults to
3219 ``"realization"``.
3220 single_plot: bool, optional
3221 If True, all postage stamp plots will be plotted in a single plot. If
3222 False, each postage stamp plot will be plotted separately. Is only valid
3223 if stamp_coordinate exists and has more than a single point.
3225 Returns
3226 -------
3227 iris.cube.Cube | iris.cube.CubeList
3228 The original Cube or CubeList (so further operations can be applied).
3229 Plotted data.
3231 Raises
3232 ------
3233 ValueError
3234 If the cube doesn't have the right dimensions.
3235 TypeError
3236 If the cube isn't a Cube or CubeList.
3237 """
3238 recipe_title = get_recipe_metadata().get("title", "Histogram")
3240 cubes = iter_maybe(cubes)
3242 # Internal plotting function.
3243 plotting_func = _plot_and_save_histogram_series
3245 num_models = get_num_models(cubes)
3247 validate_cube_shape(cubes, num_models)
3249 # If several histograms are plotted, check sequence_coordinate
3250 check_sequence_coordinate(cubes, sequence_coordinate)
3252 # Get axis minimum and maximum from levels information.
3253 # If no levels set, derive minima and maxima from data in CubeList.
3254 vmin, vmax = _set_axis_range(cubes)
3256 # Make postage stamp plots if stamp_coordinate exists and has more than a
3257 # single point. If single_plot is True:
3258 # -- all postage stamp plots will be plotted in a single plot instead of
3259 # separate postage stamp plots.
3260 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3261 # produced per single model only
3262 if num_models == 1:
3263 if ( 3263 ↛ 3267line 3263 didn't jump to line 3267 because the condition on line 3263 was never true
3264 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3265 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3266 ):
3267 if single_plot:
3268 plotting_func = (
3269 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3270 )
3271 else:
3272 plotting_func = _plot_and_save_postage_stamp_histogram_series
3273 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3274 else:
3275 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3277 plot_index = []
3278 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3279 # Create a plot for each value of the sequence coordinate. Allowing for
3280 # multiple cubes in a CubeList to be plotted in the same plot for similar
3281 # sequence values. Passing a CubeList into the internal plotting function
3282 # for similar values of the sequence coordinate. cube_slice can be an
3283 # iris.cube.Cube or an iris.cube.CubeList.
3284 for cube_slice in cube_iterables:
3285 single_cube = cube_slice
3286 if isinstance(cube_slice, iris.cube.CubeList):
3287 single_cube = cube_slice[0]
3289 # Ensure valid stamp coordinate in cube dimensions
3290 if stamp_coordinate == "realization": 3290 ↛ 3293line 3290 didn't jump to line 3293 because the condition on line 3290 was always true
3291 stamp_coordinate = check_stamp_coordinate(single_cube)
3292 # Set plot titles and filename, based on sequence coordinate
3293 seq_coord = single_cube.coord(sequence_coordinate)
3294 # Use time coordinate in title and filename if single histogram output.
3295 if sequence_coordinate == "realization" and nplot == 1: 3295 ↛ 3296line 3295 didn't jump to line 3296 because the condition on line 3295 was never true
3296 seq_coord = single_cube.coord("time")
3297 # Use station name in title and filename if model vs obs comparison
3298 if sequence_coordinate == "station": 3298 ↛ 3299line 3298 didn't jump to line 3299 because the condition on line 3298 was never true
3299 seq_coord = single_cube.coord("Station_Name")
3301 plot_title, plot_filename = _set_title_and_filename(
3302 seq_coord, nplot, recipe_title, filename
3303 )
3305 # Do the actual plotting.
3306 plotting_func(
3307 cube_slice,
3308 filename=plot_filename,
3309 stamp_coordinate=stamp_coordinate,
3310 title=plot_title,
3311 vmin=vmin,
3312 vmax=vmax,
3313 )
3314 plot_index.append(plot_filename)
3316 # Add list of plots to plot metadata.
3317 complete_plot_index = _append_to_plot_index(plot_index)
3319 # Make a page to display the plots.
3320 _make_plot_html_page(complete_plot_index)
3322 return cubes
3325def plot_scatter_series(
3326 cubes: iris.cube.Cube | iris.cube.CubeList,
3327 filename: str | None = None,
3328 sequence_coordinate: str = "time",
3329 stamp_coordinate: str = "realization",
3330 hexbin: bool = False,
3331 **kwargs,
3332) -> iris.cube.Cube | iris.cube.CubeList:
3333 """Plot a scatter plot for each sequence coordinate provided.
3335 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3336 is present then a sequence of plots will be produced using the time slider
3337 functionality to scroll through scatter against time. If a
3338 stamp_coordinate is present then postage stamp plots will be produced. If
3339 stamp_coordinate and single_plot is True, all postage stamp plots will be
3340 plotted in a single plot instead of separate postage stamp plots.
3342 Parameters
3343 ----------
3344 cubes: Cube | iris.cube.CubeList
3345 Iris cube or CubeList of the data to plot. It should have a single dimension other
3346 than the stamp coordinate.
3347 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3348 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3349 filename: str, optional
3350 Name of the plot to write, used as a prefix for plot sequences. Defaults
3351 to the recipe name.
3352 sequence_coordinate: str, optional
3353 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3354 This coordinate must exist in the cube and will be used for the time
3355 slider.
3356 stamp_coordinate: str, optional
3357 Coordinate about which to plot postage stamp plots. Defaults to
3358 ``"realization"``.
3359 hexbin: bool, optional
3360 If True, generate hexbin comparison plot.
3361 If False, generate point-by-point scatter plot.
3363 Returns
3364 -------
3365 iris.cube.Cube | iris.cube.CubeList
3366 The original Cube or CubeList (so further operations can be applied).
3367 Plotted data.
3369 Raises
3370 ------
3371 ValueError
3372 If the cube doesn't have the right dimensions.
3373 TypeError
3374 If the cube isn't a Cube or CubeList.
3375 """
3376 recipe_title = get_recipe_metadata().get("title", "Scatter")
3378 cubes = iter_maybe(cubes)
3380 # Internal plotting function.
3381 plotting_func = _plot_and_save_scatter_series
3383 num_models = get_num_models(cubes)
3385 validate_cube_shape(cubes, num_models)
3387 check_sequence_coordinate(cubes, sequence_coordinate)
3389 vmin, vmax = _set_axis_range(cubes)
3391 # Require >1 models to compare on scatter plot
3392 if num_models > 1:
3393 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3394 else:
3395 raise ValueError(
3396 "Scatter plot series requires multiple number of models in input data."
3397 )
3399 plot_index = []
3400 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3401 # Create a plot for each value of the sequence coordinate. Allowing for
3402 # multiple cubes in a CubeList to be plotted in the same plot for similar
3403 # sequence values. Passing a CubeList into the internal plotting function
3404 # for similar values of the sequence coordinate. cube_slice can be an
3405 # iris.cube.Cube or an iris.cube.CubeList.
3406 for cube_slice in cube_iterables:
3407 single_cube = cube_slice
3408 if isinstance(cube_slice, iris.cube.CubeList): 3408 ↛ 3412line 3408 didn't jump to line 3412 because the condition on line 3408 was always true
3409 single_cube = cube_slice[0]
3411 # Ensure valid stamp coordinate in cube dimensions
3412 if stamp_coordinate == "realization": 3412 ↛ 3415line 3412 didn't jump to line 3415 because the condition on line 3412 was always true
3413 stamp_coordinate = check_stamp_coordinate(single_cube)
3414 # Set plot titles and filename, based on sequence coordinate
3415 seq_coord = single_cube.coord(sequence_coordinate)
3416 # Use time coordinate in title and filename if single histogram output.
3417 if sequence_coordinate == "realization" and nplot == 1:
3418 seq_coord = single_cube.coord("time")
3419 # Use station name in title and filename if model vs obs comparison
3420 if sequence_coordinate == "station":
3421 seq_coord = single_cube.coord("Station_Name")
3423 plot_title, plot_filename = _set_title_and_filename(
3424 seq_coord, nplot, recipe_title, filename
3425 )
3427 # Do the actual plotting.
3428 plotting_func(
3429 cube_slice,
3430 filename=plot_filename,
3431 stamp_coordinate=stamp_coordinate,
3432 title=plot_title,
3433 vmin=vmin,
3434 vmax=vmax,
3435 hexbin=hexbin,
3436 )
3437 plot_index.append(plot_filename)
3439 # Add list of plots to plot metadata.
3440 complete_plot_index = _append_to_plot_index(plot_index)
3442 # Make a page to display the plots.
3443 _make_plot_html_page(complete_plot_index)
3445 return cubes
3448def _plot_and_save_postage_stamp_power_spectrum_series(
3449 cubes: iris.cube.Cube,
3450 coords: list[iris.coords.Coord],
3451 stamp_coordinate: str,
3452 filename: str,
3453 title: str,
3454 series_coordinate: str | None = None,
3455 **kwargs,
3456):
3457 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3459 Parameters
3460 ----------
3461 cubes: Cube or CubeList
3462 Cube or Cubelist of the power spectrum data.
3463 coords: list[Coord]
3464 Coordinates to plot on the x-axis, one per cube.
3465 stamp_coordinate: str
3466 Coordinate that becomes different plots.
3467 filename: str
3468 Filename of the plot to write.
3469 title: str
3470 Plot title.
3471 series_coordinate: str, optional
3472 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3474 """
3475 # Use the smallest square grid that will fit the members.
3476 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3478 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3479 model_colors_map = get_model_colors_map(cubes)
3480 # ax = plt.gca()
3481 # Make a subplot for each member.
3482 for member, subplot in zip(
3483 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3484 ):
3485 ax = plt.subplot(grid_size, grid_size, subplot)
3487 # Store min/max ranges.
3488 y_levels = []
3490 line_marker = None
3491 line_width = 1
3493 for cube in iter_maybe(member):
3494 xcoord = _select_series_coord(cube, series_coordinate)
3495 xname = xcoord.points
3497 yfield = cube.data # power spectrum
3498 label = None
3499 color = "black"
3500 if model_colors_map: 3500 ↛ 3501line 3500 didn't jump to line 3501 because the condition on line 3500 was never true
3501 label = cube.attributes.get("model_name")
3502 color = model_colors_map.get(label)
3504 if member.coord(stamp_coordinate).points == [0]:
3505 ax.plot(
3506 xname,
3507 yfield,
3508 color=color,
3509 marker=line_marker,
3510 ls="-",
3511 lw=line_width,
3512 label=f"{label} (control)"
3513 if len(cube.coord(stamp_coordinate).points) > 1
3514 else label,
3515 )
3516 # Label with member if part of an ensemble and not the control.
3517 else:
3518 ax.plot(
3519 xname,
3520 yfield,
3521 color=color,
3522 ls="-",
3523 lw=1.5,
3524 alpha=0.75,
3525 label=f"{label} (member)",
3526 )
3528 # Calculate the global min/max if multiple cubes are given.
3529 _, levels, _ = colorbar_map_levels(cube, axis="y")
3530 if levels is not None: 3530 ↛ 3531line 3530 didn't jump to line 3531 because the condition on line 3530 was never true
3531 y_levels.append(min(levels))
3532 y_levels.append(max(levels))
3534 # Add some labels and tweak the style.
3535 title = f"{title}"
3536 ax.set_title(title, fontsize=16)
3538 # Set appropriate x-axis label based on coordinate
3539 if series_coordinate == "wavelength" or ( 3539 ↛ 3542line 3539 didn't jump to line 3542 because the condition on line 3539 was never true
3540 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3541 ):
3542 ax.set_xlabel("Wavelength (km)", fontsize=14)
3543 elif series_coordinate == "physical_wavenumber" or ( 3543 ↛ 3548line 3543 didn't jump to line 3548 because the condition on line 3543 was always true
3544 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3545 ):
3546 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3547 else: # frequency or check units
3548 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3549 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3550 else:
3551 ax.set_xlabel("Wavenumber", fontsize=14)
3553 ax.set_ylabel("Power Spectral Density", fontsize=14)
3554 ax.tick_params(axis="both", labelsize=12)
3556 # Set log-log scale
3557 ax.set_xscale("log")
3558 ax.set_yscale("log")
3560 # Add gridlines
3561 ax.grid(linestyle="--", color="grey", linewidth=1)
3562 # Ientify unique labels for legend
3563 handles = list(
3564 {
3565 label: handle
3566 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3567 }.values()
3568 )
3569 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3571 ax = plt.gca()
3572 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3574 # Save plot.
3575 _save_close_figure(fig, "histogram postage stamp", filename)
3578def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3579 cubes: iris.cube.Cube,
3580 coords: list[iris.coords.Coord],
3581 stamp_coordinate: str,
3582 filename: str,
3583 title: str,
3584 series_coordinate: str | None = None,
3585 **kwargs,
3586):
3587 """Plot and save power spectra for ensemble members in single plot.
3589 Parameters
3590 ----------
3591 cubes: Cube or CubeList
3592 Cube or Cubelist of the power spectrum data.
3593 coords: list[Coord]
3594 Coordinates to plot on the x-axis, one per cube.
3595 stamp_coordinate: str
3596 Coordinate that becomes different plots.
3597 filename: str
3598 Filename of the plot to write.
3599 title: str
3600 Plot title.
3601 series_coordinate: str, optional
3602 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3604 """
3605 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3606 model_colors_map = get_model_colors_map(cubes)
3608 line_marker = None
3609 line_width = 1
3611 # Compute ensemble statistics to show spread
3612 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3613 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3614 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3616 xcoord_global = mean_cube.coord(series_coordinate)
3617 x_global = xcoord_global.points
3619 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3620 xcoord = _select_series_coord(member, series_coordinate)
3621 xname = xcoord.points
3623 yfield = member.data # power spectrum
3624 color = "black"
3625 if model_colors_map: 3625 ↛ 3629line 3625 didn't jump to line 3629 because the condition on line 3625 was always true
3626 label = member.attributes.get("model_name") if i == 0 else None
3627 color = model_colors_map.get(label)
3629 if member.coord(stamp_coordinate).points == [0]:
3630 ax.plot(
3631 xname,
3632 yfield,
3633 color=color,
3634 marker=line_marker,
3635 ls="-",
3636 lw=line_width,
3637 label=f"{label} (control)"
3638 if len(member.coord(stamp_coordinate).points) > 1
3639 else label,
3640 )
3641 # Label with member number if part of an ensemble and not the control.
3642 else:
3643 ax.plot(
3644 xname,
3645 yfield,
3646 color=color,
3647 ls="-",
3648 lw=1.5,
3649 alpha=0.75,
3650 label=label,
3651 )
3653 # Set appropriate x-axis label based on coordinate
3654 if series_coordinate == "wavelength" or ( 3654 ↛ 3657line 3654 didn't jump to line 3657 because the condition on line 3654 was never true
3655 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3656 ):
3657 ax.set_xlabel("Wavelength (km)", fontsize=14)
3658 elif series_coordinate == "physical_wavenumber" or ( 3658 ↛ 3663line 3658 didn't jump to line 3663 because the condition on line 3658 was always true
3659 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3660 ):
3661 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3662 else: # frequency or check units
3663 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3664 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3665 else:
3666 ax.set_xlabel("Wavenumber", fontsize=14)
3668 # Add ensemble spread shading
3669 ax.fill_between(
3670 x_global,
3671 min_cube.data,
3672 max_cube.data,
3673 color="grey",
3674 alpha=0.3,
3675 label="Ensemble spread",
3676 )
3678 # Add ensemble mean line
3679 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3681 ax.set_ylabel("Power Spectral Density", fontsize=14)
3682 ax.tick_params(axis="both", labelsize=12)
3684 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3685 # Set log-log scale
3686 ax.set_xscale("log")
3687 ax.set_yscale("log")
3689 # Add gridlines
3690 ax.grid(linestyle="--", color="grey", linewidth=1)
3691 # Identify unique labels for legend
3692 handles = list(
3693 {
3694 label: handle
3695 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3696 }.values()
3697 )
3698 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3700 # Figure title.
3701 ax.set_title(title, fontsize=16)
3703 # Save plot.
3704 _save_close_figure(fig, "power spectra postage stamp", filename)