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