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