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