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