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