Coverage for src/CSET/operators/plot.py: 54%
1024 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-10 13:27 +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.colors as mcolors
35import matplotlib.pyplot as plt
36import numpy as np
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._colormaps import (
50 colorbar_map_levels,
51 get_model_colors_map,
52)
53from CSET.operators._utils import (
54 check_sequence_coordinate,
55 check_single_cube,
56 check_stamp_coordinate,
57 fully_equalise_attributes,
58 get_cube_yxcoordname,
59 get_num_models,
60 is_transect,
61 slice_over_maybe,
62 validate_cube_shape,
63 validate_cubes_coords,
64)
65from CSET.operators.collapse import collapse
66from CSET.operators.misc import _extract_common_time_points
67from CSET.operators.regrid import regrid_onto_cube
69# Use a non-interactive plotting backend.
70mpl.use("agg")
73############################
74# Private helper functions #
75############################
78def _append_to_plot_index(plot_index: list) -> list:
79 """Add plots into the plot index, returning the complete plot index."""
80 with open("meta.json", "r+t", encoding="UTF-8") as fp:
81 fcntl.flock(fp, fcntl.LOCK_EX)
82 fp.seek(0)
83 meta = json.load(fp)
84 complete_plot_index = meta.get("plots", [])
85 complete_plot_index = complete_plot_index + plot_index
86 meta["plots"] = complete_plot_index
87 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
88 os.getenv("DO_CASE_AGGREGATION")
89 ):
90 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
91 fp.seek(0)
92 fp.truncate()
93 json.dump(meta, fp, indent=2)
94 return complete_plot_index
97def _make_plot_html_page(plots: list):
98 """Create a HTML page to display a plot image."""
99 # Debug check that plots actually contains some strings.
100 assert isinstance(plots[0], str)
102 # Load HTML template file.
103 operator_files = importlib.resources.files()
104 template_file = operator_files.joinpath("_plot_page_template.html")
106 # Get some metadata.
107 meta = get_recipe_metadata()
108 title = meta.get("title", "Untitled")
109 description = MarkdownIt().render(meta.get("description", "*No description.*"))
111 # Prepare template variables.
112 variables = {
113 "title": title,
114 "description": description,
115 "initial_plot": plots[0],
116 "plots": plots,
117 "title_slug": slugify(title),
118 }
120 # Render template.
121 html = render_file(template_file, **variables)
123 # Save completed HTML.
124 with open("index.html", "wt", encoding="UTF-8") as fp:
125 fp.write(html)
128def _setup_spatial_map(
129 cube: iris.cube.Cube,
130 figure,
131 cmap,
132 grid_size: tuple[int, int] | None = None,
133 subplot: int | None = None,
134):
135 """Define map projections, extent and add coastlines and borderlines for spatial plots.
137 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
138 is specified, and map extent defined based on the input data.
140 Parameters
141 ----------
142 cube: Cube
143 2 dimensional (lat and lon) Cube of the data to plot.
144 figure:
145 Matplotlib Figure object holding all plot elements.
146 cmap:
147 Matplotlib colormap.
148 grid_size: (int, int), optional
149 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
150 subplot: int, optional
151 Subplot index if multiple spatial subplots in figure.
153 Returns
154 -------
155 axes:
156 Matplotlib GeoAxes definition.
157 """
158 # Identify min/max plot bounds.
159 try:
160 lat_axis, lon_axis = get_cube_yxcoordname(cube)
161 x1 = np.min(cube.coord(lon_axis).points)
162 x2 = np.max(cube.coord(lon_axis).points)
163 y1 = np.min(cube.coord(lat_axis).points)
164 y2 = np.max(cube.coord(lat_axis).points)
166 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
167 if np.abs(x2 - x1) > 180.0:
168 x1 = x1 - 180.0
169 x2 = x2 - 180.0
170 logging.debug("Adjusting plot bounds to fit global extent.")
172 # Consider map projection orientation.
173 # Adapting orientation enables plotting across international dateline.
174 # Users can adapt the default central_longitude if alternative projections views.
175 if x2 > 180.0 or x1 < -180.0:
176 central_longitude = 180.0
177 else:
178 central_longitude = 0.0
180 # Define spatial map projection.
181 coord_system = cube.coord(lat_axis).coord_system
182 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
183 # Define rotated pole map projection for rotated pole inputs.
184 projection = ccrs.RotatedPole(
185 pole_longitude=coord_system.grid_north_pole_longitude,
186 pole_latitude=coord_system.grid_north_pole_latitude,
187 central_rotated_longitude=central_longitude,
188 )
189 crs = projection
190 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 190 ↛ 192line 190 didn't jump to line 192 because the condition on line 190 was never true
191 # Define Transverse Mercator projection for TM inputs.
192 projection = ccrs.TransverseMercator(
193 central_longitude=coord_system.longitude_of_central_meridian,
194 central_latitude=coord_system.latitude_of_projection_origin,
195 false_easting=coord_system.false_easting,
196 false_northing=coord_system.false_northing,
197 scale_factor=coord_system.scale_factor_at_central_meridian,
198 )
199 crs = projection
200 else:
201 # Define regular map projection for non-rotated pole inputs.
202 # Alternatives might include e.g. for global model outputs:
203 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
204 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
205 projection = ccrs.PlateCarree(central_longitude=central_longitude)
206 crs = ccrs.PlateCarree()
208 # Define axes for plot (or subplot) with required map projection.
209 if subplot is not None:
210 axes = figure.add_subplot(
211 grid_size[0], grid_size[1], subplot, projection=projection
212 )
213 else:
214 axes = figure.add_subplot(projection=projection)
216 # Add coastlines and borderlines if cube contains x and y map coordinates.
217 # Avoid adding lines for masked data or specific fixed ancillary spatial plots.
218 if iris.util.is_masked(cube.data) or any( 218 ↛ 221line 218 didn't jump to line 221 because the condition on line 218 was never true
219 name in cube.name() for name in ["land_", "orography", "altitude"]
220 ):
221 pass
222 else:
223 if cmap.name in ["viridis", "Greys"]:
224 coastcol = "magenta"
225 else:
226 coastcol = "black"
227 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
228 axes.coastlines(resolution="10m", color=coastcol)
229 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
231 # Add gridlines.
232 gl = axes.gridlines(
233 alpha=0.3,
234 draw_labels=True,
235 dms=False,
236 x_inline=False,
237 y_inline=False,
238 )
239 gl.top_labels = False
240 gl.right_labels = False
241 if subplot:
242 gl.bottom_labels = False
243 gl.left_labels = False
244 if subplot % grid_size[1] == 1:
245 gl.left_labels = True
246 if subplot > ((grid_size[0] - 1) * grid_size[1]): 246 ↛ 251line 246 didn't jump to line 251 because the condition on line 246 was always true
247 gl.bottom_labels = True
249 # If is lat/lon spatial map, fix extent to keep plot tight.
250 # Specifying crs within set_extent helps ensure only data region is shown.
251 if isinstance(coord_system, iris.coord_systems.GeogCS):
252 axes.set_extent([x1, x2, y1, y2], crs=crs)
254 except ValueError:
255 # Skip if not both x and y map coordinates.
256 axes = figure.gca()
257 pass
259 return axes
262def _get_plot_resolution() -> int:
263 """Get resolution of rasterised plots in pixels per inch."""
264 return get_recipe_metadata().get("plot_resolution", 100)
267def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
268 """Return title and filename based on start and end points or bounds."""
269 if use_bounds and seq_coord.has_bounds():
270 vals = seq_coord.bounds.flatten()
271 else:
272 vals = seq_coord.points
273 start = seq_coord.units.title(vals[0])
274 end = seq_coord.units.title(vals[-1])
276 if start == end:
277 sequence_title = f"\n [{start}]"
278 sequence_fname = f"_{filename_slugify(start)}"
279 else:
280 sequence_title = f"\n [{start} to {end}]"
281 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
283 if seq_coord.units == "unknown": 283 ↛ 285line 283 didn't jump to line 285 because the condition on line 283 was never true
284 # remove the "unknown" in title and filename strings if unit is unknown
285 sequence_title = sequence_title.replace("unknown", "")
286 sequence_fname = sequence_fname.replace("unknown", "")
288 # Do not include time if coord set to zero.
289 if (
290 seq_coord.units == "hours since 0001-01-01 00:00:00"
291 and vals[0] == 0
292 and vals[-1] == 0
293 ):
294 sequence_title = ""
295 sequence_fname = ""
297 return sequence_title, sequence_fname
300def _set_title_and_filename(
301 seq_coord: iris.coords.Coord,
302 nplot: int,
303 recipe_title: str,
304 filename: str,
305):
306 """Set plot title and filename based on cube coordinate.
308 Parameters
309 ----------
310 sequence_coordinate: iris.coords.Coord
311 Coordinate about which to make a plot sequence.
312 nplot: int
313 Number of output plots to generate - controls title/naming.
314 recipe_title: str
315 Default plot title, potentially to update.
316 filename: str
317 Input plot filename, potentially to update.
319 Returns
320 -------
321 plot_title: str
322 Output formatted plot title string, based on plotted data.
323 plot_filename: str
324 Output formatted plot filename string.
325 """
326 ndim = seq_coord.ndim
327 npoints = np.size(seq_coord.points)
328 sequence_title = ""
329 sequence_fname = ""
331 # Case 1: Multiple dimension sequence input - list number of aggregated cases
332 # (e.g. aggregation histogram plots)
333 if ndim > 1:
334 ncase = np.shape(seq_coord)[0]
335 sequence_title = f"\n [{ncase} cases]"
336 sequence_fname = f"_{ncase}cases"
338 # Case 2: Single dimension input
339 else:
340 # Single sequence point
341 if npoints == 1:
342 if nplot > 1:
343 # Default labels for sequence inputs
344 sequence_value = seq_coord.units.title(seq_coord.points[0])
345 sequence_title = f"\n [{sequence_value}]"
346 sequence_fname = f"_{filename_slugify(sequence_value)}"
347 else:
348 # Aggregated attribute available where input collapsed over aggregation
349 try:
350 ncase = seq_coord.attributes["number_reference_times"]
351 sequence_title = f"\n [{ncase} cases]"
352 sequence_fname = f"_{ncase}cases"
353 except KeyError:
354 sequence_title, sequence_fname = _get_start_end_strings(
355 seq_coord, use_bounds=seq_coord.has_bounds()
356 )
357 # Multiple sequence (e.g. time) points
358 else:
359 sequence_title, sequence_fname = _get_start_end_strings(
360 seq_coord, use_bounds=False
361 )
363 # Set plot title and filename
364 plot_title = f"{recipe_title}{sequence_title}"
366 # Set plot filename, defaulting to user input if provided.
367 if filename is None:
368 filename = slugify(recipe_title)
369 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
370 else:
371 if nplot > 1:
372 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
373 else:
374 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
376 return plot_title, plot_filename
379def _select_series_coord(cube, series_coordinate):
380 """Determine the grid coordinates to use to calculate grid spacing."""
381 spacing_coordinates = ("frequency", "physical_wavenumber", "wavelength")
382 if series_coordinate in spacing_coordinates: 382 ↛ 388line 382 didn't jump to line 388 because the condition on line 382 was always true
383 # Try the requested coordinate first then the fallbacks in order.
384 fallbacks = [series_coordinate] + [
385 c for c in spacing_coordinates if c != series_coordinate
386 ]
387 else:
388 fallbacks = {series_coordinate}
390 # Try each possible coordinate.
391 for coord in fallbacks:
392 try:
393 return cube.coord(coord)
394 except iris.exceptions.CoordinateNotFoundError:
395 logging.debug("Coordinate %s not found.", coord)
397 # If we get here, none of the fallback options were found.
398 raise iris.exceptions.CoordinateNotFoundError(
399 f"No valid coordinate found for '{series_coordinate}' "
400 f"or fallback options {fallbacks}"
401 )
404def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
405 """Control postage stamp plot output titles based on stamp coordinate."""
406 if stamp_coord.name() == "realization":
407 mtitle = "Member"
408 else:
409 mtitle = stamp_coord.name().capitalize()
411 if stamp_coord.name() == "time":
412 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
413 else:
414 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
416 return mtitle
419def _set_axis_range(cubes):
420 """Get minimum and maximum from levels information."""
421 levels = None
422 for cube in cubes:
423 # First check if user-specified "auto" range variable.
424 # This maintains the value of levels as None, so proceed.
425 _, levels, _ = colorbar_map_levels(cube, axis="y")
426 if levels is None:
427 break
428 # If levels is changed, recheck to use the vmin,vmax or
429 # levels-based ranges for histogram plots.
430 _, levels, _ = colorbar_map_levels(cube)
431 logging.debug("levels: %s", levels)
432 if levels is not None:
433 vmin = min(levels)
434 vmax = max(levels)
435 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
436 break
438 if levels is None:
439 vmin = min(cb.data.min() for cb in cubes)
440 vmax = max(cb.data.max() for cb in cubes)
442 return vmin, vmax
445def _find_matched_slices(cubes, sequence_coordinate):
446 """Identify matched cubes in CubeList by sequence_coordinate values.
448 Ensures common points are compared for multiple cube inputs.
449 """
450 all_points = sorted(
451 set(
452 itertools.chain.from_iterable(
453 cb.coord(sequence_coordinate).points for cb in cubes
454 )
455 )
456 )
457 all_slices = list(
458 itertools.chain.from_iterable(
459 cb.slices_over(sequence_coordinate) for cb in cubes
460 )
461 )
462 # Matched slices (matched by seq coord point; it may happen that
463 # evaluated models do not cover the same seq coord range, hence matching
464 # necessary)
465 cube_iterables = [
466 iris.cube.CubeList(
467 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
468 )
469 for point in all_points
470 ]
472 return cube_iterables
475def _plot_and_save_spatial_plot(
476 cube: iris.cube.Cube,
477 filename: str,
478 title: str,
479 method: Literal["contourf", "pcolormesh"],
480 overlay_cube: iris.cube.Cube | None = None,
481 contour_cube: iris.cube.Cube | None = None,
482 **kwargs,
483):
484 """Plot and save a spatial plot.
486 Parameters
487 ----------
488 cube: Cube
489 2 dimensional (lat and lon) Cube of the data to plot.
490 filename: str
491 Filename of the plot to write.
492 title: str
493 Plot title.
494 method: "contourf" | "pcolormesh"
495 The plotting method to use.
496 overlay_cube: Cube, optional
497 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
498 contour_cube: Cube, optional
499 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
500 """
501 # Setup plot details, size, resolution, etc.
502 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
504 # Specify the color bar
505 cmap, levels, norm = colorbar_map_levels(cube)
507 # If overplotting, set required colorbars
508 if overlay_cube:
509 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
510 if contour_cube:
511 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
513 # Setup plot map projection, extent and coastlines and borderlines.
514 axes = _setup_spatial_map(cube, fig, cmap)
516 # Set colorscale bounds
517 try:
518 vmin = min(levels)
519 vmax = max(levels)
520 except TypeError:
521 vmin, vmax = None, None
522 # Ensure to use norm and not vmin/vmax if levels are defined.
523 if norm is not None:
524 vmin = None
525 vmax = None
526 logging.debug("Plotting using defined levels.")
528 # Plot the field.
529 if method == "contourf":
530 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
531 elif method == "pcolormesh":
532 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
533 else:
534 raise ValueError(f"Unknown plotting method: {method}")
536 # Overplot overlay field, if required
537 if overlay_cube:
538 try:
539 over_vmin = min(over_levels)
540 over_vmax = max(over_levels)
541 except TypeError:
542 over_vmin, over_vmax = None, None
543 if over_norm is not None: 543 ↛ 544line 543 didn't jump to line 544 because the condition on line 543 was never true
544 over_vmin = None
545 over_vmax = None
546 overlay = iplt.pcolormesh(
547 overlay_cube,
548 cmap=over_cmap,
549 norm=over_norm,
550 alpha=0.8,
551 vmin=over_vmin,
552 vmax=over_vmax,
553 )
554 # Overplot contour field, if required, with contour labelling.
555 if contour_cube:
556 contour = iplt.contour(
557 contour_cube,
558 colors="darkgray",
559 levels=cntr_levels,
560 norm=cntr_norm,
561 alpha=0.5,
562 linestyles="--",
563 linewidths=1,
564 )
565 plt.clabel(contour)
567 # Check to see if transect, and if so, adjust y axis.
568 if is_transect(cube):
569 if "pressure" in [coord.name() for coord in cube.coords()]:
570 axes.invert_yaxis()
571 axes.set_yscale("log")
572 axes.set_ylim(1100, 100)
573 # If both model_level_number and level_height exists, iplt can construct
574 # plot as a function of height above orography (NOT sea level).
575 elif {"model_level_number", "level_height"}.issubset( 575 ↛ 580line 575 didn't jump to line 580 because the condition on line 575 was always true
576 {coord.name() for coord in cube.coords()}
577 ):
578 axes.set_yscale("log")
580 axes.set_title(
581 f"{title}\n"
582 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
583 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
584 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
585 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
586 fontsize=16,
587 )
589 # Inset code
590 axins = inset_axes(
591 axes,
592 width="20%",
593 height="20%",
594 loc="upper right",
595 axes_class=GeoAxes,
596 axes_kwargs={"map_projection": ccrs.PlateCarree()},
597 )
599 # Slightly transparent to reduce plot blocking.
600 axins.patch.set_alpha(0.4)
602 axins.coastlines(resolution="50m")
603 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
605 SLat, SLon, ELat, ELon = (
606 float(coord) for coord in cube.attributes["transect_coords"].split("_")
607 )
609 # Draw line between them
610 axins.plot(
611 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
612 )
614 # Plot points (note: lon, lat order for Cartopy)
615 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
616 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
618 lon_min, lon_max = sorted([SLon, ELon])
619 lat_min, lat_max = sorted([SLat, ELat])
621 # Midpoints
622 lon_mid = (lon_min + lon_max) / 2
623 lat_mid = (lat_min + lat_max) / 2
625 # Maximum half-range
626 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
627 if half_range == 0: # points identical → provide small default 627 ↛ 631line 627 didn't jump to line 631 because the condition on line 627 was always true
628 half_range = 1
630 # Set square extent
631 axins.set_extent(
632 [
633 lon_mid - half_range,
634 lon_mid + half_range,
635 lat_mid - half_range,
636 lat_mid + half_range,
637 ],
638 crs=ccrs.PlateCarree(),
639 )
641 # Ensure square aspect
642 axins.set_aspect("equal")
644 else:
645 # Add title.
646 axes.set_title(title, fontsize=16)
648 # Adjust padding if spatial plot or transect
649 if is_transect(cube):
650 yinfopad = -0.1
651 ycbarpad = 0.1
652 else:
653 yinfopad = 0.01
654 ycbarpad = 0.042
656 # Add watermark with min/max/mean. Currently not user togglable.
657 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
658 axes.annotate(
659 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
660 xy=(0.025, yinfopad),
661 xycoords="axes fraction",
662 xytext=(-5, 5),
663 textcoords="offset points",
664 ha="left",
665 va="bottom",
666 size=11,
667 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
668 )
670 # Add secondary colour bar for overlay_cube field if required.
671 if overlay_cube:
672 cbarB = fig.colorbar(
673 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
674 )
675 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
676 # add ticks and tick_labels for every levels if less than 20 levels exist
677 if over_levels is not None and len(over_levels) < 20: 677 ↛ 678line 677 didn't jump to line 678 because the condition on line 677 was never true
678 cbarB.set_ticks(over_levels)
679 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
680 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
681 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
682 logging.debug("Set secondary colorbar ticks and labels.")
684 # Add main colour bar.
685 cbar = fig.colorbar(
686 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
687 )
689 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
690 # add ticks and tick_labels for every levels if less than 20 levels exist
691 if levels is not None and len(levels) < 20:
692 cbar.set_ticks(levels)
693 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
694 if "rainfall" or "snowfall" or "visibility" in cube.name(): 694 ↛ 696line 694 didn't jump to line 696 because the condition on line 694 was always true
695 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
696 logging.debug("Set colorbar ticks and labels.")
698 # Save plot.
699 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
700 logging.info("Saved spatial plot to %s", filename)
701 plt.close(fig)
704def plot_dfss_contour(
705 cube: iris.cube.Cube | iris.cube.CubeList,
706 variable: str = None,
707) -> iris.cube.Cube | iris.cube.CubeList:
708 """Plot a scatter plot between two variables.
710 Both cubes must be 1D.
712 Parameters
713 ----------
714 cube: Cube | CubeList
715 1 dimensional Cube of the data to plot on y-axis.
716 filename: str, optional
717 Filename of the plot to write.
718 variable: str, optional
719 which cube variable to plot
721 Returns
722 -------
723 cubes: Cube
724 Cube of the original cubes for further processing.
726 Notes
727 -----
728 Makes a filled contour plot for dFSS/eFSS with neighbourhood lengths on the y-axis
729 and forecast lead time on the x-axis.
731 Adds a contour line at the 0.5 contour.
732 """
733 cube_copy = cube
734 if type(cube) is iris.cube.CubeList:
735 if not variable:
736 logging.Warning(
737 "CubeList given, but variable not specified. Defaulting to first cube."
738 )
739 cube = cube[0]
740 else:
741 cube = cube.extract(variable)[0]
743 recipe_title = get_recipe_metadata().get("title", "Untitled")
745 title = cube.name()
747 if cube.attributes.locals["method"] == "centile":
748 method = cube.attributes.locals["method"]
749 centile = cube.attributes.locals["centile"]
750 filename = f"{cube.name()}_{method}_{centile}.png"
751 plot_title = f"{recipe_title} \n {title} \n method={method} | centile={centile}"
753 elif cube.attributes.locals["method"] == "threshold":
754 method = cube.attributes.locals["method"]
755 threshold = cube.attributes.locals["threshold"]
756 filename = f"{cube.name()}_{method}_{threshold}.png"
757 plot_title = (
758 f"{recipe_title} \n {title} \n method={method} | threshold={threshold}"
759 )
761 cmap = plt.colormaps["viridis"]
762 levels = np.linspace(0, 1, 11)
763 norm = mcolors.BoundaryNorm(levels, cmap.N)
765 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
767 # set the contour colour for the 0.5 contour line for dfss only
769 plot = iplt.contourf(cube, cmap=cmap, norm=norm, levels=levels)
771 if cube.name() == "dfss":
772 colors = ["white"]
773 cmap_0p5_contour = mcolors.ListedColormap(colors)
775 iplt.contour(
776 cube, cmap=cmap_0p5_contour, norm=norm, levels=[0.5], linestyles="dashed"
777 )
779 plt.xlabel(cube.dim_coords[0].name(), fontsize=14)
780 plt.ylabel(cube.dim_coords[1].name(), fontsize=14)
782 cbar = fig.colorbar(
783 plot, orientation="horizontal", location="bottom", pad=0.1, shrink=0.7
784 )
786 cbar.set_label(label=f"{cube.name()}", size=14)
788 # Overall figure title.
790 fig.suptitle(plot_title, fontsize=16)
792 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
794 logging.info("Saved contour plot", filename)
795 plt.close(fig)
797 # Add list of plots to plot metadata.
798 plot_index = _append_to_plot_index([filename])
800 # Make a page to display the plots.
801 _make_plot_html_page(plot_index)
803 return cube_copy
806def _plot_and_save_postage_stamp_spatial_plot(
807 cube: iris.cube.Cube,
808 filename: str,
809 stamp_coordinate: str,
810 title: str,
811 method: Literal["contourf", "pcolormesh"],
812 overlay_cube: iris.cube.Cube | None = None,
813 contour_cube: iris.cube.Cube | None = None,
814 **kwargs,
815):
816 """Plot postage stamp spatial plots from an ensemble.
818 Parameters
819 ----------
820 cube: Cube
821 Iris cube of data to be plotted. It must have the stamp coordinate.
822 filename: str
823 Filename of the plot to write.
824 stamp_coordinate: str
825 Coordinate that becomes different plots.
826 method: "contourf" | "pcolormesh"
827 The plotting method to use.
828 overlay_cube: Cube, optional
829 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
830 contour_cube: Cube, optional
831 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
833 Raises
834 ------
835 ValueError
836 If the cube doesn't have the right dimensions.
837 """
838 # Use the smallest square grid that will fit the members.
839 nmember = len(cube.coord(stamp_coordinate).points)
840 grid_rows = int(math.sqrt(nmember))
841 grid_size = math.ceil(nmember / grid_rows)
843 fig = plt.figure(
844 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
845 )
847 # Specify the color bar
848 cmap, levels, norm = colorbar_map_levels(cube)
849 # If overplotting, set required colorbars
850 if overlay_cube: 850 ↛ 851line 850 didn't jump to line 851 because the condition on line 850 was never true
851 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
852 if contour_cube: 852 ↛ 853line 852 didn't jump to line 853 because the condition on line 852 was never true
853 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
855 # Make a subplot for each member.
856 for member, subplot in zip(
857 cube.slices_over(stamp_coordinate),
858 range(1, grid_size * grid_rows + 1),
859 strict=False,
860 ):
861 # Setup subplot map projection, extent and coastlines and borderlines.
862 axes = _setup_spatial_map(
863 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
864 )
865 if method == "contourf":
866 # Filled contour plot of the field.
867 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
868 elif method == "pcolormesh":
869 if levels is not None:
870 vmin = min(levels)
871 vmax = max(levels)
872 else:
873 raise TypeError("Unknown vmin and vmax range.")
874 vmin, vmax = None, None
875 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
876 # if levels are defined.
877 if norm is not None: 877 ↛ 878line 877 didn't jump to line 878 because the condition on line 877 was never true
878 vmin = None
879 vmax = None
880 # pcolormesh plot of the field.
881 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
882 else:
883 raise ValueError(f"Unknown plotting method: {method}")
885 # Overplot overlay field, if required
886 if overlay_cube: 886 ↛ 887line 886 didn't jump to line 887 because the condition on line 886 was never true
887 try:
888 over_vmin = min(over_levels)
889 over_vmax = max(over_levels)
890 except TypeError:
891 over_vmin, over_vmax = None, None
892 if over_norm is not None:
893 over_vmin = None
894 over_vmax = None
895 iplt.pcolormesh(
896 overlay_cube[member.coord(stamp_coordinate).points[0]],
897 cmap=over_cmap,
898 norm=over_norm,
899 alpha=0.6,
900 vmin=over_vmin,
901 vmax=over_vmax,
902 )
903 # Overplot contour field, if required
904 if contour_cube: 904 ↛ 905line 904 didn't jump to line 905 because the condition on line 904 was never true
905 iplt.contour(
906 contour_cube[member.coord(stamp_coordinate).points[0]],
907 colors="darkgray",
908 levels=cntr_levels,
909 norm=cntr_norm,
910 alpha=0.6,
911 linestyles="--",
912 linewidths=1,
913 )
914 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
915 axes.set_title(f"{mtitle}")
917 # Put the shared colorbar in its own axes.
918 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
919 colorbar = fig.colorbar(
920 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
921 )
922 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
924 # Overall figure title.
925 fig.suptitle(title, fontsize=16)
927 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
928 logging.info("Saved contour postage stamp plot to %s", filename)
929 plt.close(fig)
932def _plot_and_save_line_series(
933 cubes: iris.cube.CubeList,
934 coords: list[iris.coords.Coord],
935 filename: str,
936 title: str,
937 ensemble_coord: str = None,
938 sequence_coord: str = None,
939 **kwargs,
940):
941 """Plot and save a 1D line series.
943 Parameters
944 ----------
945 cubes: Cube or CubeList
946 Cube or CubeList containing the cubes to plot on the y-axis.
947 coords: list[Coord]
948 Coordinates to plot on the x-axis, one per cube.
949 ensemble_coord: str
950 Ensemble coordinate in the cube.
951 sequence_coord str
952 sequence coordinate for plotting, needed only if not an ensemble.
953 filename: str
954 Filename of the plot to write.
955 title: str
956 Plot title.
957 """
958 if (ensemble_coord is None) == (sequence_coord is None): 958 ↛ 959line 958 didn't jump to line 959 because the condition on line 958 was never true
959 raise ValueError(
960 "Exactly one of ensemble_coord and sequence_coord must be provided"
961 )
963 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
965 model_colors_map = get_model_colors_map(cubes)
967 # Store min/max ranges.
968 y_levels = []
970 # Check match-up across sequence coords gives consistent sizes
971 validate_cubes_coords(cubes, coords)
973 for cube, coord in zip(cubes, coords, strict=True):
974 label = None
975 color = "black"
976 if model_colors_map:
977 label = cube.attributes.get("model_name")
978 color = model_colors_map.get(label)
979 if ensemble_coord is not None: 979 ↛ 1005line 979 didn't jump to line 1005 because the condition on line 979 was always true
980 for cube_slice in cube.slices_over(ensemble_coord):
981 # Label with (control) if part of an ensemble or not otherwise.
982 if cube_slice.coord(ensemble_coord).points == [0]: 982 ↛ 996line 982 didn't jump to line 996 because the condition on line 982 was always true
983 iplt.plot(
984 coord,
985 cube_slice,
986 color=color,
987 marker="o",
988 ls="-",
989 lw=3,
990 label=f"{label} (control)"
991 if len(cube.coord(ensemble_coord).points) > 1
992 else label,
993 )
994 # Label with (perturbed) if part of an ensemble and not the control.
995 else:
996 iplt.plot(
997 coord,
998 cube_slice,
999 color=color,
1000 ls="-",
1001 lw=1.5,
1002 alpha=0.75,
1003 label=f"{label} (member)",
1004 )
1005 if sequence_coord is not None: 1005 ↛ 1006line 1005 didn't jump to line 1006 because the condition on line 1005 was never true
1006 for cube_slice in cube.slices_over(sequence_coord):
1007 iplt.plot(coord, cube_slice, color=color, ls="-", lw=1.5, alpha=0.75)
1009 # Calculate the global min/max if multiple cubes are given.
1010 _, levels, _ = colorbar_map_levels(cube, axis="y")
1011 if levels is not None: 1011 ↛ 1012line 1011 didn't jump to line 1012 because the condition on line 1011 was never true
1012 y_levels.append(min(levels))
1013 y_levels.append(max(levels))
1015 # Get the current axes.
1016 ax = plt.gca()
1018 # Add some labels and tweak the style.
1019 # check if cubes[0] works for single cube if not CubeList
1020 if coords[0].name() == "time": 1020 ↛ 1021line 1020 didn't jump to line 1021 because the condition on line 1020 was never true
1021 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1022 else:
1023 if coords[0].units == "unknown": 1023 ↛ 1024line 1023 didn't jump to line 1024 because the condition on line 1023 was never true
1024 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
1025 else:
1026 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
1028 if cubes[0].units == "unknown": 1028 ↛ 1031line 1028 didn't jump to line 1031 because the condition on line 1028 was always true
1029 ax.set_ylabel(f"{cubes[0].name()}", fontsize=14)
1030 else:
1031 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
1033 ax.set_title(title, fontsize=16)
1035 ax.ticklabel_format(axis="y", useOffset=False)
1036 ax.tick_params(axis="x", labelrotation=15)
1037 ax.tick_params(axis="both", labelsize=12)
1039 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1040 if y_levels: 1040 ↛ 1041line 1040 didn't jump to line 1041 because the condition on line 1040 was never true
1041 ax.set_ylim(min(y_levels), max(y_levels))
1042 # Add zero line.
1043 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
1044 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1045 logging.debug(
1046 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1047 )
1048 else:
1049 ax.autoscale()
1051 # Add gridlines
1052 ax.grid(linestyle="--", color="grey", linewidth=1)
1053 # Ientify unique labels for legend
1054 handles = list(
1055 {
1056 label: handle
1057 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1058 }.values()
1059 )
1060 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1062 # Save plot.
1063 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1064 logging.info("Saved line plot to %s", filename)
1065 plt.close(fig)
1068def _plot_and_save_line_power_spectrum_series(
1069 cubes: iris.cube.Cube | iris.cube.CubeList,
1070 coords: list[iris.coords.Coord],
1071 ensemble_coord: str,
1072 filename: str,
1073 title: str,
1074 series_coordinate: str,
1075 **kwargs,
1076):
1077 """Plot and save a 1D line series.
1079 Parameters
1080 ----------
1081 cubes: Cube or CubeList
1082 Cube or CubeList containing the cubes to plot on the y-axis.
1083 coords: list[Coord]
1084 Coordinates to plot on the x-axis, one per cube.
1085 ensemble_coord: str
1086 Ensemble coordinate in the cube.
1087 filename: str
1088 Filename of the plot to write.
1089 title: str
1090 Plot title.
1091 series_coordinate: str
1092 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1093 """
1094 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1095 model_colors_map = get_model_colors_map(cubes)
1096 ax = plt.gca()
1098 # Store min/max ranges.
1099 y_levels = []
1101 line_marker = None
1102 line_width = 1
1104 for cube in iter_maybe(cubes):
1105 # next 2 lines replace chunk of code.
1106 xcoord = _select_series_coord(cube, series_coordinate)
1107 xname = xcoord.points
1109 yfield = cube.data # power spectrum
1110 label = None
1111 color = "black"
1112 if model_colors_map: 1112 ↛ 1115line 1112 didn't jump to line 1115 because the condition on line 1112 was always true
1113 label = cube.attributes.get("model_name")
1114 color = model_colors_map.get(label)
1115 for cube_slice in cube.slices_over(ensemble_coord):
1116 # Label with (control) if part of an ensemble or not otherwise.
1117 if cube_slice.coord(ensemble_coord).points == [0]: 1117 ↛ 1131line 1117 didn't jump to line 1131 because the condition on line 1117 was always true
1118 ax.plot(
1119 xname,
1120 yfield,
1121 color=color,
1122 marker=line_marker,
1123 ls="-",
1124 lw=line_width,
1125 label=f"{label} (control)"
1126 if len(cube.coord(ensemble_coord).points) > 1
1127 else label,
1128 )
1129 # Label with (perturbed) if part of an ensemble and not the control.
1130 else:
1131 ax.plot(
1132 xname,
1133 yfield,
1134 color=color,
1135 ls="-",
1136 lw=1.5,
1137 alpha=0.75,
1138 label=f"{label} (member)",
1139 )
1141 # Calculate the global min/max if multiple cubes are given.
1142 _, levels, _ = colorbar_map_levels(cube, axis="y")
1143 if levels is not None: 1143 ↛ 1144line 1143 didn't jump to line 1144 because the condition on line 1143 was never true
1144 y_levels.append(min(levels))
1145 y_levels.append(max(levels))
1147 # Add some labels and tweak the style.
1149 title = f"{title}"
1150 ax.set_title(title, fontsize=16)
1152 # Set appropriate x-axis label based on coordinate
1153 if series_coordinate == "wavelength" or ( 1153 ↛ 1156line 1153 didn't jump to line 1156 because the condition on line 1153 was never true
1154 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1155 ):
1156 ax.set_xlabel("Wavelength (km)", fontsize=14)
1157 elif series_coordinate == "physical_wavenumber" or ( 1157 ↛ 1160line 1157 didn't jump to line 1160 because the condition on line 1157 was never true
1158 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1159 ):
1160 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1161 else: # frequency or check units
1162 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1162 ↛ 1163line 1162 didn't jump to line 1163 because the condition on line 1162 was never true
1163 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1164 else:
1165 ax.set_xlabel("Wavenumber", fontsize=14)
1167 ax.set_ylabel("Power Spectral Density", fontsize=14)
1168 ax.tick_params(axis="both", labelsize=12)
1170 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1172 # Set log-log scale
1173 ax.set_xscale("log")
1174 ax.set_yscale("log")
1176 # Add gridlines
1177 ax.grid(linestyle="--", color="grey", linewidth=1)
1178 # Ientify unique labels for legend
1179 handles = list(
1180 {
1181 label: handle
1182 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1183 }.values()
1184 )
1185 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1187 # Save plot.
1188 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1189 logging.info("Saved line plot to %s", filename)
1190 plt.close(fig)
1193def _plot_and_save_vertical_line_series(
1194 cubes: iris.cube.CubeList,
1195 coords: list[iris.coords.Coord],
1196 ensemble_coord: str,
1197 filename: str,
1198 series_coordinate: str,
1199 title: str,
1200 vmin: float,
1201 vmax: float,
1202 **kwargs,
1203):
1204 """Plot and save a 1D line series in vertical.
1206 Parameters
1207 ----------
1208 cubes: CubeList
1209 1 dimensional Cube or CubeList of the data to plot on x-axis.
1210 coord: list[Coord]
1211 Coordinates to plot on the y-axis, one per cube.
1212 ensemble_coord: str
1213 Ensemble coordinate in the cube.
1214 filename: str
1215 Filename of the plot to write.
1216 series_coordinate: str
1217 Coordinate to use as vertical axis.
1218 title: str
1219 Plot title.
1220 vmin: float
1221 Minimum value for the x-axis.
1222 vmax: float
1223 Maximum value for the x-axis.
1224 """
1225 # plot the vertical pressure axis using log scale
1226 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1228 model_colors_map = get_model_colors_map(cubes)
1230 # Check match-up across sequence coords gives consistent sizes
1231 validate_cubes_coords(cubes, coords)
1233 for cube, coord in zip(cubes, coords, strict=True):
1234 label = None
1235 color = "black"
1236 if model_colors_map:
1237 label = cube.attributes.get("model_name")
1238 color = model_colors_map.get(label)
1240 for cube_slice in cube.slices_over(ensemble_coord):
1241 # If ensemble data given plot control member with (control)
1242 # unless single forecast.
1243 if cube_slice.coord(ensemble_coord).points == [0]:
1244 iplt.plot(
1245 cube_slice,
1246 coord,
1247 color=color,
1248 marker="o",
1249 ls="-",
1250 lw=3,
1251 label=f"{label} (control)"
1252 if len(cube.coord(ensemble_coord).points) > 1
1253 else label,
1254 )
1255 # If ensemble data given plot perturbed members with (perturbed).
1256 else:
1257 iplt.plot(
1258 cube_slice,
1259 coord,
1260 color=color,
1261 ls="-",
1262 lw=1.5,
1263 alpha=0.75,
1264 label=f"{label} (member)",
1265 )
1267 # Get the current axis
1268 ax = plt.gca()
1270 # Special handling for pressure level data.
1271 if series_coordinate == "pressure":
1272 # Invert y-axis and set to log scale.
1273 ax.invert_yaxis()
1274 ax.set_yscale("log")
1276 # Define y-ticks and labels for pressure log axis.
1277 y_tick_labels = [
1278 "1000",
1279 "850",
1280 "700",
1281 "500",
1282 "300",
1283 "200",
1284 "100",
1285 ]
1286 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1288 # Set y-axis limits and ticks.
1289 ax.set_ylim(1100, 100)
1291 # Test if series_coordinate is model level data. The UM data uses
1292 # model_level_number and lfric uses full_levels as coordinate.
1293 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1294 # Define y-ticks and labels for vertical axis.
1295 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1296 y_tick_labels = [str(int(i)) for i in y_ticks]
1297 ax.set_ylim(min(y_ticks), max(y_ticks))
1299 ax.set_yticks(y_ticks)
1300 ax.set_yticklabels(y_tick_labels)
1302 # Set x-axis limits.
1303 ax.set_xlim(vmin, vmax)
1304 # Mark y=0 if present in plot.
1305 if vmin < 0.0 and vmax > 0.0:
1306 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1308 # Add some labels and tweak the style.
1309 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1310 ax.set_xlabel(
1311 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1312 )
1313 ax.set_title(title, fontsize=16)
1314 ax.ticklabel_format(axis="x")
1315 ax.tick_params(axis="y")
1316 ax.tick_params(axis="both", labelsize=12)
1318 # Add gridlines
1319 ax.grid(linestyle="--", color="grey", linewidth=1)
1320 # Ientify unique labels for legend
1321 handles = list(
1322 {
1323 label: handle
1324 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1325 }.values()
1326 )
1327 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
1329 # Save plot.
1330 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1331 logging.info("Saved line plot to %s", filename)
1332 plt.close(fig)
1335def _plot_and_save_scatter_plot(
1336 cube_x: iris.cube.Cube | iris.cube.CubeList,
1337 cube_y: iris.cube.Cube | iris.cube.CubeList,
1338 filename: str,
1339 title: str,
1340 one_to_one: bool,
1341 model_names: list[str] = None,
1342 **kwargs,
1343):
1344 """Plot and save a 2D scatter plot.
1346 Parameters
1347 ----------
1348 cube_x: Cube | CubeList
1349 1 dimensional Cube or CubeList of the data to plot on x-axis.
1350 cube_y: Cube | CubeList
1351 1 dimensional Cube or CubeList of the data to plot on y-axis.
1352 filename: str
1353 Filename of the plot to write.
1354 title: str
1355 Plot title.
1356 one_to_one: bool
1357 Whether a 1:1 line is plotted.
1358 """
1359 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1360 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1361 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1362 # over the pairs simultaneously.
1364 # Ensure cube_x and cube_y are iterable
1365 cube_x_iterable = iter_maybe(cube_x)
1366 cube_y_iterable = iter_maybe(cube_y)
1368 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1369 iplt.scatter(cube_x_iter, cube_y_iter)
1370 if one_to_one is True: 1370 ↛ 1383line 1370 didn't jump to line 1383 because the condition on line 1370 was always true
1371 plt.plot(
1372 [
1373 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1374 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1375 ],
1376 [
1377 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1378 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1379 ],
1380 "k",
1381 linestyle="--",
1382 )
1383 ax = plt.gca()
1385 # Add some labels and tweak the style.
1386 if model_names is None: 1386 ↛ 1387line 1386 didn't jump to line 1387 because the condition on line 1386 was never true
1387 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1388 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1389 else:
1390 # Add the model names, these should be order of base (x) and other (y).
1391 ax.set_xlabel(
1392 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1393 )
1394 ax.set_ylabel(
1395 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1396 )
1397 ax.set_title(title, fontsize=16)
1398 ax.ticklabel_format(axis="y", useOffset=False)
1399 ax.tick_params(axis="x", labelrotation=15)
1400 ax.tick_params(axis="both", labelsize=12)
1401 ax.autoscale()
1403 # Save plot.
1404 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1405 logging.info("Saved scatter plot to %s", filename)
1406 plt.close(fig)
1409def _plot_and_save_vector_plot(
1410 cube_u: iris.cube.Cube,
1411 cube_v: iris.cube.Cube,
1412 filename: str,
1413 title: str,
1414 method: Literal["contourf", "pcolormesh"],
1415 **kwargs,
1416):
1417 """Plot and save a 2D vector plot.
1419 Parameters
1420 ----------
1421 cube_u: Cube
1422 2 dimensional Cube of u component of the data.
1423 cube_v: Cube
1424 2 dimensional Cube of v component of the data.
1425 filename: str
1426 Filename of the plot to write.
1427 title: str
1428 Plot title.
1429 """
1430 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1431 # Create a cube containing the magnitude of the vector field.
1432 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1433 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1434 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1435 cube_vec_mag.rename(
1436 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1437 )
1439 # Specify the color bar
1440 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1442 # Setup plot map projection, extent and coastlines and borderlines.
1443 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1445 if method == "contourf":
1446 # Filled contour plot of the field.
1447 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1448 elif method == "pcolormesh":
1449 try:
1450 vmin = min(levels)
1451 vmax = max(levels)
1452 except TypeError:
1453 vmin, vmax = None, None
1454 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1455 # if levels are defined.
1456 if norm is not None:
1457 vmin = None
1458 vmax = None
1459 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1460 else:
1461 raise ValueError(f"Unknown plotting method: {method}")
1463 # Check to see if transect, and if so, adjust y axis.
1464 if is_transect(cube_vec_mag):
1465 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1466 axes.invert_yaxis()
1467 axes.set_yscale("log")
1468 axes.set_ylim(1100, 100)
1469 # If both model_level_number and level_height exists, iplt can construct
1470 # plot as a function of height above orography (NOT sea level).
1471 elif {"model_level_number", "level_height"}.issubset(
1472 {coord.name() for coord in cube_vec_mag.coords()}
1473 ):
1474 axes.set_yscale("log")
1476 axes.set_title(
1477 f"{title}\n"
1478 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1479 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1480 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1481 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1482 fontsize=16,
1483 )
1485 else:
1486 # Add title.
1487 axes.set_title(title, fontsize=16)
1489 # Add watermark with min/max/mean. Currently not user togglable.
1490 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1491 axes.annotate(
1492 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}",
1493 xy=(0.05, -0.05),
1494 xycoords="axes fraction",
1495 xytext=(-5, 5),
1496 textcoords="offset points",
1497 ha="right",
1498 va="bottom",
1499 size=11,
1500 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1501 )
1503 # Add colour bar.
1504 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1505 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1506 # add ticks and tick_labels for every levels if less than 20 levels exist
1507 if levels is not None and len(levels) < 20:
1508 cbar.set_ticks(levels)
1509 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1511 # 30 barbs along the longest axis of the plot, or a barb per point for data
1512 # with less than 30 points.
1513 step = max(max(cube_u.shape) // 30, 1)
1514 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1516 # Save plot.
1517 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1518 logging.info("Saved vector plot to %s", filename)
1519 plt.close(fig)
1522def _plot_and_save_histogram_series(
1523 cubes: iris.cube.Cube | iris.cube.CubeList,
1524 filename: str,
1525 title: str,
1526 vmin: float,
1527 vmax: float,
1528 **kwargs,
1529):
1530 """Plot and save a histogram series.
1532 Parameters
1533 ----------
1534 cubes: Cube or CubeList
1535 2 dimensional Cube or CubeList of the data to plot as histogram.
1536 filename: str
1537 Filename of the plot to write.
1538 title: str
1539 Plot title.
1540 vmin: float
1541 minimum for colorbar
1542 vmax: float
1543 maximum for colorbar
1544 """
1545 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1546 ax = plt.gca()
1548 model_colors_map = get_model_colors_map(cubes)
1550 # Set default that histograms will produce probability density function
1551 # at each bin (integral over range sums to 1).
1552 density = True
1554 for cube in iter_maybe(cubes):
1555 # Easier to check title (where var name originates)
1556 # than seeing if long names exist etc.
1557 # Exception case, where distribution better fits log scales/bins.
1558 if "surface_microphysical" in title:
1559 if "amount" in title:
1560 # Compute histogram following Klingaman et al. (2017): ASoP
1561 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1562 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1563 density = False
1564 else:
1565 bins = 10.0 ** (
1566 np.arange(-10, 27, 1) / 10.0
1567 ) # Suggestion from RMED toolbox.
1568 bins = np.insert(bins, 0, 0)
1569 ax.set_yscale("log")
1570 vmin = bins[1]
1571 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1572 ax.set_xscale("log")
1573 elif "lightning" in title:
1574 bins = [0, 1, 2, 3, 4, 5]
1575 else:
1576 bins = np.linspace(vmin, vmax, 51)
1577 logging.debug(
1578 "Plotting histogram with %s bins %s - %s.",
1579 np.size(bins),
1580 np.min(bins),
1581 np.max(bins),
1582 )
1584 # Reshape cube data into a single array to allow for a single histogram.
1585 # Otherwise we plot xdim histograms stacked.
1586 cube_data_1d = (cube.data).flatten()
1588 label = None
1589 color = "black"
1590 if model_colors_map:
1591 label = cube.attributes.get("model_name")
1592 color = model_colors_map[label]
1593 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1595 # Compute area under curve.
1596 if "surface_microphysical" in title and "amount" in title:
1597 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1598 x = x * bin_mean / x.sum()
1599 x = x[1:]
1600 y = y[1:]
1602 ax.plot(
1603 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1604 )
1606 # Add some labels and tweak the style.
1607 ax.set_title(title, fontsize=16)
1608 ax.set_xlabel(
1609 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1610 )
1611 ax.set_ylabel("Normalised probability density", fontsize=14)
1612 if "surface_microphysical" in title and "amount" in title:
1613 ax.set_ylabel(
1614 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1615 )
1616 ax.set_xlim(vmin, vmax)
1617 ax.tick_params(axis="both", labelsize=12)
1619 # Overlay grid-lines onto histogram plot.
1620 ax.grid(linestyle="--", color="grey", linewidth=1)
1621 if model_colors_map:
1622 ax.legend(loc="best", ncol=1, frameon=False, fontsize=16)
1624 # Save plot.
1625 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1626 logging.info("Saved histogram plot to %s", filename)
1627 plt.close(fig)
1630def _plot_and_save_postage_stamp_histogram_series(
1631 cube: iris.cube.Cube,
1632 filename: str,
1633 title: str,
1634 stamp_coordinate: str,
1635 vmin: float,
1636 vmax: float,
1637 **kwargs,
1638):
1639 """Plot and save postage (ensemble members) stamps for a histogram series.
1641 Parameters
1642 ----------
1643 cube: Cube
1644 2 dimensional Cube of the data to plot as histogram.
1645 filename: str
1646 Filename of the plot to write.
1647 title: str
1648 Plot title.
1649 stamp_coordinate: str
1650 Coordinate that becomes different plots.
1651 vmin: float
1652 minimum for pdf x-axis
1653 vmax: float
1654 maximum for pdf x-axis
1655 """
1656 # Use the smallest square grid that will fit the members.
1657 nmember = len(cube.coord(stamp_coordinate).points)
1658 grid_rows = int(math.sqrt(nmember))
1659 grid_size = math.ceil(nmember / grid_rows)
1661 fig = plt.figure(
1662 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1663 )
1664 # Make a subplot for each member.
1665 for member, subplot in zip(
1666 cube.slices_over(stamp_coordinate),
1667 range(1, grid_size * grid_rows + 1),
1668 strict=False,
1669 ):
1670 # Implicit interface is much easier here, due to needing to have the
1671 # cartopy GeoAxes generated.
1672 plt.subplot(grid_rows, grid_size, subplot)
1673 # Reshape cube data into a single array to allow for a single histogram.
1674 # Otherwise we plot xdim histograms stacked.
1675 member_data_1d = (member.data).flatten()
1676 plt.hist(member_data_1d, density=True, stacked=True)
1677 axes = plt.gca()
1678 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1679 axes.set_title(f"{mtitle}")
1680 axes.set_xlim(vmin, vmax)
1682 # Overall figure title.
1683 fig.suptitle(title, fontsize=16)
1685 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1686 logging.info("Saved histogram postage stamp plot to %s", filename)
1687 plt.close(fig)
1690def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1691 cube: iris.cube.Cube,
1692 filename: str,
1693 title: str,
1694 stamp_coordinate: str,
1695 vmin: float,
1696 vmax: float,
1697 **kwargs,
1698):
1699 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1700 ax.set_title(title, fontsize=16)
1701 ax.set_xlim(vmin, vmax)
1702 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1703 ax.set_ylabel("normalised probability density", fontsize=14)
1704 # Loop over all slices along the stamp_coordinate
1705 for member in cube.slices_over(stamp_coordinate):
1706 # Flatten the member data to 1D
1707 member_data_1d = member.data.flatten()
1708 # Plot the histogram using plt.hist
1709 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1710 plt.hist(
1711 member_data_1d,
1712 density=True,
1713 stacked=True,
1714 label=f"{mtitle}",
1715 )
1717 # Add a legend
1718 ax.legend(fontsize=16)
1720 # Save the figure to a file
1721 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1722 logging.info("Saved histogram postage stamp plot to %s", filename)
1724 # Close the figure
1725 plt.close(fig)
1728def _plot_and_save_scattermap_plot(
1729 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1730):
1731 """Plot and save a geographical scatter plot.
1733 Parameters
1734 ----------
1735 cube: Cube
1736 1 dimensional Cube of the data points with auxiliary latitude and
1737 longitude coordinates,
1738 filename: str
1739 Filename of the plot to write.
1740 title: str
1741 Plot title.
1742 projection: str
1743 Mapping projection to be used by cartopy.
1744 """
1745 # Setup plot details, size, resolution, etc.
1746 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1747 if projection is not None:
1748 # Apart from the default, the only projection we currently support is
1749 # a stereographic projection over the North Pole.
1750 if projection == "NP_Stereo":
1751 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1752 else:
1753 raise ValueError(f"Unknown projection: {projection}")
1754 else:
1755 axes = plt.axes(projection=ccrs.PlateCarree())
1757 # Scatter plot of the field. The marker size is chosen to give
1758 # symbols that decrease in size as the number of observations
1759 # increases, although the fraction of the figure covered by
1760 # symbols increases roughly as N^(1/2), disregarding overlaps,
1761 # and has been selected for the default figure size of (10, 10).
1762 # Should this be changed, the marker size should be adjusted in
1763 # proportion to the area of the figure.
1764 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1765 klon = None
1766 klat = None
1767 for kc in range(len(cube.aux_coords)):
1768 if cube.aux_coords[kc].standard_name == "latitude":
1769 klat = kc
1770 elif cube.aux_coords[kc].standard_name == "longitude":
1771 klon = kc
1772 scatter_map = iplt.scatter(
1773 cube.aux_coords[klon],
1774 cube.aux_coords[klat],
1775 c=cube.data[:],
1776 s=mrk_size,
1777 cmap="jet",
1778 edgecolors="k",
1779 )
1781 # Add coastlines and borderlines.
1782 try:
1783 axes.coastlines(resolution="10m")
1784 axes.add_feature(cfeature.BORDERS)
1785 except AttributeError:
1786 pass
1788 # Add title.
1789 axes.set_title(title, fontsize=16)
1791 # Add colour bar.
1792 cbar = fig.colorbar(scatter_map)
1793 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1795 # Save plot.
1796 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1797 logging.info("Saved geographical scatter plot to %s", filename)
1798 plt.close(fig)
1801def _spatial_plot(
1802 method: Literal["contourf", "pcolormesh"],
1803 cube: iris.cube.Cube,
1804 filename: str | None,
1805 sequence_coordinate: str,
1806 stamp_coordinate: str,
1807 overlay_cube: iris.cube.Cube | None = None,
1808 contour_cube: iris.cube.Cube | None = None,
1809 **kwargs,
1810):
1811 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1813 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1814 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1815 is present then postage stamp plots will be produced.
1817 If an overlay_cube and/or contour_cube are specified, multiple variables can
1818 be overplotted on the same figure.
1820 Parameters
1821 ----------
1822 method: "contourf" | "pcolormesh"
1823 The plotting method to use.
1824 cube: Cube
1825 Iris cube of the data to plot. It should have two spatial dimensions,
1826 such as lat and lon, and may also have a another two dimension to be
1827 plotted sequentially and/or as postage stamp plots.
1828 filename: str | None
1829 Name of the plot to write, used as a prefix for plot sequences. If None
1830 uses the recipe name.
1831 sequence_coordinate: str
1832 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1833 This coordinate must exist in the cube.
1834 stamp_coordinate: str
1835 Coordinate about which to plot postage stamp plots. Defaults to
1836 ``"realization"``.
1837 overlay_cube: Cube | None, optional
1838 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1839 contour_cube: Cube | None, optional
1840 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1842 Raises
1843 ------
1844 ValueError
1845 If the cube doesn't have the right dimensions.
1846 TypeError
1847 If the cube isn't a single cube.
1848 """
1849 recipe_title = get_recipe_metadata().get("title", "Untitled")
1851 # Ensure we've got a single cube.
1852 cube = check_single_cube(cube)
1854 # Check if there is a valid stamp coordinate in cube dimensions.
1855 if stamp_coordinate == "realization": 1855 ↛ 1860line 1855 didn't jump to line 1860 because the condition on line 1855 was always true
1856 stamp_coordinate = check_stamp_coordinate(cube)
1858 # Make postage stamp plots if stamp_coordinate exists and has more than a
1859 # single point.
1860 plotting_func = _plot_and_save_spatial_plot
1861 try:
1862 if cube.coord(stamp_coordinate).shape[0] > 1:
1863 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1864 except iris.exceptions.CoordinateNotFoundError:
1865 pass
1867 # Produce a geographical scatter plot if the data have a
1868 # dimension called observation or model_obs_error
1869 if any( 1869 ↛ 1873line 1869 didn't jump to line 1873 because the condition on line 1869 was never true
1870 crd.var_name == "station" or crd.var_name == "model_obs_error"
1871 for crd in cube.coords()
1872 ):
1873 plotting_func = _plot_and_save_scattermap_plot
1875 # Must have a sequence coordinate.
1876 try:
1877 cube.coord(sequence_coordinate)
1878 except iris.exceptions.CoordinateNotFoundError as err:
1879 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1881 # Create a plot for each value of the sequence coordinate.
1882 plot_index = []
1883 nplot = np.size(cube.coord(sequence_coordinate).points)
1885 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1886 # Set plot titles and filename
1887 seq_coord = cube_slice.coord(sequence_coordinate)
1888 plot_title, plot_filename = _set_title_and_filename(
1889 seq_coord, nplot, recipe_title, filename
1890 )
1892 # Extract sequence slice for overlay_cube and contour_cube if required.
1893 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1894 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1896 # Do the actual plotting.
1897 plotting_func(
1898 cube_slice,
1899 filename=plot_filename,
1900 stamp_coordinate=stamp_coordinate,
1901 title=plot_title,
1902 method=method,
1903 overlay_cube=overlay_slice,
1904 contour_cube=contour_slice,
1905 **kwargs,
1906 )
1907 plot_index.append(plot_filename)
1909 # Add list of plots to plot metadata.
1910 complete_plot_index = _append_to_plot_index(plot_index)
1912 # Make a page to display the plots.
1913 _make_plot_html_page(complete_plot_index)
1916####################
1917# Public functions #
1918####################
1921def spatial_contour_plot(
1922 cube: iris.cube.Cube,
1923 filename: str = None,
1924 sequence_coordinate: str = "time",
1925 stamp_coordinate: str = "realization",
1926 **kwargs,
1927) -> iris.cube.Cube:
1928 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1930 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1931 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1932 is present then postage stamp plots will be produced.
1934 Parameters
1935 ----------
1936 cube: Cube
1937 Iris cube of the data to plot. It should have two spatial dimensions,
1938 such as lat and lon, and may also have a another two dimension to be
1939 plotted sequentially and/or as postage stamp plots.
1940 filename: str, optional
1941 Name of the plot to write, used as a prefix for plot sequences. Defaults
1942 to the recipe name.
1943 sequence_coordinate: str, optional
1944 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1945 This coordinate must exist in the cube.
1946 stamp_coordinate: str, optional
1947 Coordinate about which to plot postage stamp plots. Defaults to
1948 ``"realization"``.
1950 Returns
1951 -------
1952 Cube
1953 The original cube (so further operations can be applied).
1955 Raises
1956 ------
1957 ValueError
1958 If the cube doesn't have the right dimensions.
1959 TypeError
1960 If the cube isn't a single cube.
1961 """
1962 _spatial_plot(
1963 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1964 )
1965 return cube
1968def spatial_pcolormesh_plot(
1969 cube: iris.cube.Cube,
1970 filename: str = None,
1971 sequence_coordinate: str = "time",
1972 stamp_coordinate: str = "realization",
1973 **kwargs,
1974) -> iris.cube.Cube:
1975 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1977 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1978 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1979 is present then postage stamp plots will be produced.
1981 This function is significantly faster than ``spatial_contour_plot``,
1982 especially at high resolutions, and should be preferred unless contiguous
1983 contour areas are important.
1985 Parameters
1986 ----------
1987 cube: Cube
1988 Iris cube of the data to plot. It should have two spatial dimensions,
1989 such as lat and lon, and may also have a another two dimension to be
1990 plotted sequentially and/or as postage stamp plots.
1991 filename: str, optional
1992 Name of the plot to write, used as a prefix for plot sequences. Defaults
1993 to the recipe name.
1994 sequence_coordinate: str, optional
1995 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1996 This coordinate must exist in the cube.
1997 stamp_coordinate: str, optional
1998 Coordinate about which to plot postage stamp plots. Defaults to
1999 ``"realization"``.
2001 Returns
2002 -------
2003 Cube
2004 The original cube (so further operations can be applied).
2006 Raises
2007 ------
2008 ValueError
2009 If the cube doesn't have the right dimensions.
2010 TypeError
2011 If the cube isn't a single cube.
2012 """
2013 _spatial_plot(
2014 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
2015 )
2016 return cube
2019def spatial_multi_pcolormesh_plot(
2020 cube: iris.cube.Cube,
2021 overlay_cube: iris.cube.Cube | None = None,
2022 contour_cube: iris.cube.Cube | None = None,
2023 filename: str = None,
2024 sequence_coordinate: str = "time",
2025 stamp_coordinate: str = "realization",
2026 **kwargs,
2027) -> iris.cube.Cube:
2028 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2030 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2031 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2032 is present then postage stamp plots will be produced.
2034 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2036 If specified, contours of a contour_cube can be overplotted on top of those.
2038 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2040 This function is significantly faster than ``spatial_contour_plot``,
2041 especially at high resolutions, and should be preferred unless contiguous
2042 contour areas are important.
2044 Parameters
2045 ----------
2046 cube: Cube
2047 Iris cube of the data to plot. It should have two spatial dimensions,
2048 such as lat and lon, and may also have a another two dimension to be
2049 plotted sequentially and/or as postage stamp plots.
2050 overlay_cube: Cube, optional
2051 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2052 such as lat and lon, and may also have a another two dimension to be
2053 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.
2054 If not provided, output plot generated without overlay cube.
2055 contour_cube: Cube, optional
2056 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,
2057 such as lat and lon, and may also have a another two dimension to be
2058 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
2059 filename: str, optional
2060 Name of the plot to write, used as a prefix for plot sequences. Defaults
2061 to the recipe name.
2062 sequence_coordinate: str, optional
2063 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2064 This coordinate must exist in the cube.
2065 stamp_coordinate: str, optional
2066 Coordinate about which to plot postage stamp plots. Defaults to
2067 ``"realization"``.
2069 Returns
2070 -------
2071 Cube
2072 The original cube (so further operations can be applied).
2074 Raises
2075 ------
2076 ValueError
2077 If the cube doesn't have the right dimensions.
2078 TypeError
2079 If the cube isn't a single cube.
2080 """
2081 _spatial_plot(
2082 "pcolormesh",
2083 cube,
2084 filename,
2085 sequence_coordinate,
2086 stamp_coordinate,
2087 overlay_cube=overlay_cube,
2088 contour_cube=contour_cube,
2089 )
2090 return cube, overlay_cube, contour_cube
2093# TODO: Expand function to handle ensemble data.
2094# line_coordinate: str, optional
2095# Coordinate about which to plot multiple lines. Defaults to
2096# ``"realization"``.
2097def plot_line_series(
2098 cube: iris.cube.Cube | iris.cube.CubeList,
2099 filename: str = None,
2100 series_coordinate: str = "time",
2101 sequence_coordinate: str = "time",
2102 # add the following for ensembles
2103 stamp_coordinate: str = "realization",
2104 single_plot: bool = False,
2105 **kwargs,
2106) -> iris.cube.Cube | iris.cube.CubeList:
2107 """Plot a line plot for the specified coordinate.
2109 The Cube or CubeList must be 1D.
2111 Parameters
2112 ----------
2113 iris.cube | iris.cube.CubeList
2114 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2115 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2116 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2117 filename: str, optional
2118 Name of the plot to write, used as a prefix for plot sequences. Defaults
2119 to the recipe name.
2120 series_coordinate: str, optional
2121 Coordinate about which to make a series. Defaults to ``"time"``. This
2122 coordinate must exist in the cube.
2124 Returns
2125 -------
2126 iris.cube.Cube | iris.cube.CubeList
2127 The original Cube or CubeList (so further operations can be applied).
2128 plotted data.
2130 Raises
2131 ------
2132 ValueError
2133 If the cubes don't have the right dimensions.
2134 TypeError
2135 If the cube isn't a Cube or CubeList.
2136 """
2137 # Ensure we have a name for the plot file.
2138 recipe_title = get_recipe_metadata().get("title", "Untitled")
2140 num_models = get_num_models(cube)
2142 validate_cube_shape(cube, num_models)
2144 # Iterate over all cubes and extract coordinate to plot.
2145 cubes = iter_maybe(cube)
2146 coords = []
2147 for cube in cubes:
2148 try:
2149 coords.append(cube.coord(series_coordinate))
2150 except iris.exceptions.CoordinateNotFoundError as err:
2151 raise ValueError(
2152 f"Cube must have a {series_coordinate} coordinate."
2153 ) from err
2154 if cube.coords("realization"): 2154 ↛ 2158line 2154 didn't jump to line 2158 because the condition on line 2154 was always true
2155 if cube.ndim > 3: 2155 ↛ 2156line 2155 didn't jump to line 2156 because the condition on line 2155 was never true
2156 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2157 else:
2158 raise ValueError("Cube must have a realization coordinate.")
2160 plot_index = []
2162 # Check if this is a spectral plot by looking for spectral coordinates
2163 is_spectral_plot = series_coordinate in [
2164 "frequency",
2165 "physical_wavenumber",
2166 "wavelength",
2167 ]
2169 if is_spectral_plot:
2170 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2171 # coordinate frequency/wavenumber.
2172 # If several power spectra are plotted with time as sequence_coordinate for the
2173 # time slider option.
2175 # Internal plotting function.
2176 plotting_func = _plot_and_save_line_power_spectrum_series
2178 for cube in cubes:
2179 try:
2180 cube.coord(sequence_coordinate)
2181 except iris.exceptions.CoordinateNotFoundError as err:
2182 raise ValueError(
2183 f"Cube must have a {sequence_coordinate} coordinate."
2184 ) from err
2186 if num_models == 1: 2186 ↛ 2200line 2186 didn't jump to line 2200 because the condition on line 2186 was always true
2187 # check for ensembles
2188 if ( 2188 ↛ 2192line 2188 didn't jump to line 2192 because the condition on line 2188 was never true
2189 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2190 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2191 ):
2192 if single_plot:
2193 # Plot spectra, mean and ensemble spread on 1 plot
2194 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2195 else:
2196 # Plot postage stamps
2197 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2198 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2199 else:
2200 all_points = sorted(
2201 set(
2202 itertools.chain.from_iterable(
2203 cb.coord(sequence_coordinate).points for cb in cubes
2204 )
2205 )
2206 )
2207 all_slices = list(
2208 itertools.chain.from_iterable(
2209 cb.slices_over(sequence_coordinate) for cb in cubes
2210 )
2211 )
2212 # Matched slices (matched by seq coord point; it may happen that
2213 # evaluated models do not cover the same seq coord range, hence matching
2214 # necessary)
2215 cube_iterables = [
2216 iris.cube.CubeList(
2217 s
2218 for s in all_slices
2219 if s.coord(sequence_coordinate).points[0] == point
2220 )
2221 for point in all_points
2222 ]
2224 nplot = np.size(cube.coord(sequence_coordinate).points)
2226 # Create a plot for each value of the sequence coordinate. Allowing for
2227 # multiple cubes in a CubeList to be plotted in the same plot for similar
2228 # sequence values. Passing a CubeList into the internal plotting function
2229 # for similar values of the sequence coordinate. cube_slice can be an
2230 # iris.cube.Cube or an iris.cube.CubeList.
2232 for cube_slice in cube_iterables:
2233 # Normalize cube_slice to a list of cubes
2234 if isinstance(cube_slice, iris.cube.CubeList): 2234 ↛ 2235line 2234 didn't jump to line 2235 because the condition on line 2234 was never true
2235 cubes = list(cube_slice)
2236 elif isinstance(cube_slice, iris.cube.Cube): 2236 ↛ 2239line 2236 didn't jump to line 2239 because the condition on line 2236 was always true
2237 cubes = [cube_slice]
2238 else:
2239 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2241 # Use sequence value so multiple sequences can merge.
2242 seq_coord = cube_slice[0].coord(sequence_coordinate)
2243 plot_title, plot_filename = _set_title_and_filename(
2244 seq_coord, nplot, recipe_title, filename
2245 )
2247 # Format the coordinate value in a unit appropriate way.
2248 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2250 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2251 if nplot == 1 and seq_coord.has_bounds: 2251 ↛ 2256line 2251 didn't jump to line 2256 because the condition on line 2251 was always true
2252 if np.size(seq_coord.bounds) > 1: 2252 ↛ 2253line 2252 didn't jump to line 2253 because the condition on line 2252 was never true
2253 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2255 # Do the actual plotting.
2256 plotting_func(
2257 cube_slice,
2258 coords,
2259 stamp_coordinate,
2260 plot_filename,
2261 title,
2262 series_coordinate,
2263 )
2265 plot_index.append(plot_filename)
2266 else:
2267 # Format the title and filename using plotted series coordinate
2268 nplot = 1
2269 seq_coord = coords[0]
2270 plot_title, plot_filename = _set_title_and_filename(
2271 seq_coord, nplot, recipe_title, filename
2272 )
2273 # Do the actual plotting for all other series coordinate options.
2274 _plot_and_save_line_series(
2275 cubes, coords, stamp_coordinate, plot_filename, plot_title
2276 )
2278 plot_index.append(plot_filename)
2279 # Do the actual plotting.
2280 _plot_and_save_line_series(
2281 cubes, coords, plot_filename, plot_title, ensemble_coord="realization"
2282 )
2284 # append plot to list of plots
2285 complete_plot_index = _append_to_plot_index(plot_index)
2287 # Make a page to display the plots.
2288 _make_plot_html_page(complete_plot_index)
2290 return cube
2293def plot_dfss_line_series_sequence(
2294 cube: iris.cube.Cube | iris.cube.CubeList,
2295 filename: str = None,
2296 variable: str = None,
2297 series_coordinate: str = None,
2298 sequence_coordinate: str = None,
2299 **kwargs,
2300) -> iris.cube.Cube | iris.cube.CubeList:
2301 """Plot a line plot."""
2302 cube_copy = cube
2303 if type(cube) is iris.cube.CubeList:
2304 if not variable:
2305 logging.Warning(
2306 "CubeList given, but variable not specified. Defaulting to first cube."
2307 )
2308 cube = cube[0]
2309 else:
2310 cube = cube.extract(variable)
2312 recipe_title = get_recipe_metadata().get("title", "Untitled")
2314 num_models = get_num_models(cube)
2316 validate_cube_shape(cube, num_models)
2318 # Iterate over all cubes and extract coordinate to plot.
2319 cubes = iter_maybe(cube)
2321 coords = []
2322 for cube in cubes:
2323 try:
2324 coords.append(cube.coord(series_coordinate))
2325 except iris.exceptions.CoordinateNotFoundError as err:
2326 raise ValueError(
2327 f"Cube must have a {series_coordinate} coordinate."
2328 ) from err
2329 if cube.ndim > 2 or not cube.coords(sequence_coordinate):
2330 raise ValueError(
2331 f"Cube must be 1D or 2D with a {sequence_coordinate} coordinate."
2332 )
2334 # Format the title and filename using plotted series coordinate
2335 nplot = 1
2336 seq_coord = coords[0]
2337 plot_title, plot_filename = _set_title_and_filename(
2338 seq_coord, nplot, recipe_title, filename
2339 )
2340 # Do the actual plotting
2342 for i, cubes in enumerate(cube.slices_over(sequence_coordinate)):
2343 if cubes.coord(sequence_coordinate).units == "unknown":
2344 sequence_point = cubes.coord(sequence_coordinate).points[0]
2345 else:
2346 sequence_point = cubes.coord(sequence_coordinate).units.title(
2347 cubes.coord(sequence_coordinate).points[0]
2348 )
2350 if cube.attributes.locals["method"] == "centile":
2351 method = cube.attributes.locals["method"]
2352 centile = cube.attributes.locals["centile"]
2353 plot_filename_with_sequence_coord = f"{cube.name()}_{sequence_coordinate}_point_{str(i)}_{method}_{centile}_{plot_filename}"
2354 plot_title_with_time = f"{cubes.name()} vs {series_coordinate} ({sequence_coordinate}: {sequence_point}) \n method: Centile | centile: {centile}"
2355 elif cube.attributes.locals["method"] == "threshold":
2356 method = cube.attributes.locals["method"]
2357 threshold = cube.attributes.locals["threshold"]
2358 plot_filename_with_sequence_coord = f"{cube.name()}_{sequence_coordinate}_point_{str(i)}_{method}_{threshold}_{plot_filename}"
2359 plot_title_with_time = f"{cubes.name()} vs {series_coordinate} ({sequence_coordinate}: {sequence_point}) \n method: Threshold | threshold: {threshold}"
2361 cubes_in = iter_maybe(cubes)
2362 _plot_and_save_line_series(
2363 cubes_in,
2364 coords,
2365 plot_filename_with_sequence_coord,
2366 plot_title_with_time,
2367 sequence_coord=sequence_coordinate,
2368 )
2370 # Add list of plots to plot metadata.
2371 plot_index = _append_to_plot_index([plot_filename_with_sequence_coord])
2373 # Make a page to display the plots.
2374 _make_plot_html_page(plot_index)
2376 return cube_copy
2379def plot_vertical_line_series(
2380 cubes: iris.cube.Cube | iris.cube.CubeList,
2381 filename: str = None,
2382 series_coordinate: str = "model_level_number",
2383 sequence_coordinate: str = "time",
2384 # line_coordinate: str = "realization",
2385 **kwargs,
2386) -> iris.cube.Cube | iris.cube.CubeList:
2387 """Plot a line plot against a type of vertical coordinate.
2389 The Cube or CubeList must be 1D.
2391 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2392 then a sequence of plots will be produced.
2394 Parameters
2395 ----------
2396 iris.cube | iris.cube.CubeList
2397 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2398 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2399 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2400 filename: str, optional
2401 Name of the plot to write, used as a prefix for plot sequences. Defaults
2402 to the recipe name.
2403 series_coordinate: str, optional
2404 Coordinate to plot on the y-axis. Can be ``pressure`` or
2405 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2406 for LFRic. Defaults to ``model_level_number``.
2407 This coordinate must exist in the cube.
2408 sequence_coordinate: str, optional
2409 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2410 This coordinate must exist in the cube.
2412 Returns
2413 -------
2414 iris.cube.Cube | iris.cube.CubeList
2415 The original Cube or CubeList (so further operations can be applied).
2416 Plotted data.
2418 Raises
2419 ------
2420 ValueError
2421 If the cubes doesn't have the right dimensions.
2422 TypeError
2423 If the cube isn't a Cube or CubeList.
2424 """
2425 # Ensure we have a name for the plot file.
2426 recipe_title = get_recipe_metadata().get("title", "Untitled")
2428 cubes = iter_maybe(cubes)
2429 # Initialise empty list to hold all data from all cubes in a CubeList
2430 all_data = []
2432 # Store min/max ranges for x range.
2433 x_levels = []
2435 num_models = get_num_models(cubes)
2437 validate_cube_shape(cubes, num_models)
2439 # Iterate over all cubes in cube or CubeList and plot.
2440 coords = []
2441 for cube in cubes: 2441 ↛ 2466line 2441 didn't jump to line 2466 because the loop on line 2441 didn't complete
2442 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2443 try:
2444 coords.append(cube.coord(series_coordinate))
2445 except iris.exceptions.CoordinateNotFoundError as err:
2446 raise ValueError(
2447 f"Cube must have a {series_coordinate} coordinate."
2448 ) from err
2450 try:
2451 if cube.ndim > 1 or not cube.coords("realization"):
2452 cube.coord(sequence_coordinate)
2453 except iris.exceptions.CoordinateNotFoundError as err:
2454 raise ValueError(
2455 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2456 ) from err
2458 # Get minimum and maximum from levels information.
2459 _, levels, _ = colorbar_map_levels(cube, axis="x")
2460 if levels is not None:
2461 x_levels.append(min(levels))
2462 x_levels.append(max(levels))
2463 else:
2464 all_data.append(cube.data)
2466 if len(x_levels) == 0:
2467 # Combine all data into a single NumPy array
2468 combined_data = np.concatenate(all_data)
2470 # Set the lower and upper limit for the x-axis to ensure all plots have
2471 # same range. This needs to read the whole cube over the range of the
2472 # sequence and if applicable postage stamp coordinate.
2473 vmin = np.floor(combined_data.min())
2474 vmax = np.ceil(combined_data.max())
2475 else:
2476 vmin = min(x_levels)
2477 vmax = max(x_levels)
2479 # Matching the slices (matching by seq coord point; it may happen that
2480 # evaluated models do not cover the same seq coord range, hence matching
2481 # necessary)
2482 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2484 # Create a plot for each value of the sequence coordinate.
2485 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2486 # similar sequence values. Passing a CubeList into the internal plotting function
2487 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2488 # or an iris.cube.CubeList.
2489 plot_index = []
2490 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2491 for cubes_slice in cube_iterables:
2492 # Format the coordinate value in a unit appropriate way.
2493 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2494 plot_title, plot_filename = _set_title_and_filename(
2495 seq_coord, nplot, recipe_title, filename
2496 )
2498 # Do the actual plotting.
2499 _plot_and_save_vertical_line_series(
2500 cubes_slice,
2501 coords,
2502 "realization",
2503 plot_filename,
2504 series_coordinate,
2505 title=plot_title,
2506 vmin=vmin,
2507 vmax=vmax,
2508 )
2509 plot_index.append(plot_filename)
2511 # Add list of plots to plot metadata.
2512 complete_plot_index = _append_to_plot_index(plot_index)
2514 # Make a page to display the plots.
2515 _make_plot_html_page(complete_plot_index)
2517 return cubes
2520def qq_plot(
2521 cubes: iris.cube.CubeList,
2522 coordinates: list[str],
2523 percentiles: list[float],
2524 model_names: list[str],
2525 filename: str = None,
2526 one_to_one: bool = True,
2527 **kwargs,
2528) -> iris.cube.CubeList:
2529 """Plot a Quantile-Quantile plot between two models for common time points.
2531 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2532 collapsed within the operator over all specified coordinates such as
2533 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2534 ensemble members to ensure a 1D cube (array).
2536 Parameters
2537 ----------
2538 cubes: iris.cube.CubeList
2539 Two cubes of the same variable with different models.
2540 coordinate: list[str]
2541 The list of coordinates to collapse over. This list should be
2542 every coordinate within the cube to result in a 1D cube around
2543 the percentile coordinate.
2544 percent: list[float]
2545 A list of percentiles to appear in the plot.
2546 model_names: list[str]
2547 A list of model names to appear on the axis of the plot.
2548 filename: str, optional
2549 Filename of the plot to write.
2550 one_to_one: bool, optional
2551 If True a 1:1 line is plotted; if False it is not. Default is True.
2553 Raises
2554 ------
2555 ValueError
2556 When the cubes are not compatible.
2558 Notes
2559 -----
2560 The quantile-quantile plot is a variant on the scatter plot representing
2561 two datasets by their quantiles (percentiles) for common time points.
2562 This plot does not use a theoretical distribution to compare against, but
2563 compares percentiles of two datasets. This plot does
2564 not use all raw data points, but plots the selected percentiles (quantiles) of
2565 each variable instead for the two datasets, thereby normalising the data for a
2566 direct comparison between the selected percentiles of the two dataset distributions.
2568 Quantile-quantile plots are valuable for comparing against
2569 observations and other models. Identical percentiles between the variables
2570 will lie on the one-to-one line implying the values correspond well to each
2571 other. Where there is a deviation from the one-to-one line a range of
2572 possibilities exist depending on how and where the data is shifted (e.g.,
2573 Wilks 2011 [Wilks2011]_).
2575 For distributions above the one-to-one line the distribution is left-skewed;
2576 below is right-skewed. A distinct break implies a bimodal distribution, and
2577 closer values/values further apart at the tails imply poor representation of
2578 the extremes.
2580 References
2581 ----------
2582 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2583 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2584 """
2585 # Check cubes using same functionality as the difference operator.
2586 if len(cubes) != 2:
2587 raise ValueError("cubes should contain exactly 2 cubes.")
2588 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2589 other: Cube = cubes.extract_cube(
2590 iris.Constraint(
2591 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2592 )
2593 )
2595 # Get spatial coord names.
2596 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2597 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2599 # Ensure cubes to compare are on common differencing grid.
2600 # This is triggered if either
2601 # i) latitude and longitude shapes are not the same. Note grid points
2602 # are not compared directly as these can differ through rounding
2603 # errors.
2604 # ii) or variables are known to often sit on different grid staggering
2605 # in different models (e.g. cell center vs cell edge), as is the case
2606 # for UM and LFRic comparisons.
2607 # In future greater choice of regridding method might be applied depending
2608 # on variable type. Linear regridding can in general be appropriate for smooth
2609 # variables. Care should be taken with interpretation of differences
2610 # given this dependency on regridding.
2611 if (
2612 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2613 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2614 ) or (
2615 base.long_name
2616 in [
2617 "eastward_wind_at_10m",
2618 "northward_wind_at_10m",
2619 "northward_wind_at_cell_centres",
2620 "eastward_wind_at_cell_centres",
2621 "zonal_wind_at_pressure_levels",
2622 "meridional_wind_at_pressure_levels",
2623 "potential_vorticity_at_pressure_levels",
2624 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2625 ]
2626 ):
2627 logging.debug(
2628 "Linear regridding base cube to other grid to compute differences"
2629 )
2630 base = regrid_onto_cube(base, other, method="Linear")
2632 # Extract just common time points.
2633 base, other = _extract_common_time_points(base, other)
2635 # Equalise attributes so we can merge.
2636 fully_equalise_attributes([base, other])
2637 logging.debug("Base: %s\nOther: %s", base, other)
2639 # Collapse cubes.
2640 base = collapse(
2641 base,
2642 coordinate=coordinates,
2643 method="PERCENTILE",
2644 additional_percent=percentiles,
2645 )
2646 other = collapse(
2647 other,
2648 coordinate=coordinates,
2649 method="PERCENTILE",
2650 additional_percent=percentiles,
2651 )
2653 # Ensure we have a name for the plot file.
2654 recipe_title = get_recipe_metadata().get("title", "Untitled")
2655 title = f"{recipe_title}"
2657 if filename is None:
2658 filename = slugify(recipe_title)
2660 # Add file extension.
2661 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2663 # Do the actual plotting on a scatter plot
2664 _plot_and_save_scatter_plot(
2665 base, other, plot_filename, title, one_to_one, model_names
2666 )
2668 # Add list of plots to plot metadata.
2669 plot_index = _append_to_plot_index([plot_filename])
2671 # Make a page to display the plots.
2672 _make_plot_html_page(plot_index)
2674 return iris.cube.CubeList([base, other])
2677def scatter_plot(
2678 cube_x: iris.cube.Cube | iris.cube.CubeList,
2679 cube_y: iris.cube.Cube | iris.cube.CubeList,
2680 filename: str = None,
2681 one_to_one: bool = True,
2682 **kwargs,
2683) -> iris.cube.CubeList:
2684 """Plot a scatter plot between two variables.
2686 Both cubes must be 1D.
2688 Parameters
2689 ----------
2690 cube_x: Cube | CubeList
2691 1 dimensional Cube of the data to plot on y-axis.
2692 cube_y: Cube | CubeList
2693 1 dimensional Cube of the data to plot on x-axis.
2694 filename: str, optional
2695 Filename of the plot to write.
2696 one_to_one: bool, optional
2697 If True a 1:1 line is plotted; if False it is not. Default is True.
2699 Returns
2700 -------
2701 cubes: CubeList
2702 CubeList of the original x and y cubes for further processing.
2704 Raises
2705 ------
2706 ValueError
2707 If the cube doesn't have the right dimensions and cubes not the same
2708 size.
2709 TypeError
2710 If the cube isn't a single cube.
2712 Notes
2713 -----
2714 Scatter plots are used for determining if there is a relationship between
2715 two variables. Positive relations have a slope going from bottom left to top
2716 right; Negative relations have a slope going from top left to bottom right.
2717 """
2718 # Iterate over all cubes in cube or CubeList and plot.
2719 for cube_iter in iter_maybe(cube_x):
2720 # Check cubes are correct shape.
2721 cube_iter = check_single_cube(cube_iter)
2722 if cube_iter.ndim > 1:
2723 raise ValueError("cube_x must be 1D.")
2725 # Iterate over all cubes in cube or CubeList and plot.
2726 for cube_iter in iter_maybe(cube_y):
2727 # Check cubes are correct shape.
2728 cube_iter = check_single_cube(cube_iter)
2729 if cube_iter.ndim > 1:
2730 raise ValueError("cube_y must be 1D.")
2732 # Ensure we have a name for the plot file.
2733 recipe_title = get_recipe_metadata().get("title", "Untitled")
2734 title = f"{recipe_title}"
2736 if filename is None:
2737 filename = slugify(recipe_title)
2739 # Add file extension.
2740 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2742 # Do the actual plotting.
2743 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2745 # Add list of plots to plot metadata.
2746 plot_index = _append_to_plot_index([plot_filename])
2748 # Make a page to display the plots.
2749 _make_plot_html_page(plot_index)
2751 return iris.cube.CubeList([cube_x, cube_y])
2754def vector_plot(
2755 cube_u: iris.cube.Cube,
2756 cube_v: iris.cube.Cube,
2757 filename: str = None,
2758 sequence_coordinate: str = "time",
2759 **kwargs,
2760) -> iris.cube.CubeList:
2761 """Plot a vector plot based on the input u and v components."""
2762 recipe_title = get_recipe_metadata().get("title", "Untitled")
2764 # Cubes must have a matching sequence coordinate.
2765 try:
2766 # Check that the u and v cubes have the same sequence coordinate.
2767 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2767 ↛ anywhereline 2767 didn't jump anywhere: it always raised an exception.
2768 raise ValueError("Coordinates do not match.")
2769 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2770 raise ValueError(
2771 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2772 ) from err
2774 # Create a plot for each value of the sequence coordinate.
2775 plot_index = []
2776 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2777 for cube_u_slice, cube_v_slice in zip(
2778 cube_u.slices_over(sequence_coordinate),
2779 cube_v.slices_over(sequence_coordinate),
2780 strict=True,
2781 ):
2782 # Format the coordinate value in a unit appropriate way.
2783 seq_coord = cube_u_slice.coord(sequence_coordinate)
2784 plot_title, plot_filename = _set_title_and_filename(
2785 seq_coord, nplot, recipe_title, filename
2786 )
2788 # Do the actual plotting.
2789 _plot_and_save_vector_plot(
2790 cube_u_slice,
2791 cube_v_slice,
2792 filename=plot_filename,
2793 title=plot_title,
2794 method="pcolormesh",
2795 )
2796 plot_index.append(plot_filename)
2798 # Add list of plots to plot metadata.
2799 complete_plot_index = _append_to_plot_index(plot_index)
2801 # Make a page to display the plots.
2802 _make_plot_html_page(complete_plot_index)
2804 return iris.cube.CubeList([cube_u, cube_v])
2807def plot_histogram_series(
2808 cubes: iris.cube.Cube | iris.cube.CubeList,
2809 filename: str = None,
2810 sequence_coordinate: str = "time",
2811 stamp_coordinate: str = "realization",
2812 single_plot: bool = False,
2813 **kwargs,
2814) -> iris.cube.Cube | iris.cube.CubeList:
2815 """Plot a histogram plot for each vertical level provided.
2817 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2818 is present then a sequence of plots will be produced using the time slider
2819 functionality to scroll through histograms against time. If a
2820 stamp_coordinate is present then postage stamp plots will be produced. If
2821 stamp_coordinate and single_plot is True, all postage stamp plots will be
2822 plotted in a single plot instead of separate postage stamp plots.
2824 Parameters
2825 ----------
2826 cubes: Cube | iris.cube.CubeList
2827 Iris cube or CubeList of the data to plot. It should have a single dimension other
2828 than the stamp coordinate.
2829 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2830 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2831 filename: str, optional
2832 Name of the plot to write, used as a prefix for plot sequences. Defaults
2833 to the recipe name.
2834 sequence_coordinate: str, optional
2835 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2836 This coordinate must exist in the cube and will be used for the time
2837 slider.
2838 stamp_coordinate: str, optional
2839 Coordinate about which to plot postage stamp plots. Defaults to
2840 ``"realization"``.
2841 single_plot: bool, optional
2842 If True, all postage stamp plots will be plotted in a single plot. If
2843 False, each postage stamp plot will be plotted separately. Is only valid
2844 if stamp_coordinate exists and has more than a single point.
2846 Returns
2847 -------
2848 iris.cube.Cube | iris.cube.CubeList
2849 The original Cube or CubeList (so further operations can be applied).
2850 Plotted data.
2852 Raises
2853 ------
2854 ValueError
2855 If the cube doesn't have the right dimensions.
2856 TypeError
2857 If the cube isn't a Cube or CubeList.
2858 """
2859 recipe_title = get_recipe_metadata().get("title", "Untitled")
2861 cubes = iter_maybe(cubes)
2862 # Ensure we have a name for the plot file.
2863 if filename is None:
2864 filename = slugify(recipe_title)
2866 # Internal plotting function.
2867 plotting_func = _plot_and_save_histogram_series
2869 num_models = get_num_models(cubes)
2871 validate_cube_shape(cubes, num_models)
2873 # If several histograms are plotted, check sequence_coordinate
2874 check_sequence_coordinate(cubes, sequence_coordinate)
2876 # Get axis minimum and maximum from levels information.
2877 # If no levels set, derive minima and maxima from data in CubeList.
2878 vmin, vmax = _set_axis_range(cubes)
2880 # Make postage stamp plots if stamp_coordinate exists and has more than a
2881 # single point. If single_plot is True:
2882 # -- all postage stamp plots will be plotted in a single plot instead of
2883 # separate postage stamp plots.
2884 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2885 # produced per single model only
2886 if num_models == 1:
2887 if (
2888 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2889 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2890 ):
2891 if single_plot:
2892 plotting_func = (
2893 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2894 )
2895 else:
2896 plotting_func = _plot_and_save_postage_stamp_histogram_series
2897 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2898 else:
2899 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2901 plot_index = []
2902 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2903 # Create a plot for each value of the sequence coordinate. Allowing for
2904 # multiple cubes in a CubeList to be plotted in the same plot for similar
2905 # sequence values. Passing a CubeList into the internal plotting function
2906 # for similar values of the sequence coordinate. cube_slice can be an
2907 # iris.cube.Cube or an iris.cube.CubeList.
2908 for cube_slice in cube_iterables:
2909 single_cube = cube_slice
2910 if isinstance(cube_slice, iris.cube.CubeList):
2911 single_cube = cube_slice[0]
2913 # Ensure valid stamp coordinate in cube dimensions
2914 if stamp_coordinate == "realization":
2915 stamp_coordinate = check_stamp_coordinate(single_cube)
2916 # Set plot titles and filename, based on sequence coordinate
2917 seq_coord = single_cube.coord(sequence_coordinate)
2918 # Use time coordinate in title and filename if single histogram output.
2919 if sequence_coordinate == "realization" and nplot == 1:
2920 seq_coord = single_cube.coord("time")
2921 plot_title, plot_filename = _set_title_and_filename(
2922 seq_coord, nplot, recipe_title, filename
2923 )
2925 # Do the actual plotting.
2926 plotting_func(
2927 cube_slice,
2928 filename=plot_filename,
2929 stamp_coordinate=stamp_coordinate,
2930 title=plot_title,
2931 vmin=vmin,
2932 vmax=vmax,
2933 )
2934 plot_index.append(plot_filename)
2936 # Add list of plots to plot metadata.
2937 complete_plot_index = _append_to_plot_index(plot_index)
2939 # Make a page to display the plots.
2940 _make_plot_html_page(complete_plot_index)
2942 return cubes
2945def _plot_and_save_postage_stamp_power_spectrum_series(
2946 cubes: iris.cube.Cube,
2947 coords: list[iris.coords.Coord],
2948 stamp_coordinate: str,
2949 filename: str,
2950 title: str,
2951 series_coordinate: str = None,
2952 **kwargs,
2953):
2954 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2956 Parameters
2957 ----------
2958 cubes: Cube or CubeList
2959 Cube or Cubelist of the power spectrum data.
2960 coords: list[Coord]
2961 Coordinates to plot on the x-axis, one per cube.
2962 stamp_coordinate: str
2963 Coordinate that becomes different plots.
2964 filename: str
2965 Filename of the plot to write.
2966 title: str
2967 Plot title.
2968 series_coordinate: str, optional
2969 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2971 """
2972 # Use the smallest square grid that will fit the members.
2973 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
2975 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
2976 model_colors_map = get_model_colors_map(cubes)
2977 # ax = plt.gca()
2978 # Make a subplot for each member.
2979 for member, subplot in zip(
2980 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
2981 ):
2982 ax = plt.subplot(grid_size, grid_size, subplot)
2984 # Store min/max ranges.
2985 y_levels = []
2987 line_marker = None
2988 line_width = 1
2990 for cube in iter_maybe(member):
2991 xcoord = _select_series_coord(cube, series_coordinate)
2992 xname = xcoord.points
2994 yfield = cube.data # power spectrum
2995 label = None
2996 color = "black"
2997 if model_colors_map: 2997 ↛ 2998line 2997 didn't jump to line 2998 because the condition on line 2997 was never true
2998 label = cube.attributes.get("model_name")
2999 color = model_colors_map.get(label)
3001 if member.coord(stamp_coordinate).points == [0]:
3002 ax.plot(
3003 xname,
3004 yfield,
3005 color=color,
3006 marker=line_marker,
3007 ls="-",
3008 lw=line_width,
3009 label=f"{label} (control)"
3010 if len(cube.coord(stamp_coordinate).points) > 1
3011 else label,
3012 )
3013 # Label with member if part of an ensemble and not the control.
3014 else:
3015 ax.plot(
3016 xname,
3017 yfield,
3018 color=color,
3019 ls="-",
3020 lw=1.5,
3021 alpha=0.75,
3022 label=f"{label} (member)",
3023 )
3025 # Calculate the global min/max if multiple cubes are given.
3026 _, levels, _ = colorbar_map_levels(cube, axis="y")
3027 if levels is not None: 3027 ↛ 3028line 3027 didn't jump to line 3028 because the condition on line 3027 was never true
3028 y_levels.append(min(levels))
3029 y_levels.append(max(levels))
3031 # Add some labels and tweak the style.
3032 title = f"{title}"
3033 ax.set_title(title, fontsize=16)
3035 # Set appropriate x-axis label based on coordinate
3036 if series_coordinate == "wavelength" or ( 3036 ↛ 3039line 3036 didn't jump to line 3039 because the condition on line 3036 was never true
3037 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3038 ):
3039 ax.set_xlabel("Wavelength (km)", fontsize=14)
3040 elif series_coordinate == "physical_wavenumber" or ( 3040 ↛ 3045line 3040 didn't jump to line 3045 because the condition on line 3040 was always true
3041 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3042 ):
3043 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3044 else: # frequency or check units
3045 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3046 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3047 else:
3048 ax.set_xlabel("Wavenumber", fontsize=14)
3050 ax.set_ylabel("Power Spectral Density", fontsize=14)
3051 ax.tick_params(axis="both", labelsize=12)
3053 # Set log-log scale
3054 ax.set_xscale("log")
3055 ax.set_yscale("log")
3057 # Add gridlines
3058 ax.grid(linestyle="--", color="grey", linewidth=1)
3059 # Ientify unique labels for legend
3060 handles = list(
3061 {
3062 label: handle
3063 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3064 }.values()
3065 )
3066 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
3068 ax = plt.gca()
3069 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3071 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3072 logging.info("Saved histogram postage stamp plot to %s", filename)
3073 plt.close(fig)
3076def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3077 cubes: iris.cube.Cube,
3078 coords: list[iris.coords.Coord],
3079 stamp_coordinate: str,
3080 filename: str,
3081 title: str,
3082 series_coordinate: str = None,
3083 **kwargs,
3084):
3085 """Plot and save power spectra for ensemble members in single plot.
3087 Parameters
3088 ----------
3089 cubes: Cube or CubeList
3090 Cube or Cubelist of the power spectrum data.
3091 coords: list[Coord]
3092 Coordinates to plot on the x-axis, one per cube.
3093 stamp_coordinate: str
3094 Coordinate that becomes different plots.
3095 filename: str
3096 Filename of the plot to write.
3097 title: str
3098 Plot title.
3099 series_coordinate: str, optional
3100 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3102 """
3103 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3104 model_colors_map = get_model_colors_map(cubes)
3106 line_marker = None
3107 line_width = 1
3109 # Compute ensemble statistics to show spread
3110 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3111 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3112 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3114 xcoord_global = mean_cube.coord(series_coordinate)
3115 x_global = xcoord_global.points
3117 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3118 xcoord = _select_series_coord(member, series_coordinate)
3119 xname = xcoord.points
3121 yfield = member.data # power spectrum
3122 color = "black"
3123 if model_colors_map: 3123 ↛ 3127line 3123 didn't jump to line 3127 because the condition on line 3123 was always true
3124 label = member.attributes.get("model_name") if i == 0 else None
3125 color = model_colors_map.get(label)
3127 if member.coord(stamp_coordinate).points == [0]:
3128 ax.plot(
3129 xname,
3130 yfield,
3131 color=color,
3132 marker=line_marker,
3133 ls="-",
3134 lw=line_width,
3135 label=f"{label} (control)"
3136 if len(member.coord(stamp_coordinate).points) > 1
3137 else label,
3138 )
3139 # Label with member number if part of an ensemble and not the control.
3140 else:
3141 ax.plot(
3142 xname,
3143 yfield,
3144 color=color,
3145 ls="-",
3146 lw=1.5,
3147 alpha=0.75,
3148 label=label,
3149 )
3151 # Set appropriate x-axis label based on coordinate
3152 if series_coordinate == "wavelength" or ( 3152 ↛ 3155line 3152 didn't jump to line 3155 because the condition on line 3152 was never true
3153 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3154 ):
3155 ax.set_xlabel("Wavelength (km)", fontsize=14)
3156 elif series_coordinate == "physical_wavenumber" or ( 3156 ↛ 3161line 3156 didn't jump to line 3161 because the condition on line 3156 was always true
3157 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3158 ):
3159 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3160 else: # frequency or check units
3161 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3162 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3163 else:
3164 ax.set_xlabel("Wavenumber", fontsize=14)
3166 # Add ensemble spread shading
3167 ax.fill_between(
3168 x_global,
3169 min_cube.data,
3170 max_cube.data,
3171 color="grey",
3172 alpha=0.3,
3173 label="Ensemble spread",
3174 )
3176 # Add ensemble mean line
3177 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3179 ax.set_ylabel("Power Spectral Density", fontsize=14)
3180 ax.tick_params(axis="both", labelsize=12)
3182 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3183 # Set log-log scale
3184 ax.set_xscale("log")
3185 ax.set_yscale("log")
3187 # Add gridlines
3188 ax.grid(linestyle="--", color="grey", linewidth=1)
3189 # Identify unique labels for legend
3190 handles = list(
3191 {
3192 label: handle
3193 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3194 }.values()
3195 )
3196 ax.legend(handles=handles, loc="best", ncol=1, frameon=False, fontsize=16)
3198 # Add a legend
3199 ax.legend(fontsize=16)
3201 # Figure title.
3202 ax.set_title(title, fontsize=16)
3204 # Save the figure to a file
3205 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3207 # Close the figure
3208 plt.close(fig)