Coverage for src/CSET/operators/plot.py: 83%
1174 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-20 14:18 +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: {np.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(cube.data):.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 strict_title: bool = False,
1881 overlay_cube: iris.cube.Cube | None = None,
1882 contour_cube: iris.cube.Cube | None = None,
1883 point_cube: iris.cube.Cube | None = None,
1884 **kwargs,
1885):
1886 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1888 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1889 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1890 is present then postage stamp plots will be produced.
1892 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1893 be overplotted on the same figure.
1895 Parameters
1896 ----------
1897 method: "contourf" | "pcolormesh" | "scatter"
1898 The plotting method to use.
1899 Select choice of "contourf" or "pcolormesh" for gridded data.
1900 Use "scatter" for point-based data.
1901 cube: Cube
1902 Iris cube of the data to plot. It should have two spatial dimensions,
1903 such as lat and lon, and may also have a another two dimension to be
1904 plotted sequentially and/or as postage stamp plots.
1905 filename: str | None
1906 Name of the plot to write, used as a prefix for plot sequences. If None
1907 uses the recipe name.
1908 sequence_coordinate: str
1909 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1910 This coordinate must exist in the cube.
1911 stamp_coordinate: str
1912 Coordinate about which to plot postage stamp plots. Defaults to
1913 ``"realization"``.
1914 strict_title: bool, optional
1915 Logical switch that if set to True will ensure that the MODEL_NAME
1916 string is not prepended to the plot title. The default is False.
1917 overlay_cube: Cube | None, optional
1918 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1919 contour_cube: Cube | None, optional
1920 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1921 point_cube: Cube | None, optional
1922 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
1924 Raises
1925 ------
1926 ValueError
1927 If the cube doesn't have the right dimensions.
1928 TypeError
1929 If the cube isn't a single cube.
1930 """
1931 # Ensure we've got a single cube.
1932 cube = check_single_cube(cube)
1934 # Set title based on recipe metadata or use cube name
1935 recipe_title = get_recipe_metadata().get("title", cube.name())
1937 # Check if there is a valid stamp coordinate in cube dimensions.
1938 if stamp_coordinate == "realization": 1938 ↛ 1943line 1938 didn't jump to line 1943 because the condition on line 1938 was always true
1939 stamp_coordinate = check_stamp_coordinate(cube)
1941 # Make postage stamp plots if stamp_coordinate exists and has more than a
1942 # single point.
1943 plotting_func = _plot_and_save_spatial_plot
1944 try:
1945 if cube.coord(stamp_coordinate).shape[0] > 1:
1946 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1947 except iris.exceptions.CoordinateNotFoundError:
1948 pass
1950 # Produce a geographical scatter plot if the data have a
1951 # dimension called observation or model_obs_error
1952 if any(
1953 crd.var_name == "station"
1954 or crd.var_name == "Station_Name"
1955 or crd.var_name == "model_obs_error"
1956 for crd in cube.coords()
1957 ):
1958 plotting_func = _plot_and_save_spatial_plot
1959 method = "scatter"
1961 # Must have a sequence coordinate.
1962 try:
1963 cube.coord(sequence_coordinate)
1964 except iris.exceptions.CoordinateNotFoundError as err:
1965 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1967 # Create a plot for each value of the sequence coordinate.
1968 plot_index = []
1969 nplot = np.size(cube.coord(sequence_coordinate).points)
1971 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1972 # Set plot titles and filename
1973 seq_coord = cube_slice.coord(sequence_coordinate)
1975 if "model_name" in cube.attributes: 1975 ↛ 1976line 1975 didn't jump to line 1976 because the condition on line 1975 was never true
1976 model_name = cube.attributes["model_name"]
1977 else:
1978 model_name = None
1979 if strict_title: 1979 ↛ 1980line 1979 didn't jump to line 1980 because the condition on line 1979 was never true
1980 model_name = None
1982 plot_title, plot_filename = _set_title_and_filename(
1983 seq_coord, nplot, recipe_title, filename, model_name=model_name
1984 )
1986 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1987 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1988 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1989 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1991 # Do the actual plotting.
1992 plotting_func(
1993 cube_slice,
1994 filename=plot_filename,
1995 stamp_coordinate=stamp_coordinate,
1996 title=plot_title,
1997 method=method,
1998 overlay_cube=overlay_slice,
1999 contour_cube=contour_slice,
2000 point_cube=point_slice,
2001 **kwargs,
2002 )
2003 plot_index.append(plot_filename)
2005 # Add list of plots to plot metadata.
2006 complete_plot_index = _append_to_plot_index(plot_index)
2008 # Make a page to display the plots.
2009 _make_plot_html_page(complete_plot_index)
2012####################
2013# Public functions #
2014####################
2017def spatial_contour_plot(
2018 cube: iris.cube.Cube,
2019 filename: str | None = None,
2020 sequence_coordinate: str = "time",
2021 stamp_coordinate: str = "realization",
2022 **kwargs,
2023) -> iris.cube.Cube:
2024 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2026 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2027 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2028 is present then postage stamp plots will be produced.
2030 Parameters
2031 ----------
2032 cube: Cube
2033 Iris cube of the data to plot. It should have two spatial dimensions,
2034 such as lat and lon, and may also have a another two dimension to be
2035 plotted sequentially and/or as postage stamp plots.
2036 filename: str, optional
2037 Name of the plot to write, used as a prefix for plot sequences. Defaults
2038 to the recipe name.
2039 sequence_coordinate: str, optional
2040 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2041 This coordinate must exist in the cube.
2042 stamp_coordinate: str, optional
2043 Coordinate about which to plot postage stamp plots. Defaults to
2044 ``"realization"``.
2046 Returns
2047 -------
2048 Cube
2049 The original cube (so further operations can be applied).
2051 Raises
2052 ------
2053 ValueError
2054 If the cube doesn't have the right dimensions.
2055 TypeError
2056 If the cube isn't a single cube.
2057 """
2058 _spatial_plot(
2059 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2060 )
2061 return cube
2064def spatial_pcolormesh_plot(
2065 cubes: iris.cube.Cube | iris.cube.CubeList,
2066 filename: str | None = None,
2067 sequence_coordinate: str = "time",
2068 stamp_coordinate: str = "realization",
2069 strict_title: bool = False,
2070 **kwargs,
2071) -> iris.cube.Cube:
2072 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
2074 A 2D spatial field can be plotted, but if the sequence_coordinate is present
2075 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2076 is present then postage stamp plots will be produced.
2078 This function is significantly faster than ``spatial_contour_plot``,
2079 especially at high resolutions, and should be preferred unless contiguous
2080 contour areas are important.
2082 Parameters
2083 ----------
2084 cube: Cubes
2085 Iris cube or cubelist of the data to plot. Each cube should have two spatial dimensions,
2086 such as lat and lon, and may also have a another two dimension to be
2087 plotted sequentially and/or as postage stamp plots.
2088 filename: str, optional
2089 Name of the plot to write, used as a prefix for plot sequences. Defaults
2090 to the recipe name.
2091 sequence_coordinate: str, optional
2092 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2093 This coordinate must exist in the cube.
2094 stamp_coordinate: str, optional
2095 Coordinate about which to plot postage stamp plots. Defaults to
2096 ``"realization"``.
2097 strict_title: bool, optional
2098 Logical switch that if set to True will ensure that the MODEL_NAME
2099 string is not prepended to the plot title. The default is False.
2101 Returns
2102 -------
2103 Cubes
2104 The original cube/cubelist (so further operations can be applied).
2106 Raises
2107 ------
2108 ValueError
2109 If the cube doesn't have the right dimensions.
2110 """
2111 if isinstance(cubes, iris.cube.CubeList): 2111 ↛ 2112line 2111 didn't jump to line 2112 because the condition on line 2111 was never true
2112 for model_cube in cubes:
2113 _spatial_plot(
2114 "pcolormesh",
2115 model_cube,
2116 filename,
2117 sequence_coordinate,
2118 stamp_coordinate,
2119 strict_title,
2120 **kwargs,
2121 )
2122 elif isinstance(cubes, iris.cube.Cube): 2122 ↛ 2132line 2122 didn't jump to line 2132 because the condition on line 2122 was always true
2123 _spatial_plot(
2124 "pcolormesh",
2125 cubes,
2126 filename,
2127 sequence_coordinate,
2128 stamp_coordinate,
2129 strict_title,
2130 **kwargs,
2131 )
2132 return cubes
2135def spatial_multi_pcolormesh_plot(
2136 cube: iris.cube.Cube,
2137 overlay_cube: iris.cube.Cube | None = None,
2138 contour_cube: iris.cube.Cube | None = None,
2139 point_cube: iris.cube.Cube | None = None,
2140 filename: str | None = None,
2141 sequence_coordinate: str = "time",
2142 stamp_coordinate: str = "realization",
2143 **kwargs,
2144) -> iris.cube.Cube:
2145 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2147 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2148 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2149 is present then postage stamp plots will be produced.
2151 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2153 If specified, contours of a contour_cube can be overplotted on top of those.
2155 If specified, a spatial scatter map of point_cube can be overplotted.
2157 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2159 This function is significantly faster than ``spatial_contour_plot``,
2160 especially at high resolutions, and should be preferred unless contiguous
2161 contour areas are important.
2163 Parameters
2164 ----------
2165 cube: Cube
2166 Iris cube of the data to plot. It should have two spatial dimensions,
2167 such as lat and lon, and may also have two additional dimensions to be
2168 plotted sequentially and/or as postage stamp plots.
2169 overlay_cube: Cube, optional
2170 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2171 such as lat and lon, and may also have two additional dimensions to be
2172 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.
2173 If not provided, output plot generated without overlay cube.
2174 contour_cube: Cube, optional
2175 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,
2176 such as lat and lon, and may also have two additional dimensions to be
2177 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2178 point_cube: Cube, optional
2179 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
2180 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
2181 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
2182 postage stamp plots. If not provided, output plot generated without point-based layer.
2183 filename: str, optional
2184 Name of the plot to write, used as a prefix for plot sequences. Defaults
2185 to the recipe name.
2186 sequence_coordinate: str, optional
2187 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2188 This coordinate must exist in the cube.
2189 stamp_coordinate: str, optional
2190 Coordinate about which to plot postage stamp plots. Defaults to
2191 ``"realization"``.
2193 Returns
2194 -------
2195 Cube
2196 The original cube (so further operations can be applied).
2198 Raises
2199 ------
2200 ValueError
2201 If the cube doesn't have the right dimensions.
2202 TypeError
2203 If the cube isn't a single cube.
2204 """
2205 _spatial_plot(
2206 "pcolormesh",
2207 cube,
2208 filename,
2209 sequence_coordinate,
2210 stamp_coordinate,
2211 overlay_cube=overlay_cube,
2212 contour_cube=contour_cube,
2213 point_cube=point_cube,
2214 )
2215 return cube, overlay_cube, contour_cube, point_cube
2218# TODO: Expand function to handle ensemble data.
2219# line_coordinate: str, optional
2220# Coordinate about which to plot multiple lines. Defaults to
2221# ``"realization"``.
2222def plot_line_series(
2223 cube: iris.cube.Cube | iris.cube.CubeList,
2224 filename: str | None = None,
2225 series_coordinate: str = "time",
2226 sequence_coordinate: str = "time",
2227 # add the following for ensembles
2228 stamp_coordinate: str = "realization",
2229 single_plot: bool = False,
2230 **kwargs,
2231) -> iris.cube.Cube | iris.cube.CubeList:
2232 """Plot a line plot for the specified coordinate.
2234 The Cube or CubeList must be 1D.
2236 Parameters
2237 ----------
2238 iris.cube | iris.cube.CubeList
2239 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2240 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2241 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2242 filename: str, optional
2243 Name of the plot to write, used as a prefix for plot sequences. Defaults
2244 to the recipe name.
2245 series_coordinate: str, optional
2246 Coordinate about which to make a series. Defaults to ``"time"``. This
2247 coordinate must exist in the cube.
2249 Returns
2250 -------
2251 iris.cube.Cube | iris.cube.CubeList
2252 The original Cube or CubeList (so further operations can be applied).
2254 Raises
2255 ------
2256 ValueError
2257 If the cubes don't have the right dimensions.
2258 TypeError
2259 If the cube isn't a Cube or CubeList.
2260 """
2261 # Ensure we have a name for the plot file.
2262 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2264 num_models = get_num_models(cube)
2266 validate_cube_shape(cube, num_models)
2268 # Iterate over all cubes and extract coordinate to plot.
2269 cubes = iris.cube.CubeList(iter_maybe(cube))
2271 coords = []
2272 for model_cube in cubes:
2273 try:
2274 coords.append(model_cube.coord(series_coordinate))
2275 except iris.exceptions.CoordinateNotFoundError as err:
2276 raise ValueError(
2277 f"Cube must have a {series_coordinate} coordinate."
2278 ) from err
2279 # Count cube dimensions and exclude realization and
2280 # forecast_reference_time if they exist.
2281 ndim = model_cube.ndim
2283 if model_cube.coords("realization"): 2283 ↛ 2291line 2283 didn't jump to line 2291 because the condition on line 2283 was always true
2284 # returns coord dimension
2285 realization_dims = model_cube.coord_dims("realization")
2287 # Only subtract if realization is a dimension coordinate
2288 if realization_dims:
2289 ndim -= len(realization_dims)
2291 if model_cube.coords("forecast_reference_time"):
2292 frt_dims = model_cube.coord_dims("forecast_reference_time")
2294 # Only subtract if frt is a dimension coordinate
2295 if frt_dims: 2295 ↛ 2296line 2295 didn't jump to line 2296 because the condition on line 2295 was never true
2296 ndim -= len(frt_dims)
2298 if ndim > 2:
2299 raise ValueError(
2300 "Cube must be 1D or 2D (excluding any realization or forecast_reference_time dimensions)."
2301 )
2303 plot_index = []
2305 # Check if this is a spectral plot by looking for spectral coordinates
2306 is_spectral_plot = series_coordinate in [
2307 "frequency",
2308 "physical_wavenumber",
2309 "wavelength",
2310 ]
2312 if is_spectral_plot:
2313 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2314 # coordinate frequency/wavenumber.
2315 # If several power spectra are plotted with time as sequence_coordinate for the
2316 # time slider option.
2318 # Internal plotting function.
2319 plotting_func = _plot_and_save_line_power_spectrum_series
2321 for model_cube in cubes:
2322 try:
2323 model_cube.coord(sequence_coordinate)
2324 except iris.exceptions.CoordinateNotFoundError as err:
2325 raise ValueError(
2326 f"Cube must have a {sequence_coordinate} coordinate."
2327 ) from err
2329 if num_models == 1: 2329 ↛ 2344line 2329 didn't jump to line 2344 because the condition on line 2329 was always true
2330 # check for ensembles
2331 if ( 2331 ↛ 2335line 2331 didn't jump to line 2335 because the condition on line 2331 was never true
2332 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2333 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2334 ):
2335 if single_plot:
2336 # Plot spectra, mean and ensemble spread on 1 plot
2337 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2338 else:
2339 # Plot postage stamps
2340 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2341 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2342 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2343 else:
2344 all_points = sorted(
2345 set(
2346 itertools.chain.from_iterable(
2347 cb.coord(sequence_coordinate).points for cb in cubes
2348 )
2349 )
2350 )
2351 all_slices = list(
2352 itertools.chain.from_iterable(
2353 cb.slices_over(sequence_coordinate) for cb in cubes
2354 )
2355 )
2356 # Matched slices (matched by seq coord point; it may happen that
2357 # evaluated models do not cover the same seq coord range, hence matching
2358 # necessary)
2359 cube_iterables = [
2360 iris.cube.CubeList(
2361 s
2362 for s in all_slices
2363 if s.coord(sequence_coordinate).points[0] == point
2364 )
2365 for point in all_points
2366 ]
2367 nplot = len(all_points)
2369 # Create a plot for each value of the sequence coordinate. Allowing for
2370 # multiple cubes in a CubeList to be plotted in the same plot for similar
2371 # sequence values. Passing a CubeList into the internal plotting function
2372 # for similar values of the sequence coordinate. cube_slice can be an
2373 # iris.cube.Cube or an iris.cube.CubeList.
2375 for cube_slice in cube_iterables:
2376 # Normalize cube_slice to a list of cubes
2377 if isinstance(cube_slice, iris.cube.CubeList): 2377 ↛ 2378line 2377 didn't jump to line 2378 because the condition on line 2377 was never true
2378 cubes = list(cube_slice)
2379 elif isinstance(cube_slice, iris.cube.Cube): 2379 ↛ 2382line 2379 didn't jump to line 2382 because the condition on line 2379 was always true
2380 cubes = [cube_slice]
2381 else:
2382 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2384 # Use sequence value so multiple sequences can merge.
2385 seq_coord = cube_slice[0].coord(sequence_coordinate)
2386 plot_title, plot_filename = _set_title_and_filename(
2387 seq_coord, nplot, recipe_title, filename
2388 )
2390 # Format the coordinate value in a unit appropriate way.
2391 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2393 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2394 if nplot == 1 and seq_coord.has_bounds and np.size(seq_coord.bounds) > 1: 2394 ↛ 2395line 2394 didn't jump to line 2395 because the condition on line 2394 was never true
2395 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2397 # Do the actual plotting.
2398 plotting_func(
2399 cube_slice,
2400 coords,
2401 stamp_coordinate,
2402 plot_filename,
2403 title,
2404 series_coordinate,
2405 )
2407 plot_index.append(plot_filename)
2408 else:
2409 # Format the title and filename using plotted series coordinate
2410 nplot = 1
2411 seq_coord = coords[0]
2412 plot_title, plot_filename = _set_title_and_filename(
2413 seq_coord, nplot, recipe_title, filename
2414 )
2416 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2417 if (
2418 "station" in [c.name() for c in cubes[0].coords()]
2419 and len(cubes[0].coord("station").points) > 1
2420 ):
2421 for station in cubes[0].coord("station").points:
2422 station_cubes = cubes.extract(iris.Constraint(station=station))
2423 station_name = station_cubes[0].coord("Station_Name").points[0]
2424 station_plotname = plot_filename.replace(
2425 ".png", "_" + station_name + ".png"
2426 )
2427 _plot_and_save_line_series(
2428 station_cubes,
2429 coords,
2430 "realization",
2431 station_plotname,
2432 f"{plot_title} {station_name}",
2433 )
2434 plot_index.append(station_plotname)
2436 else:
2437 # Do the actual plotting for all other series coordinate options.
2438 _plot_and_save_line_series(
2439 cubes, coords, stamp_coordinate, plot_filename, plot_title
2440 )
2442 plot_index.append(plot_filename)
2444 # append plot to list of plots
2445 complete_plot_index = _append_to_plot_index(plot_index)
2447 # Make a page to display the plots.
2448 _make_plot_html_page(complete_plot_index)
2450 return cube
2453def plot_vertical_line_series(
2454 cubes: iris.cube.Cube | iris.cube.CubeList,
2455 filename: str | None = None,
2456 series_coordinate: str = "model_level_number",
2457 sequence_coordinate: str = "time",
2458 # line_coordinate: str = "realization",
2459 **kwargs,
2460) -> iris.cube.Cube | iris.cube.CubeList:
2461 """Plot a line plot against a type of vertical coordinate.
2463 The Cube or CubeList must be 1D.
2465 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2466 then a sequence of plots will be produced.
2468 Parameters
2469 ----------
2470 iris.cube | iris.cube.CubeList
2471 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2472 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2473 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2474 filename: str, optional
2475 Name of the plot to write, used as a prefix for plot sequences. Defaults
2476 to the recipe name.
2477 series_coordinate: str, optional
2478 Coordinate to plot on the y-axis. Can be ``pressure`` or
2479 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2480 for LFRic. Defaults to ``model_level_number``.
2481 This coordinate must exist in the cube.
2482 sequence_coordinate: str, optional
2483 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2484 This coordinate must exist in the cube.
2486 Returns
2487 -------
2488 iris.cube.Cube | iris.cube.CubeList
2489 The original Cube or CubeList (so further operations can be applied).
2490 Plotted data.
2492 Raises
2493 ------
2494 ValueError
2495 If the cubes doesn't have the right dimensions.
2496 TypeError
2497 If the cube isn't a Cube or CubeList.
2498 """
2499 # Ensure we have a name for the plot file.
2500 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2502 cubes = iter_maybe(cubes)
2503 # Initialise empty list to hold all data from all cubes in a CubeList
2504 all_data = []
2506 # Store min/max ranges for x range.
2507 x_levels = []
2509 num_models = get_num_models(cubes)
2511 validate_cube_shape(cubes, num_models)
2513 # Iterate over all cubes in cube or CubeList and plot.
2514 coords = []
2515 for cube in cubes:
2516 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2517 try:
2518 coords.append(cube.coord(series_coordinate))
2519 except iris.exceptions.CoordinateNotFoundError as err:
2520 raise ValueError(
2521 f"Cube must have a {series_coordinate} coordinate."
2522 ) from err
2524 try:
2525 if cube.ndim > 1 or not cube.coords("realization"): 2525 ↛ 2533line 2525 didn't jump to line 2533 because the condition on line 2525 was always true
2526 cube.coord(sequence_coordinate)
2527 except iris.exceptions.CoordinateNotFoundError as err:
2528 raise ValueError(
2529 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2530 ) from err
2532 # Get minimum and maximum from levels information.
2533 _, levels, _ = colorbar_map_levels(cube, axis="x")
2534 if levels is not None: 2534 ↛ 2538line 2534 didn't jump to line 2538 because the condition on line 2534 was always true
2535 x_levels.append(min(levels))
2536 x_levels.append(max(levels))
2537 else:
2538 all_data.append(cube.data)
2540 if len(x_levels) == 0: 2540 ↛ 2542line 2540 didn't jump to line 2542 because the condition on line 2540 was never true
2541 # Combine all data into a single NumPy array
2542 combined_data = np.concatenate(all_data)
2544 # Set the lower and upper limit for the x-axis to ensure all plots have
2545 # same range. This needs to read the whole cube over the range of the
2546 # sequence and if applicable postage stamp coordinate.
2547 vmin = np.floor(combined_data.min())
2548 vmax = np.ceil(combined_data.max())
2549 else:
2550 vmin = min(x_levels)
2551 vmax = max(x_levels)
2553 # Check if the cube has a sequence coordinate (e.g. time). If not, plot
2554 # a single profile directly without iterating over a sequence.
2555 sequence_coords = [
2556 cube.coord(sequence_coordinate)
2557 for cube in cubes
2558 if cube.coords(sequence_coordinate)
2559 ]
2560 has_sequence_coord = len(sequence_coords) == len(cubes) and all(
2561 np.size(coord.points) > 1 for coord in sequence_coords
2562 )
2563 has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all(
2564 np.size(coord.points) == 1 for coord in sequence_coords
2565 )
2567 plot_index = []
2568 if has_sequence_coord: 2568 ↛ 2593line 2568 didn't jump to line 2593 because the condition on line 2568 was always true
2569 # Matching the slices (matching by seq coord point; it may happen that
2570 # evaluated models do not cover the same seq coord range, hence matching
2571 # necessary)
2572 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2573 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2574 for cubes_slice in cube_iterables:
2575 # Format the coordinate value in a unit appropriate way.
2576 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2577 plot_title, plot_filename = _set_title_and_filename(
2578 seq_coord, nplot, recipe_title, filename
2579 )
2581 # Do the actual plotting.
2582 _plot_and_save_vertical_line_series(
2583 cubes_slice,
2584 coords,
2585 "realization",
2586 plot_filename,
2587 series_coordinate,
2588 title=plot_title,
2589 vmin=vmin,
2590 vmax=vmax,
2591 )
2592 plot_index.append(plot_filename)
2593 elif has_scalar_sequence_coord:
2594 # Scalar sequence coordinate (typically aggregated time bounds):
2595 # make one plot and include sequence period in title/filename.
2596 plot_title, plot_filename = _set_title_and_filename(
2597 sequence_coords[0], 1, recipe_title, filename
2598 )
2600 _plot_and_save_vertical_line_series(
2601 cubes,
2602 coords,
2603 "realization",
2604 plot_filename,
2605 series_coordinate,
2606 title=plot_title,
2607 vmin=vmin,
2608 vmax=vmax,
2609 )
2610 plot_index.append(plot_filename)
2611 else:
2612 # 1D case: no sequence coordinate, plot a single profile.
2613 plot_title = recipe_title
2614 if filename:
2615 plot_filename = filename
2616 else:
2617 plot_filename = f"{slugify(plot_title)}.png"
2619 _plot_and_save_vertical_line_series(
2620 cubes,
2621 coords,
2622 "realization",
2623 plot_filename,
2624 series_coordinate,
2625 title=plot_title,
2626 vmin=vmin,
2627 vmax=vmax,
2628 )
2629 plot_index.append(plot_filename)
2631 # Add list of plots to plot metadata.
2632 complete_plot_index = _append_to_plot_index(plot_index)
2634 # Make a page to display the plots.
2635 _make_plot_html_page(complete_plot_index)
2637 return cubes
2640def qq_plot(
2641 cubes: iris.cube.CubeList,
2642 coordinates: list[str],
2643 percentiles: list[float],
2644 model_names: list[str],
2645 filename: str | None = None,
2646 one_to_one: bool = True,
2647 **kwargs,
2648) -> iris.cube.CubeList:
2649 """Plot a Quantile-Quantile plot between two models for common time points.
2651 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2652 collapsed within the operator over all specified coordinates such as
2653 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2654 ensemble members to ensure a 1D cube (array).
2656 Parameters
2657 ----------
2658 cubes: iris.cube.CubeList
2659 Two cubes of the same variable with different models.
2660 coordinate: list[str]
2661 The list of coordinates to collapse over. This list should be
2662 every coordinate within the cube to result in a 1D cube around
2663 the percentile coordinate.
2664 percent: list[float]
2665 A list of percentiles to appear in the plot.
2666 model_names: list[str]
2667 A list of model names to appear on the axis of the plot.
2668 filename: str, optional
2669 Filename of the plot to write.
2670 one_to_one: bool, optional
2671 If True a 1:1 line is plotted; if False it is not. Default is True.
2673 Raises
2674 ------
2675 ValueError
2676 When the cubes are not compatible.
2678 Notes
2679 -----
2680 The quantile-quantile plot is a variant on the scatter plot representing
2681 two datasets by their quantiles (percentiles) for common time points.
2682 This plot does not use a theoretical distribution to compare against, but
2683 compares percentiles of two datasets. This plot does
2684 not use all raw data points, but plots the selected percentiles (quantiles) of
2685 each variable instead for the two datasets, thereby normalising the data for a
2686 direct comparison between the selected percentiles of the two dataset distributions.
2688 Quantile-quantile plots are valuable for comparing against
2689 observations and other models. Identical percentiles between the variables
2690 will lie on the one-to-one line implying the values correspond well to each
2691 other. Where there is a deviation from the one-to-one line a range of
2692 possibilities exist depending on how and where the data is shifted (e.g.,
2693 Wilks 2011 [Wilks2011]_).
2695 For distributions above the one-to-one line the distribution is left-skewed;
2696 below is right-skewed. A distinct break implies a bimodal distribution, and
2697 closer values/values further apart at the tails imply poor representation of
2698 the extremes.
2700 """
2701 # Check cubes using same functionality as the difference operator.
2702 if len(cubes) != 2:
2703 raise ValueError("cubes should contain exactly 2 cubes.")
2704 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2705 other: Cube = cubes.extract_cube(
2706 iris.Constraint(
2707 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2708 )
2709 )
2711 # Get spatial coord names.
2712 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2713 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2715 # Ensure cubes to compare are on common differencing grid.
2716 # This is triggered if either
2717 # i) latitude and longitude shapes are not the same. Note grid points
2718 # are not compared directly as these can differ through rounding
2719 # errors.
2720 # ii) or variables are known to often sit on different grid staggering
2721 # in different models (e.g. cell center vs cell edge), as is the case
2722 # for UM and LFRic comparisons.
2723 # In future greater choice of regridding method might be applied depending
2724 # on variable type. Linear regridding can in general be appropriate for smooth
2725 # variables. Care should be taken with interpretation of differences
2726 # given this dependency on regridding.
2727 if (
2728 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2729 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2730 ) or (
2731 base.long_name
2732 in [
2733 "eastward_wind_at_10m",
2734 "northward_wind_at_10m",
2735 "northward_wind_at_cell_centres",
2736 "eastward_wind_at_cell_centres",
2737 "zonal_wind_at_pressure_levels",
2738 "meridional_wind_at_pressure_levels",
2739 "potential_vorticity_at_pressure_levels",
2740 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2741 ]
2742 ):
2743 logger.debug("Linear regridding base cube to other grid to compute differences")
2744 base = regrid_onto_cube(base, other, method="Linear")
2746 # Extract just common time points.
2747 base, other = _extract_common_time_points(base, other)
2749 # Equalise attributes so we can merge.
2750 fully_equalise_attributes([base, other])
2751 logger.debug("Base: %s\nOther: %s", base, other)
2753 # Collapse cubes.
2754 base = collapse(
2755 base,
2756 coordinate=coordinates,
2757 method="PERCENTILE",
2758 additional_percent=percentiles,
2759 )
2760 other = collapse(
2761 other,
2762 coordinate=coordinates,
2763 method="PERCENTILE",
2764 additional_percent=percentiles,
2765 )
2767 # Ensure we have a name for the plot file.
2768 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2769 title = f"{recipe_title}"
2771 if filename is None:
2772 filename = slugify(recipe_title)
2774 # Add file extension.
2775 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2777 # Do the actual plotting on a scatter plot
2778 _plot_and_save_scatter_plot(
2779 base, other, plot_filename, title, one_to_one, model_names
2780 )
2782 # Add list of plots to plot metadata.
2783 plot_index = _append_to_plot_index([plot_filename])
2785 # Make a page to display the plots.
2786 _make_plot_html_page(plot_index)
2788 return iris.cube.CubeList([base, other])
2791def hinton_plot(
2792 cubes: iris.cube.CubeList, base_name: str, other_name: str, magnitude: bool = False
2793) -> None:
2794 """
2795 Plot a Hinton style triangle/scorecard plot.
2797 This plot type can be useful for summarising high level information, such as comparing
2798 how 'skillful' two models are when verified against observations for a variety of metrics,
2799 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2800 than customisable by the user as input arguments; many have been designed to automatically
2801 scale the plot depending on the number of x and y components.
2803 Parameters
2804 ----------
2805 cubes: iris.cube.CubeList
2806 A iris cubelist, containing at least two cubes of a skill metric to plot (model vs obs). Can
2807 include multiple variables, in which the plot will automatically scale for. If cubes containing
2808 a name significance_<var> exist, containing a bool array, then it will also plot whether each
2809 triangle is significant by using a thick black outline. Each cube should be 1D, with
2810 forecast_period as the only dimension.
2811 base_name: str
2812 The name of the base model to use as the control in the Hinton plot, as a string.
2813 other_name: str
2814 The name of the other model to use as the test in the Hinton plot, as a string.
2815 magnitude: bool
2816 Option bool when if True, then plot the numerical value difference in two values under each
2817 triangle.
2818 """
2819 # Ensure we have a name for the plot file.
2820 recipe_title = get_recipe_metadata().get("title", "Hinton")
2821 title = f"{recipe_title}"
2822 filename = slugify(recipe_title)
2824 # Check that all cubes only have one dimension called forecast_period
2825 for cube in cubes:
2826 if len(cube.dim_coords) > 1:
2827 raise ValueError(f"Should only have one dimension coord, {cube}")
2828 if cube.dim_coords[0].name() != "forecast_period":
2829 raise ValueError(
2830 f"Single coord should be forecast_period, not {cube.dim_coords[0].name()}"
2831 )
2833 # Separate out base cubes and other cubes.
2834 base_cubes = iris.cube.CubeList()
2835 other_cubes = iris.cube.CubeList()
2836 for c in cubes:
2837 if c.attributes["model_name"] == base_name:
2838 base_cubes.append(c)
2839 elif c.attributes["model_name"] == other_name:
2840 other_cubes.append(c)
2842 # base cubes should be the same length as other cubes, otherwise one is missing a variable.
2843 if len(base_cubes) != len(other_cubes):
2844 raise ValueError(
2845 f"base cubes {base_cubes} are not same number as {other_cubes}"
2846 )
2848 # Find common variable names in the two groups.
2849 base_vars = {cube.long_name for cube in base_cubes if cube.long_name is not None}
2850 other_vars = {cube.long_name for cube in other_cubes if cube.long_name is not None}
2851 common_vars = sorted(base_vars & other_vars)
2853 # Iterate over each variable (row)
2854 rows = []
2855 for var in common_vars:
2856 # Extract cube with matching variable name
2857 base_cube = next(
2858 (c for c in base_cubes if c.long_name == var),
2859 None,
2860 )
2862 other_cube = next(
2863 (c for c in other_cubes if c.long_name == var),
2864 None,
2865 )
2867 # If we can't find a variable in both cubes, then skip
2868 if base_cube is None or other_cube is None: 2868 ↛ 2869line 2868 didn't jump to line 2869 because the condition on line 2868 was never true
2869 continue
2871 # Compute difference (1D array)
2872 # We can already make assumption both on same forecast_periods as checked
2873 # prior to computing metric.
2874 diff = other_cube.data - base_cube.data
2876 # See if there is a significance cube present, if not, set as None.
2877 sig_cube = next(
2878 (cube for cube in cubes if cube.long_name == f"significance_{var}"),
2879 None,
2880 )
2882 # Append row information.
2883 rows.append(
2884 {
2885 "name": var,
2886 "forecast_periods": base_cube.coord("forecast_period").points,
2887 "change": diff,
2888 "significance": sig_cube.data.astype(bool)
2889 if sig_cube is not None
2890 else None,
2891 }
2892 )
2894 # For each row, compute standardised anomalies
2895 for row in rows:
2896 change = np.asarray(row["change"])
2898 anoms = change - np.mean(change)
2900 scale = np.max(np.abs(anoms))
2902 if scale > 0:
2903 scaled = anoms / scale
2904 else:
2905 scaled = np.zeros_like(anoms)
2907 row["anoms"] = anoms
2908 row["scaled"] = scaled
2910 # Setup colors of triangles
2911 color_pos = "#7CAE00"
2912 color_neg = "#7B68EE"
2914 # Setup cell/text size ratios
2915 figsize = None
2916 cell_size_in = 1.5
2917 text_row_ratio = 0.25
2919 # Get the number of x and y elements
2920 ny = len(rows)
2921 nx = max(len(row["forecast_periods"]) for row in rows)
2923 # Build non-uniform y coordinates
2924 tri_height = 1.0
2925 txt_height = text_row_ratio
2927 tri_y = []
2928 txt_y = []
2929 y_edges = [0.0]
2931 y = 0.0
2932 for _j in range(ny):
2933 tri_y.append(y + tri_height / 2)
2934 y += tri_height
2935 y_edges.append(y)
2937 if magnitude:
2938 txt_y.append(y + txt_height / 2)
2939 y += txt_height
2940 y_edges.append(y)
2942 total_height = y
2944 # Dynamic figure size
2945 if figsize is None: 2945 ↛ 2950line 2945 didn't jump to line 2950 because the condition on line 2945 was always true
2946 width = nx * cell_size_in
2947 height = total_height * cell_size_in + 2
2948 figsize = (width, height)
2950 fig, ax = plt.subplots(figsize=figsize)
2952 # Setup axes and grid.
2953 ax.set_aspect("equal", adjustable="box")
2954 ax.set_xlim(-0.5, nx - 0.5)
2955 ax.set_ylim(0, total_height)
2957 longest_row = max(rows, key=lambda row: len(row["forecast_periods"]))
2959 ax.set_xticks(np.arange(nx))
2960 ax.set_xticklabels(
2961 longest_row["forecast_periods"],
2962 rotation=90,
2963 )
2965 ax.set_yticks(tri_y)
2966 ax.set_yticklabels([row["name"] for row in rows])
2968 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2969 ax.set_yticks(y_edges, minor=True)
2971 ax.set_axisbelow(True)
2972 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2973 ax.grid(False, which="major")
2974 ax.tick_params(which="minor", length=0)
2976 ax.invert_yaxis()
2978 # Compute marker scaling (fixed overlap)
2979 fig.canvas.draw()
2981 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2982 width_in, height_in = bbox.width, bbox.height
2984 cell_w = (width_in * fig.dpi) / nx
2985 cell_h = (height_in * fig.dpi) / total_height
2986 cell_pixels = min(cell_w, cell_h)
2988 max_marker_size = (0.6 * cell_pixels) ** 2
2990 text_fontsize = cell_pixels * 0.15
2992 # Plot triangles + text
2993 for j, row in enumerate(rows):
2994 scaled = row["scaled"]
2995 anoms = row["anoms"]
2996 signif = row["significance"]
2998 for i in range(len(scaled)):
2999 val = scaled[i]
3001 if np.isnan(val): 3001 ↛ 3002line 3001 didn't jump to line 3002 because the condition on line 3001 was never true
3002 continue
3004 if abs(val) < 0.01:
3005 continue
3007 size = max_marker_size * abs(val)
3009 if val >= 0:
3010 marker = "^"
3011 color = color_pos
3012 else:
3013 marker = "v"
3014 color = color_neg
3016 if signif is not None:
3017 sig = bool(signif[i])
3018 edgecolor = "black" if sig else "none"
3019 linewidth = 0.6 if sig else 0.0
3020 else:
3021 edgecolor = "none"
3022 linewidth = 0.0
3024 ax.scatter(
3025 i,
3026 tri_y[j],
3027 s=size,
3028 marker=marker,
3029 c=color,
3030 edgecolors=edgecolor,
3031 linewidths=linewidth,
3032 zorder=3,
3033 clip_on=True,
3034 )
3036 # Text row
3037 if magnitude:
3038 mag_val = anoms[i]
3040 if not np.isnan(mag_val): 3040 ↛ 2998line 3040 didn't jump to line 2998 because the condition on line 3040 was always true
3041 ax.text(
3042 i,
3043 txt_y[j],
3044 f"{mag_val:.1f}",
3045 ha="center",
3046 va="center",
3047 fontsize=text_fontsize,
3048 color="black",
3049 zorder=4,
3050 )
3052 ax.set_title(title)
3053 plt.tight_layout()
3055 # Save plot.
3056 _save_close_figure(fig, "hinton", filename)
3058 # Add file extension.
3059 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3061 # Add list of plots to plot metadata.
3062 plot_index = _append_to_plot_index([plot_filename])
3064 # Make a page to display the plots.
3065 _make_plot_html_page(plot_index)
3068def scatter_plot(
3069 cube_x: iris.cube.Cube | iris.cube.CubeList,
3070 cube_y: iris.cube.Cube | iris.cube.CubeList,
3071 filename: str | None = None,
3072 one_to_one: bool = True,
3073 **kwargs,
3074) -> iris.cube.CubeList:
3075 """Plot a scatter plot between two variables.
3077 Both cubes must be 1D.
3079 Parameters
3080 ----------
3081 cube_x: Cube | CubeList
3082 1 dimensional Cube of the data to plot on y-axis.
3083 cube_y: Cube | CubeList
3084 1 dimensional Cube of the data to plot on x-axis.
3085 filename: str, optional
3086 Filename of the plot to write.
3087 one_to_one: bool, optional
3088 If True a 1:1 line is plotted; if False it is not. Default is True.
3090 Returns
3091 -------
3092 cubes: CubeList
3093 CubeList of the original x and y cubes for further processing.
3095 Raises
3096 ------
3097 ValueError
3098 If the cube doesn't have the right dimensions and cubes not the same
3099 size.
3100 TypeError
3101 If the cube isn't a single cube.
3103 Notes
3104 -----
3105 Scatter plots are used for determining if there is a relationship between
3106 two variables. Positive relations have a slope going from bottom left to top
3107 right; Negative relations have a slope going from top left to bottom right.
3108 """
3109 # Iterate over all cubes in cube or CubeList and plot.
3110 for cube_iter in iter_maybe(cube_x):
3111 # Check cubes are correct shape.
3112 cube_iter = check_single_cube(cube_iter)
3113 if cube_iter.ndim > 1:
3114 raise ValueError("cube_x must be 1D.")
3116 # Iterate over all cubes in cube or CubeList and plot.
3117 for cube_iter in iter_maybe(cube_y):
3118 # Check cubes are correct shape.
3119 cube_iter = check_single_cube(cube_iter)
3120 if cube_iter.ndim > 1:
3121 raise ValueError("cube_y must be 1D.")
3123 # Ensure we have a name for the plot file.
3124 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
3125 title = f"{recipe_title}"
3127 if filename is None:
3128 filename = slugify(recipe_title)
3130 # Add file extension.
3131 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
3133 # Do the actual plotting.
3134 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
3136 # Add list of plots to plot metadata.
3137 plot_index = _append_to_plot_index([plot_filename])
3139 # Make a page to display the plots.
3140 _make_plot_html_page(plot_index)
3142 return iris.cube.CubeList([cube_x, cube_y])
3145def vector_plot(
3146 cube_u: iris.cube.Cube,
3147 cube_v: iris.cube.Cube,
3148 filename: str | None = None,
3149 sequence_coordinate: str = "time",
3150 **kwargs,
3151) -> iris.cube.CubeList:
3152 """Plot a vector plot based on the input u and v components."""
3153 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
3155 # Cubes must have a matching sequence coordinate.
3156 try:
3157 # Check that the u and v cubes have the same sequence coordinate.
3158 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 3158 ↛ anywhereline 3158 didn't jump anywhere: it always raised an exception.
3159 raise ValueError("Coordinates do not match.")
3160 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
3161 raise ValueError(
3162 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
3163 ) from err
3165 # Create a plot for each value of the sequence coordinate.
3166 plot_index = []
3167 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
3168 for cube_u_slice, cube_v_slice in zip(
3169 cube_u.slices_over(sequence_coordinate),
3170 cube_v.slices_over(sequence_coordinate),
3171 strict=True,
3172 ):
3173 # Format the coordinate value in a unit appropriate way.
3174 seq_coord = cube_u_slice.coord(sequence_coordinate)
3175 plot_title, plot_filename = _set_title_and_filename(
3176 seq_coord, nplot, recipe_title, filename
3177 )
3179 # Do the actual plotting.
3180 _plot_and_save_vector_plot(
3181 cube_u_slice,
3182 cube_v_slice,
3183 filename=plot_filename,
3184 title=plot_title,
3185 method="pcolormesh",
3186 )
3187 plot_index.append(plot_filename)
3189 # Add list of plots to plot metadata.
3190 complete_plot_index = _append_to_plot_index(plot_index)
3192 # Make a page to display the plots.
3193 _make_plot_html_page(complete_plot_index)
3195 return iris.cube.CubeList([cube_u, cube_v])
3198def plot_histogram_series(
3199 cubes: iris.cube.Cube | iris.cube.CubeList,
3200 filename: str | None = None,
3201 sequence_coordinate: str = "time",
3202 stamp_coordinate: str = "realization",
3203 single_plot: bool = False,
3204 **kwargs,
3205) -> iris.cube.Cube | iris.cube.CubeList:
3206 """Plot a histogram plot for each vertical level provided.
3208 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
3209 is present then a sequence of plots will be produced using the time slider
3210 functionality to scroll through histograms against time. If a
3211 stamp_coordinate is present then postage stamp plots will be produced. If
3212 stamp_coordinate and single_plot is True, all postage stamp plots will be
3213 plotted in a single plot instead of separate postage stamp plots.
3215 Parameters
3216 ----------
3217 cubes: Cube | iris.cube.CubeList
3218 Iris cube or CubeList of the data to plot. It should have a single dimension other
3219 than the stamp coordinate.
3220 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3221 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3222 filename: str, optional
3223 Name of the plot to write, used as a prefix for plot sequences. Defaults
3224 to the recipe name.
3225 sequence_coordinate: str, optional
3226 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3227 This coordinate must exist in the cube and will be used for the time
3228 slider.
3229 stamp_coordinate: str, optional
3230 Coordinate about which to plot postage stamp plots. Defaults to
3231 ``"realization"``.
3232 single_plot: bool, optional
3233 If True, all postage stamp plots will be plotted in a single plot. If
3234 False, each postage stamp plot will be plotted separately. Is only valid
3235 if stamp_coordinate exists and has more than a single point.
3237 Returns
3238 -------
3239 iris.cube.Cube | iris.cube.CubeList
3240 The original Cube or CubeList (so further operations can be applied).
3241 Plotted data.
3243 Raises
3244 ------
3245 ValueError
3246 If the cube doesn't have the right dimensions.
3247 TypeError
3248 If the cube isn't a Cube or CubeList.
3249 """
3250 recipe_title = get_recipe_metadata().get("title", "Histogram")
3252 cubes = iter_maybe(cubes)
3254 # Internal plotting function.
3255 plotting_func = _plot_and_save_histogram_series
3257 num_models = get_num_models(cubes)
3259 validate_cube_shape(cubes, num_models)
3261 # If several histograms are plotted, check sequence_coordinate
3262 check_sequence_coordinate(cubes, sequence_coordinate)
3264 # Get axis minimum and maximum from levels information.
3265 # If no levels set, derive minima and maxima from data in CubeList.
3266 vmin, vmax = _set_axis_range(cubes)
3268 # Make postage stamp plots if stamp_coordinate exists and has more than a
3269 # single point. If single_plot is True:
3270 # -- all postage stamp plots will be plotted in a single plot instead of
3271 # separate postage stamp plots.
3272 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3273 # produced per single model only
3274 if num_models == 1:
3275 if ( 3275 ↛ 3279line 3275 didn't jump to line 3279 because the condition on line 3275 was never true
3276 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3277 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3278 ):
3279 if single_plot:
3280 plotting_func = (
3281 _plot_and_save_postage_stamps_in_single_plot_histogram_series
3282 )
3283 else:
3284 plotting_func = _plot_and_save_postage_stamp_histogram_series
3285 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3286 else:
3287 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3289 plot_index = []
3290 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3291 # Create a plot for each value of the sequence coordinate. Allowing for
3292 # multiple cubes in a CubeList to be plotted in the same plot for similar
3293 # sequence values. Passing a CubeList into the internal plotting function
3294 # for similar values of the sequence coordinate. cube_slice can be an
3295 # iris.cube.Cube or an iris.cube.CubeList.
3296 for cube_slice in cube_iterables:
3297 single_cube = cube_slice
3298 if isinstance(cube_slice, iris.cube.CubeList):
3299 single_cube = cube_slice[0]
3301 # Ensure valid stamp coordinate in cube dimensions
3302 if stamp_coordinate == "realization": 3302 ↛ 3305line 3302 didn't jump to line 3305 because the condition on line 3302 was always true
3303 stamp_coordinate = check_stamp_coordinate(single_cube)
3304 # Set plot titles and filename, based on sequence coordinate
3305 seq_coord = single_cube.coord(sequence_coordinate)
3306 # Use time coordinate in title and filename if single histogram output.
3307 if sequence_coordinate == "realization" and nplot == 1: 3307 ↛ 3308line 3307 didn't jump to line 3308 because the condition on line 3307 was never true
3308 seq_coord = single_cube.coord("time")
3309 # Use station name in title and filename if model vs obs comparison
3310 if sequence_coordinate == "station": 3310 ↛ 3311line 3310 didn't jump to line 3311 because the condition on line 3310 was never true
3311 seq_coord = single_cube.coord("Station_Name")
3313 plot_title, plot_filename = _set_title_and_filename(
3314 seq_coord, nplot, recipe_title, filename
3315 )
3317 # Do the actual plotting.
3318 plotting_func(
3319 cube_slice,
3320 filename=plot_filename,
3321 stamp_coordinate=stamp_coordinate,
3322 title=plot_title,
3323 vmin=vmin,
3324 vmax=vmax,
3325 )
3326 plot_index.append(plot_filename)
3328 # Add list of plots to plot metadata.
3329 complete_plot_index = _append_to_plot_index(plot_index)
3331 # Make a page to display the plots.
3332 _make_plot_html_page(complete_plot_index)
3334 return cubes
3337def plot_scatter_series(
3338 cubes: iris.cube.Cube | iris.cube.CubeList,
3339 filename: str | None = None,
3340 sequence_coordinate: str = "time",
3341 stamp_coordinate: str = "realization",
3342 hexbin: bool = False,
3343 **kwargs,
3344) -> iris.cube.Cube | iris.cube.CubeList:
3345 """Plot a scatter plot for each sequence coordinate provided.
3347 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
3348 is present then a sequence of plots will be produced using the time slider
3349 functionality to scroll through scatter against time. If a
3350 stamp_coordinate is present then postage stamp plots will be produced. If
3351 stamp_coordinate and single_plot is True, all postage stamp plots will be
3352 plotted in a single plot instead of separate postage stamp plots.
3354 Parameters
3355 ----------
3356 cubes: Cube | iris.cube.CubeList
3357 Iris cube or CubeList of the data to plot. It should have a single dimension other
3358 than the stamp coordinate.
3359 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
3360 We do not support different data such as temperature and humidity in the same CubeList for plotting.
3361 filename: str, optional
3362 Name of the plot to write, used as a prefix for plot sequences. Defaults
3363 to the recipe name.
3364 sequence_coordinate: str, optional
3365 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
3366 This coordinate must exist in the cube and will be used for the time
3367 slider.
3368 stamp_coordinate: str, optional
3369 Coordinate about which to plot postage stamp plots. Defaults to
3370 ``"realization"``.
3371 hexbin: bool, optional
3372 If True, generate hexbin comparison plot.
3373 If False, generate point-by-point scatter plot.
3375 Returns
3376 -------
3377 iris.cube.Cube | iris.cube.CubeList
3378 The original Cube or CubeList (so further operations can be applied).
3379 Plotted data.
3381 Raises
3382 ------
3383 ValueError
3384 If the cube doesn't have the right dimensions.
3385 TypeError
3386 If the cube isn't a Cube or CubeList.
3387 """
3388 recipe_title = get_recipe_metadata().get("title", "Scatter")
3390 cubes = iter_maybe(cubes)
3392 # Internal plotting function.
3393 plotting_func = _plot_and_save_scatter_series
3395 num_models = get_num_models(cubes)
3397 validate_cube_shape(cubes, num_models)
3399 check_sequence_coordinate(cubes, sequence_coordinate)
3401 vmin, vmax = _set_axis_range(cubes)
3403 # Require >1 models to compare on scatter plot
3404 if num_models > 1:
3405 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3406 else:
3407 raise ValueError(
3408 "Scatter plot series requires multiple number of models in input data."
3409 )
3411 plot_index = []
3412 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3413 # Create a plot for each value of the sequence coordinate. Allowing for
3414 # multiple cubes in a CubeList to be plotted in the same plot for similar
3415 # sequence values. Passing a CubeList into the internal plotting function
3416 # for similar values of the sequence coordinate. cube_slice can be an
3417 # iris.cube.Cube or an iris.cube.CubeList.
3418 for cube_slice in cube_iterables:
3419 single_cube = cube_slice
3420 if isinstance(cube_slice, iris.cube.CubeList): 3420 ↛ 3424line 3420 didn't jump to line 3424 because the condition on line 3420 was always true
3421 single_cube = cube_slice[0]
3423 # Ensure valid stamp coordinate in cube dimensions
3424 if stamp_coordinate == "realization": 3424 ↛ 3427line 3424 didn't jump to line 3427 because the condition on line 3424 was always true
3425 stamp_coordinate = check_stamp_coordinate(single_cube)
3426 # Set plot titles and filename, based on sequence coordinate
3427 seq_coord = single_cube.coord(sequence_coordinate)
3428 # Use time coordinate in title and filename if single histogram output.
3429 if sequence_coordinate == "realization" and nplot == 1:
3430 seq_coord = single_cube.coord("time")
3431 # Use station name in title and filename if model vs obs comparison
3432 if sequence_coordinate == "station":
3433 seq_coord = single_cube.coord("Station_Name")
3435 plot_title, plot_filename = _set_title_and_filename(
3436 seq_coord, nplot, recipe_title, filename
3437 )
3439 # Do the actual plotting.
3440 plotting_func(
3441 cube_slice,
3442 filename=plot_filename,
3443 stamp_coordinate=stamp_coordinate,
3444 title=plot_title,
3445 vmin=vmin,
3446 vmax=vmax,
3447 hexbin=hexbin,
3448 )
3449 plot_index.append(plot_filename)
3451 # Add list of plots to plot metadata.
3452 complete_plot_index = _append_to_plot_index(plot_index)
3454 # Make a page to display the plots.
3455 _make_plot_html_page(complete_plot_index)
3457 return cubes
3460def _plot_and_save_postage_stamp_power_spectrum_series(
3461 cubes: iris.cube.Cube,
3462 coords: list[iris.coords.Coord],
3463 stamp_coordinate: str,
3464 filename: str,
3465 title: str,
3466 series_coordinate: str | None = None,
3467 **kwargs,
3468):
3469 """Plot and save postage (ensemble members) stamps for a power spectrum series.
3471 Parameters
3472 ----------
3473 cubes: Cube or CubeList
3474 Cube or Cubelist of the power spectrum data.
3475 coords: list[Coord]
3476 Coordinates to plot on the x-axis, one per cube.
3477 stamp_coordinate: str
3478 Coordinate that becomes different plots.
3479 filename: str
3480 Filename of the plot to write.
3481 title: str
3482 Plot title.
3483 series_coordinate: str, optional
3484 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3486 """
3487 # Use the smallest square grid that will fit the members.
3488 grid_size = math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points)))
3490 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
3491 model_colors_map = get_model_colors_map(cubes)
3492 # ax = plt.gca()
3493 # Make a subplot for each member.
3494 for member, subplot in zip(
3495 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
3496 ):
3497 ax = plt.subplot(grid_size, grid_size, subplot)
3499 # Store min/max ranges.
3500 y_levels = []
3502 line_marker = None
3503 line_width = 1
3505 for cube in iter_maybe(member):
3506 xcoord = _select_series_coord(cube, series_coordinate)
3507 xname = xcoord.points
3509 yfield = cube.data # power spectrum
3510 label = None
3511 color = "black"
3512 if model_colors_map: 3512 ↛ 3513line 3512 didn't jump to line 3513 because the condition on line 3512 was never true
3513 label = cube.attributes.get("model_name")
3514 color = model_colors_map.get(label)
3516 if member.coord(stamp_coordinate).points == [0]:
3517 ax.plot(
3518 xname,
3519 yfield,
3520 color=color,
3521 marker=line_marker,
3522 ls="-",
3523 lw=line_width,
3524 label=f"{label} (control)"
3525 if len(cube.coord(stamp_coordinate).points) > 1
3526 else label,
3527 )
3528 # Label with member if part of an ensemble and not the control.
3529 else:
3530 ax.plot(
3531 xname,
3532 yfield,
3533 color=color,
3534 ls="-",
3535 lw=1.5,
3536 alpha=0.75,
3537 label=f"{label} (member)",
3538 )
3540 # Calculate the global min/max if multiple cubes are given.
3541 _, levels, _ = colorbar_map_levels(cube, axis="y")
3542 if levels is not None: 3542 ↛ 3543line 3542 didn't jump to line 3543 because the condition on line 3542 was never true
3543 y_levels.append(min(levels))
3544 y_levels.append(max(levels))
3546 # Add some labels and tweak the style.
3547 title = f"{title}"
3548 ax.set_title(title, fontsize=16)
3550 # Set appropriate x-axis label based on coordinate
3551 if series_coordinate == "wavelength" or ( 3551 ↛ 3554line 3551 didn't jump to line 3554 because the condition on line 3551 was never true
3552 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3553 ):
3554 ax.set_xlabel("Wavelength (km)", fontsize=14)
3555 elif series_coordinate == "physical_wavenumber" or ( 3555 ↛ 3560line 3555 didn't jump to line 3560 because the condition on line 3555 was always true
3556 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3557 ):
3558 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3559 else: # frequency or check units
3560 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3561 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3562 else:
3563 ax.set_xlabel("Wavenumber", fontsize=14)
3565 ax.set_ylabel("Power Spectral Density", fontsize=14)
3566 ax.tick_params(axis="both", labelsize=12)
3568 # Set log-log scale
3569 ax.set_xscale("log")
3570 ax.set_yscale("log")
3572 # Add gridlines
3573 ax.grid(linestyle="--", color="grey", linewidth=1)
3574 # Ientify unique labels for legend
3575 handles = list(
3576 {
3577 label: handle
3578 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3579 }.values()
3580 )
3581 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3583 ax = plt.gca()
3584 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3586 # Save plot.
3587 _save_close_figure(fig, "histogram postage stamp", filename)
3590def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3591 cubes: iris.cube.Cube,
3592 coords: list[iris.coords.Coord],
3593 stamp_coordinate: str,
3594 filename: str,
3595 title: str,
3596 series_coordinate: str | None = None,
3597 **kwargs,
3598):
3599 """Plot and save power spectra for ensemble members in single plot.
3601 Parameters
3602 ----------
3603 cubes: Cube or CubeList
3604 Cube or Cubelist of the power spectrum data.
3605 coords: list[Coord]
3606 Coordinates to plot on the x-axis, one per cube.
3607 stamp_coordinate: str
3608 Coordinate that becomes different plots.
3609 filename: str
3610 Filename of the plot to write.
3611 title: str
3612 Plot title.
3613 series_coordinate: str, optional
3614 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3616 """
3617 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3618 model_colors_map = get_model_colors_map(cubes)
3620 line_marker = None
3621 line_width = 1
3623 # Compute ensemble statistics to show spread
3624 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3625 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3626 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3628 xcoord_global = mean_cube.coord(series_coordinate)
3629 x_global = xcoord_global.points
3631 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3632 xcoord = _select_series_coord(member, series_coordinate)
3633 xname = xcoord.points
3635 yfield = member.data # power spectrum
3636 color = "black"
3637 if model_colors_map: 3637 ↛ 3641line 3637 didn't jump to line 3641 because the condition on line 3637 was always true
3638 label = member.attributes.get("model_name") if i == 0 else None
3639 color = model_colors_map.get(label)
3641 if member.coord(stamp_coordinate).points == [0]:
3642 ax.plot(
3643 xname,
3644 yfield,
3645 color=color,
3646 marker=line_marker,
3647 ls="-",
3648 lw=line_width,
3649 label=f"{label} (control)"
3650 if len(member.coord(stamp_coordinate).points) > 1
3651 else label,
3652 )
3653 # Label with member number if part of an ensemble and not the control.
3654 else:
3655 ax.plot(
3656 xname,
3657 yfield,
3658 color=color,
3659 ls="-",
3660 lw=1.5,
3661 alpha=0.75,
3662 label=label,
3663 )
3665 # Set appropriate x-axis label based on coordinate
3666 if series_coordinate == "wavelength" or ( 3666 ↛ 3669line 3666 didn't jump to line 3669 because the condition on line 3666 was never true
3667 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3668 ):
3669 ax.set_xlabel("Wavelength (km)", fontsize=14)
3670 elif series_coordinate == "physical_wavenumber" or ( 3670 ↛ 3675line 3670 didn't jump to line 3675 because the condition on line 3670 was always true
3671 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3672 ):
3673 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3674 else: # frequency or check units
3675 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3676 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3677 else:
3678 ax.set_xlabel("Wavenumber", fontsize=14)
3680 # Add ensemble spread shading
3681 ax.fill_between(
3682 x_global,
3683 min_cube.data,
3684 max_cube.data,
3685 color="grey",
3686 alpha=0.3,
3687 label="Ensemble spread",
3688 )
3690 # Add ensemble mean line
3691 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3693 ax.set_ylabel("Power Spectral Density", fontsize=14)
3694 ax.tick_params(axis="both", labelsize=12)
3696 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3697 # Set log-log scale
3698 ax.set_xscale("log")
3699 ax.set_yscale("log")
3701 # Add gridlines
3702 ax.grid(linestyle="--", color="grey", linewidth=1)
3703 # Identify unique labels for legend
3704 handles = list(
3705 {
3706 label: handle
3707 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3708 }.values()
3709 )
3710 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3712 # Figure title.
3713 ax.set_title(title, fontsize=16)
3715 # Save plot.
3716 _save_close_figure(fig, "power spectra postage stamp", filename)