Coverage for src/CSET/operators/plot.py: 78%
1011 statements
« prev ^ index » next coverage.py v7.15.1, created at 2026-07-15 14:27 +0000
1# © Crown copyright, Met Office (2022-2025) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators to produce various kinds of plots."""
17import fcntl
18import importlib.resources
19import itertools
20import json
21import logging
22import math
23import os
24from typing import Literal
26import cartopy.crs as ccrs
27import cartopy.feature as cfeature
28import iris
29import iris.coords
30import iris.cube
31import iris.exceptions
32import iris.plot as iplt
33import matplotlib as mpl
34import matplotlib.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: 201 ↛ 202line 201 didn't jump to line 202 because the condition on line 201 was never true
202 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
203 elif ymin < -80.0 and ymax < -20.0: 203 ↛ 204line 203 didn't jump to line 204 because the condition on line 203 was never true
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( 228 ↛ 231line 228 didn't jump to line 231 because the condition on line 228 was never true
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"],
485 overlay_cube: iris.cube.Cube | None = None,
486 contour_cube: iris.cube.Cube | None = None,
487 **kwargs,
488):
489 """Plot and save a spatial plot.
491 Parameters
492 ----------
493 cube: Cube
494 2 dimensional (lat and lon) Cube of the data to plot.
495 filename: str
496 Filename of the plot to write.
497 title: str
498 Plot title.
499 method: "contourf" | "pcolormesh"
500 The plotting method to use.
501 overlay_cube: Cube, optional
502 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
503 contour_cube: Cube, optional
504 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
505 """
506 # Setup plot details, size, resolution, etc.
507 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
509 # Specify the color bar
510 cmap, levels, norm = colorbar_map_levels(cube)
512 # If overplotting, set required colorbars
513 if overlay_cube:
514 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
515 if contour_cube:
516 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
518 # Setup plot map projection, extent and coastlines and borderlines.
519 axes = _setup_spatial_map(cube, fig, cmap)
521 # Set colorscale bounds
522 try:
523 vmin = min(levels)
524 vmax = max(levels)
525 except TypeError:
526 vmin, vmax = None, None
527 # Ensure to use norm and not vmin/vmax if levels are defined.
528 if norm is not None:
529 vmin = None
530 vmax = None
531 logging.debug("Plotting using defined levels.")
533 # Plot the field.
534 if method == "contourf":
535 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
536 elif method == "pcolormesh":
537 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
538 else:
539 raise ValueError(f"Unknown plotting method: {method}")
541 # Overplot overlay field, if required
542 if overlay_cube:
543 try:
544 over_vmin = min(over_levels)
545 over_vmax = max(over_levels)
546 except TypeError:
547 over_vmin, over_vmax = None, None
548 if over_norm is not None: 548 ↛ 549line 548 didn't jump to line 549 because the condition on line 548 was never true
549 over_vmin = None
550 over_vmax = None
551 overlay = iplt.pcolormesh(
552 overlay_cube,
553 cmap=over_cmap,
554 norm=over_norm,
555 alpha=0.8,
556 vmin=over_vmin,
557 vmax=over_vmax,
558 )
559 # Overplot contour field, if required, with contour labelling.
560 if contour_cube:
561 contour = iplt.contour(
562 contour_cube,
563 colors="darkgray",
564 levels=cntr_levels,
565 norm=cntr_norm,
566 alpha=0.5,
567 linestyles="--",
568 linewidths=1,
569 )
570 plt.clabel(contour)
572 # Check to see if transect, and if so, adjust y axis.
573 if is_transect(cube):
574 if "pressure" in [coord.name() for coord in cube.coords()]:
575 axes.invert_yaxis()
576 axes.set_yscale("log")
577 axes.set_ylim(1100, 100)
578 # If both model_level_number and level_height exists, iplt can construct
579 # plot as a function of height above orography (NOT sea level).
580 elif {"model_level_number", "level_height"}.issubset( 580 ↛ 585line 580 didn't jump to line 585 because the condition on line 580 was always true
581 {coord.name() for coord in cube.coords()}
582 ):
583 axes.set_yscale("log")
585 axes.set_title(
586 f"{title}\n"
587 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
588 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
589 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
590 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
591 fontsize=16,
592 )
594 # Inset code
595 axins = inset_axes(
596 axes,
597 width="20%",
598 height="20%",
599 loc="upper right",
600 axes_class=GeoAxes,
601 axes_kwargs={"map_projection": ccrs.PlateCarree()},
602 )
604 # Slightly transparent to reduce plot blocking.
605 axins.patch.set_alpha(0.4)
607 axins.coastlines(resolution="50m")
608 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
610 SLat, SLon, ELat, ELon = (
611 float(coord) for coord in cube.attributes["transect_coords"].split("_")
612 )
614 # Draw line between them
615 axins.plot(
616 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
617 )
619 # Plot points (note: lon, lat order for Cartopy)
620 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
621 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
623 lon_min, lon_max = sorted([SLon, ELon])
624 lat_min, lat_max = sorted([SLat, ELat])
626 # Midpoints
627 lon_mid = (lon_min + lon_max) / 2
628 lat_mid = (lat_min + lat_max) / 2
630 # Maximum half-range
631 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
632 if half_range == 0: # points identical → provide small default 632 ↛ 636line 632 didn't jump to line 636 because the condition on line 632 was always true
633 half_range = 1
635 # Set square extent
636 axins.set_extent(
637 [
638 lon_mid - half_range,
639 lon_mid + half_range,
640 lat_mid - half_range,
641 lat_mid + half_range,
642 ],
643 crs=ccrs.PlateCarree(),
644 )
646 # Ensure square aspect
647 axins.set_aspect("equal")
649 else:
650 # Add title.
651 axes.set_title(title, fontsize=16)
653 # Adjust padding if spatial plot or transect
654 if is_transect(cube):
655 yinfopad = -0.1
656 ycbarpad = 0.1
657 else:
658 yinfopad = 0.01
659 ycbarpad = 0.042
661 # Add watermark with min/max/mean. Currently not user togglable.
662 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
663 axes.annotate(
664 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
665 xy=(0.025, yinfopad),
666 xycoords="axes fraction",
667 xytext=(-5, 5),
668 textcoords="offset points",
669 ha="left",
670 va="bottom",
671 size=11,
672 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
673 )
675 # Add secondary colour bar for overlay_cube field if required.
676 if overlay_cube:
677 cbarB = fig.colorbar(
678 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
679 )
680 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
681 # add ticks and tick_labels for every levels if less than 20 levels exist
682 if over_levels is not None and len(over_levels) < 20: 682 ↛ 683line 682 didn't jump to line 683 because the condition on line 682 was never true
683 cbarB.set_ticks(over_levels)
684 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
685 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
686 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
687 logging.debug("Set secondary colorbar ticks and labels.")
689 # Add main colour bar.
690 cbar = fig.colorbar(
691 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
692 )
694 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
695 # add ticks and tick_labels for every levels if less than 20 levels exist
696 if levels is not None and len(levels) < 20:
697 cbar.set_ticks(levels)
698 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
699 if "rainfall" or "snowfall" or "visibility" in cube.name(): 699 ↛ 701line 699 didn't jump to line 701 because the condition on line 699 was always true
700 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
701 logging.debug("Set colorbar ticks and labels.")
703 # Save plot.
704 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
705 logging.info("Saved spatial plot to %s", filename)
706 plt.close(fig)
709def _plot_and_save_postage_stamp_spatial_plot(
710 cube: iris.cube.Cube,
711 filename: str,
712 stamp_coordinate: str,
713 title: str,
714 method: Literal["contourf", "pcolormesh"],
715 overlay_cube: iris.cube.Cube | None = None,
716 contour_cube: iris.cube.Cube | None = None,
717 **kwargs,
718):
719 """Plot postage stamp spatial plots from an ensemble.
721 Parameters
722 ----------
723 cube: Cube
724 Iris cube of data to be plotted. It must have the stamp coordinate.
725 filename: str
726 Filename of the plot to write.
727 stamp_coordinate: str
728 Coordinate that becomes different plots.
729 method: "contourf" | "pcolormesh"
730 The plotting method to use.
731 overlay_cube: Cube, optional
732 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
733 contour_cube: Cube, optional
734 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
736 Raises
737 ------
738 ValueError
739 If the cube doesn't have the right dimensions.
740 """
741 # Use the smallest square grid that will fit the members.
742 nmember = len(cube.coord(stamp_coordinate).points)
743 grid_rows = int(math.sqrt(nmember))
744 grid_size = math.ceil(nmember / grid_rows)
746 fig = plt.figure(
747 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
748 )
750 # Specify the color bar
751 cmap, levels, norm = colorbar_map_levels(cube)
752 # If overplotting, set required colorbars
753 if overlay_cube: 753 ↛ 754line 753 didn't jump to line 754 because the condition on line 753 was never true
754 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
755 if contour_cube: 755 ↛ 756line 755 didn't jump to line 756 because the condition on line 755 was never true
756 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
758 # Make a subplot for each member.
759 for member, subplot in zip(
760 cube.slices_over(stamp_coordinate),
761 range(1, grid_size * grid_rows + 1),
762 strict=False,
763 ):
764 # Setup subplot map projection, extent and coastlines and borderlines.
765 axes = _setup_spatial_map(
766 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
767 )
768 if method == "contourf":
769 # Filled contour plot of the field.
770 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
771 elif method == "pcolormesh":
772 if levels is not None:
773 vmin = min(levels)
774 vmax = max(levels)
775 else:
776 raise TypeError("Unknown vmin and vmax range.")
777 vmin, vmax = None, None
778 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
779 # if levels are defined.
780 if norm is not None: 780 ↛ 781line 780 didn't jump to line 781 because the condition on line 780 was never true
781 vmin = None
782 vmax = None
783 # pcolormesh plot of the field.
784 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
785 else:
786 raise ValueError(f"Unknown plotting method: {method}")
788 # Overplot overlay field, if required
789 if overlay_cube: 789 ↛ 790line 789 didn't jump to line 790 because the condition on line 789 was never true
790 try:
791 over_vmin = min(over_levels)
792 over_vmax = max(over_levels)
793 except TypeError:
794 over_vmin, over_vmax = None, None
795 if over_norm is not None:
796 over_vmin = None
797 over_vmax = None
798 iplt.pcolormesh(
799 overlay_cube[member.coord(stamp_coordinate).points[0]],
800 cmap=over_cmap,
801 norm=over_norm,
802 alpha=0.6,
803 vmin=over_vmin,
804 vmax=over_vmax,
805 )
806 # Overplot contour field, if required
807 if contour_cube: 807 ↛ 808line 807 didn't jump to line 808 because the condition on line 807 was never true
808 iplt.contour(
809 contour_cube[member.coord(stamp_coordinate).points[0]],
810 colors="darkgray",
811 levels=cntr_levels,
812 norm=cntr_norm,
813 alpha=0.6,
814 linestyles="--",
815 linewidths=1,
816 )
817 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
818 axes.set_title(f"{mtitle}")
820 # Put the shared colorbar in its own axes.
821 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
822 colorbar = fig.colorbar(
823 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
824 )
825 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
827 # Overall figure title.
828 fig.suptitle(title, fontsize=16)
830 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
831 logging.info("Saved contour postage stamp plot to %s", filename)
832 plt.close(fig)
835def _plot_and_save_line_series(
836 cubes: iris.cube.CubeList,
837 coords: list[iris.coords.Coord],
838 ensemble_coord: str,
839 filename: str,
840 title: str,
841 **kwargs,
842):
843 """Plot and save a 1D line series.
845 Parameters
846 ----------
847 cubes: Cube or CubeList
848 Cube or CubeList containing the cubes to plot on the y-axis.
849 coords: list[Coord]
850 Coordinates to plot on the x-axis, one per cube.
851 ensemble_coord: str
852 Ensemble coordinate in the cube.
853 filename: str
854 Filename of the plot to write.
855 title: str
856 Plot title.
857 """
858 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
860 model_colors_map = get_model_colors_map(cubes)
862 # Store min/max ranges.
863 y_levels = []
865 # Check match-up across sequence coords gives consistent sizes
866 validate_cubes_coords(cubes, coords)
868 for cube, coord in zip(cubes, coords, strict=True):
869 label = None
870 color = "black"
871 if model_colors_map:
872 label = cube.attributes.get("model_name")
873 color = model_colors_map.get(label)
874 for cube_slice in cube.slices_over(ensemble_coord):
875 # Label with (control) if part of an ensemble or not otherwise.
876 if cube_slice.coord(ensemble_coord).points == [0]:
877 iplt.plot(
878 coord,
879 cube_slice,
880 color=color,
881 marker="o",
882 ls="-",
883 lw=3,
884 label=f"{label} (control)"
885 if len(cube.coord(ensemble_coord).points) > 1
886 else label,
887 )
888 # Label with (perturbed) if part of an ensemble and not the control.
889 else:
890 iplt.plot(
891 coord,
892 cube_slice,
893 color=color,
894 ls="-",
895 lw=1.5,
896 alpha=0.75,
897 label=f"{label} (member)",
898 )
900 # Calculate the global min/max if multiple cubes are given.
901 _, levels, _ = colorbar_map_levels(cube, axis="y")
902 if levels is not None: 902 ↛ 903line 902 didn't jump to line 903 because the condition on line 902 was never true
903 y_levels.append(min(levels))
904 y_levels.append(max(levels))
906 # Get the current axes.
907 ax = plt.gca()
909 # Add some labels and tweak the style.
910 # check if cubes[0] works for single cube if not CubeList
911 if coords[0].name() == "time":
912 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
913 else:
914 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
915 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
916 ax.set_title(title, fontsize=16)
918 ax.ticklabel_format(axis="y", useOffset=False)
919 ax.tick_params(axis="x", labelrotation=15)
920 ax.tick_params(axis="both", labelsize=12)
922 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
923 if y_levels: 923 ↛ 924line 923 didn't jump to line 924 because the condition on line 923 was never true
924 ax.set_ylim(min(y_levels), max(y_levels))
925 # Add zero line.
926 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
927 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
928 logging.debug(
929 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
930 )
931 else:
932 ax.autoscale()
934 # Add gridlines
935 ax.grid(linestyle="--", color="grey", linewidth=1)
936 # Ientify unique labels for legend
937 handles = list(
938 {
939 label: handle
940 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
941 }.values()
942 )
943 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
945 # Save plot.
946 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
947 logging.info("Saved line plot to %s", filename)
948 plt.close(fig)
951def _plot_and_save_line_power_spectrum_series(
952 cubes: iris.cube.Cube | iris.cube.CubeList,
953 coords: list[iris.coords.Coord],
954 ensemble_coord: str,
955 filename: str,
956 title: str,
957 series_coordinate: str,
958 **kwargs,
959):
960 """Plot and save a 1D line series.
962 Parameters
963 ----------
964 cubes: Cube or CubeList
965 Cube or CubeList containing the cubes to plot on the y-axis.
966 coords: list[Coord]
967 Coordinates to plot on the x-axis, one per cube.
968 ensemble_coord: str
969 Ensemble coordinate in the cube.
970 filename: str
971 Filename of the plot to write.
972 title: str
973 Plot title.
974 series_coordinate: str
975 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
976 """
977 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
978 model_colors_map = get_model_colors_map(cubes)
979 ax = plt.gca()
981 # Store min/max ranges.
982 y_levels = []
984 line_marker = None
985 line_width = 1
987 for cube in iter_maybe(cubes):
988 # next 2 lines replace chunk of code.
989 xcoord = _select_series_coord(cube, series_coordinate)
990 xname = xcoord.points
992 yfield = cube.data # power spectrum
993 label = None
994 color = "black"
995 if model_colors_map: 995 ↛ 998line 995 didn't jump to line 998 because the condition on line 995 was always true
996 label = cube.attributes.get("model_name")
997 color = model_colors_map.get(label)
998 for cube_slice in cube.slices_over(ensemble_coord):
999 # Label with (control) if part of an ensemble or not otherwise.
1000 if cube_slice.coord(ensemble_coord).points == [0]: 1000 ↛ 1014line 1000 didn't jump to line 1014 because the condition on line 1000 was always true
1001 ax.plot(
1002 xname,
1003 yfield,
1004 color=color,
1005 marker=line_marker,
1006 ls="-",
1007 lw=line_width,
1008 label=f"{label} (control)"
1009 if len(cube.coord(ensemble_coord).points) > 1
1010 else label,
1011 )
1012 # Label with (perturbed) if part of an ensemble and not the control.
1013 else:
1014 ax.plot(
1015 xname,
1016 yfield,
1017 color=color,
1018 ls="-",
1019 lw=1.5,
1020 alpha=0.75,
1021 label=f"{label} (member)",
1022 )
1024 # Calculate the global min/max if multiple cubes are given.
1025 _, levels, _ = colorbar_map_levels(cube, axis="y")
1026 if levels is not None: 1026 ↛ 1027line 1026 didn't jump to line 1027 because the condition on line 1026 was never true
1027 y_levels.append(min(levels))
1028 y_levels.append(max(levels))
1030 # Add some labels and tweak the style.
1032 title = f"{title}"
1033 ax.set_title(title, fontsize=16)
1035 # Set appropriate x-axis label based on coordinate
1036 if series_coordinate == "wavelength" or ( 1036 ↛ 1039line 1036 didn't jump to line 1039 because the condition on line 1036 was never true
1037 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1038 ):
1039 ax.set_xlabel("Wavelength (km)", fontsize=14)
1040 elif series_coordinate == "physical_wavenumber" or ( 1040 ↛ 1043line 1040 didn't jump to line 1043 because the condition on line 1040 was never true
1041 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1042 ):
1043 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1044 else: # frequency or check units
1045 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1045 ↛ 1046line 1045 didn't jump to line 1046 because the condition on line 1045 was never true
1046 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1047 else:
1048 ax.set_xlabel("Wavenumber", fontsize=14)
1050 ax.set_ylabel("Power Spectral Density", fontsize=14)
1051 ax.tick_params(axis="both", labelsize=12)
1053 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1055 # Set log-log scale
1056 ax.set_xscale("log")
1057 ax.set_yscale("log")
1059 # Add gridlines
1060 ax.grid(linestyle="--", color="grey", linewidth=1)
1061 # Ientify unique labels for legend
1062 handles = list(
1063 {
1064 label: handle
1065 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1066 }.values()
1067 )
1068 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1070 # Save plot.
1071 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1072 logging.info("Saved line plot to %s", filename)
1073 plt.close(fig)
1076def _plot_and_save_vertical_line_series(
1077 cubes: iris.cube.CubeList,
1078 coords: list[iris.coords.Coord],
1079 ensemble_coord: str,
1080 filename: str,
1081 series_coordinate: str,
1082 title: str,
1083 vmin: float,
1084 vmax: float,
1085 **kwargs,
1086):
1087 """Plot and save a 1D line series in vertical.
1089 Parameters
1090 ----------
1091 cubes: CubeList
1092 1 dimensional Cube or CubeList of the data to plot on x-axis.
1093 coord: list[Coord]
1094 Coordinates to plot on the y-axis, one per cube.
1095 ensemble_coord: str
1096 Ensemble coordinate in the cube.
1097 filename: str
1098 Filename of the plot to write.
1099 series_coordinate: str
1100 Coordinate to use as vertical axis.
1101 title: str
1102 Plot title.
1103 vmin: float
1104 Minimum value for the x-axis.
1105 vmax: float
1106 Maximum value for the x-axis.
1107 """
1108 # plot the vertical pressure axis using log scale
1109 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1111 model_colors_map = get_model_colors_map(cubes)
1113 # Check match-up across sequence coords gives consistent sizes
1114 validate_cubes_coords(cubes, coords)
1116 for cube, coord in zip(cubes, coords, strict=True):
1117 label = None
1118 color = "black"
1119 if model_colors_map: 1119 ↛ 1120line 1119 didn't jump to line 1120 because the condition on line 1119 was never true
1120 label = cube.attributes.get("model_name")
1121 color = model_colors_map.get(label)
1123 for cube_slice in cube.slices_over(ensemble_coord):
1124 # If ensemble data given plot control member with (control)
1125 # unless single forecast.
1126 if cube_slice.coord(ensemble_coord).points == [0]:
1127 iplt.plot(
1128 cube_slice,
1129 coord,
1130 color=color,
1131 marker="o",
1132 ls="-",
1133 lw=3,
1134 label=f"{label} (control)"
1135 if len(cube.coord(ensemble_coord).points) > 1
1136 else label,
1137 )
1138 # If ensemble data given plot perturbed members with (perturbed).
1139 else:
1140 iplt.plot(
1141 cube_slice,
1142 coord,
1143 color=color,
1144 ls="-",
1145 lw=1.5,
1146 alpha=0.75,
1147 label=f"{label} (member)",
1148 )
1150 # Get the current axis
1151 ax = plt.gca()
1153 # Special handling for pressure level data.
1154 if series_coordinate == "pressure": 1154 ↛ 1176line 1154 didn't jump to line 1176 because the condition on line 1154 was always true
1155 # Invert y-axis and set to log scale.
1156 ax.invert_yaxis()
1157 ax.set_yscale("log")
1159 # Define y-ticks and labels for pressure log axis.
1160 y_tick_labels = [
1161 "1000",
1162 "850",
1163 "700",
1164 "500",
1165 "300",
1166 "200",
1167 "100",
1168 ]
1169 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1171 # Set y-axis limits and ticks.
1172 ax.set_ylim(1100, 100)
1174 # Test if series_coordinate is model level data. The UM data uses
1175 # model_level_number and lfric uses full_levels as coordinate.
1176 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1177 # Define y-ticks and labels for vertical axis.
1178 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1179 y_tick_labels = [str(int(i)) for i in y_ticks]
1180 ax.set_ylim(min(y_ticks), max(y_ticks))
1182 ax.set_yticks(y_ticks)
1183 ax.set_yticklabels(y_tick_labels)
1185 # Set x-axis limits.
1186 ax.set_xlim(vmin, vmax)
1187 # Mark y=0 if present in plot.
1188 if vmin < 0.0 and vmax > 0.0: 1188 ↛ 1189line 1188 didn't jump to line 1189 because the condition on line 1188 was never true
1189 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1191 # Add some labels and tweak the style.
1192 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1193 ax.set_xlabel(
1194 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1195 )
1196 ax.set_title(title, fontsize=16)
1197 ax.ticklabel_format(axis="x")
1198 ax.tick_params(axis="y")
1199 ax.tick_params(axis="both", labelsize=12)
1201 # Add gridlines
1202 ax.grid(linestyle="--", color="grey", linewidth=1)
1203 # Ientify unique labels for legend
1204 handles = list(
1205 {
1206 label: handle
1207 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1208 }.values()
1209 )
1210 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1212 # Save plot.
1213 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1214 logging.info("Saved line plot to %s", filename)
1215 plt.close(fig)
1218def _plot_and_save_scatter_plot(
1219 cube_x: iris.cube.Cube | iris.cube.CubeList,
1220 cube_y: iris.cube.Cube | iris.cube.CubeList,
1221 filename: str,
1222 title: str,
1223 one_to_one: bool,
1224 model_names: list[str] = None,
1225 **kwargs,
1226):
1227 """Plot and save a 2D scatter plot.
1229 Parameters
1230 ----------
1231 cube_x: Cube | CubeList
1232 1 dimensional Cube or CubeList of the data to plot on x-axis.
1233 cube_y: Cube | CubeList
1234 1 dimensional Cube or CubeList of the data to plot on y-axis.
1235 filename: str
1236 Filename of the plot to write.
1237 title: str
1238 Plot title.
1239 one_to_one: bool
1240 Whether a 1:1 line is plotted.
1241 """
1242 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1243 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1244 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1245 # over the pairs simultaneously.
1247 # Ensure cube_x and cube_y are iterable
1248 cube_x_iterable = iter_maybe(cube_x)
1249 cube_y_iterable = iter_maybe(cube_y)
1251 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1252 iplt.scatter(cube_x_iter, cube_y_iter)
1253 if one_to_one is True:
1254 plt.plot(
1255 [
1256 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1257 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1258 ],
1259 [
1260 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1261 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1262 ],
1263 "k",
1264 linestyle="--",
1265 )
1266 ax = plt.gca()
1268 # Add some labels and tweak the style.
1269 if model_names is None:
1270 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1271 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1272 else:
1273 # Add the model names, these should be order of base (x) and other (y).
1274 ax.set_xlabel(
1275 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1276 )
1277 ax.set_ylabel(
1278 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1279 )
1280 ax.set_title(title, fontsize=16)
1281 ax.ticklabel_format(axis="y", useOffset=False)
1282 ax.tick_params(axis="x", labelrotation=15)
1283 ax.tick_params(axis="both", labelsize=12)
1284 ax.autoscale()
1286 # Save plot.
1287 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1288 logging.info("Saved scatter plot to %s", filename)
1289 plt.close(fig)
1292def _plot_and_save_vector_plot(
1293 cube_u: iris.cube.Cube,
1294 cube_v: iris.cube.Cube,
1295 filename: str,
1296 title: str,
1297 method: Literal["contourf", "pcolormesh"],
1298 **kwargs,
1299):
1300 """Plot and save a 2D vector plot.
1302 Parameters
1303 ----------
1304 cube_u: Cube
1305 2 dimensional Cube of u component of the data.
1306 cube_v: Cube
1307 2 dimensional Cube of v component of the data.
1308 filename: str
1309 Filename of the plot to write.
1310 title: str
1311 Plot title.
1312 """
1313 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1314 # Create a cube containing the magnitude of the vector field.
1315 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1316 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1317 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1318 cube_vec_mag.rename(
1319 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1320 )
1322 # Specify the color bar
1323 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1325 # Setup plot map projection, extent and coastlines and borderlines.
1326 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1328 if method == "contourf":
1329 # Filled contour plot of the field.
1330 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1331 elif method == "pcolormesh":
1332 try:
1333 vmin = min(levels)
1334 vmax = max(levels)
1335 except TypeError:
1336 vmin, vmax = None, None
1337 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1338 # if levels are defined.
1339 if norm is not None:
1340 vmin = None
1341 vmax = None
1342 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1343 else:
1344 raise ValueError(f"Unknown plotting method: {method}")
1346 # Check to see if transect, and if so, adjust y axis.
1347 if is_transect(cube_vec_mag):
1348 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1349 axes.invert_yaxis()
1350 axes.set_yscale("log")
1351 axes.set_ylim(1100, 100)
1352 # If both model_level_number and level_height exists, iplt can construct
1353 # plot as a function of height above orography (NOT sea level).
1354 elif {"model_level_number", "level_height"}.issubset(
1355 {coord.name() for coord in cube_vec_mag.coords()}
1356 ):
1357 axes.set_yscale("log")
1359 axes.set_title(
1360 f"{title}\n"
1361 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1362 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1363 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1364 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1365 fontsize=16,
1366 )
1368 else:
1369 # Add title.
1370 axes.set_title(title, fontsize=16)
1372 # Add watermark with min/max/mean. Currently not user togglable.
1373 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1374 axes.annotate(
1375 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}",
1376 xy=(0.05, -0.05),
1377 xycoords="axes fraction",
1378 xytext=(-5, 5),
1379 textcoords="offset points",
1380 ha="right",
1381 va="bottom",
1382 size=11,
1383 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1384 )
1386 # Add colour bar.
1387 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1388 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1389 # add ticks and tick_labels for every levels if less than 20 levels exist
1390 if levels is not None and len(levels) < 20:
1391 cbar.set_ticks(levels)
1392 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1394 # 30 barbs along the longest axis of the plot, or a barb per point for data
1395 # with less than 30 points.
1396 step = max(max(cube_u.shape) // 30, 1)
1397 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1399 # Save plot.
1400 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1401 logging.info("Saved vector plot to %s", filename)
1402 plt.close(fig)
1405def _plot_and_save_histogram_series(
1406 cubes: iris.cube.Cube | iris.cube.CubeList,
1407 filename: str,
1408 title: str,
1409 vmin: float,
1410 vmax: float,
1411 **kwargs,
1412):
1413 """Plot and save a histogram series.
1415 Parameters
1416 ----------
1417 cubes: Cube or CubeList
1418 2 dimensional Cube or CubeList of the data to plot as histogram.
1419 filename: str
1420 Filename of the plot to write.
1421 title: str
1422 Plot title.
1423 vmin: float
1424 minimum for colorbar
1425 vmax: float
1426 maximum for colorbar
1427 """
1428 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1429 ax = plt.gca()
1431 model_colors_map = get_model_colors_map(cubes)
1433 # Set default that histograms will produce probability density function
1434 # at each bin (integral over range sums to 1).
1435 density = True
1437 for cube in iter_maybe(cubes):
1438 # Easier to check title (where var name originates)
1439 # than seeing if long names exist etc.
1440 # Exception case, where distribution better fits log scales/bins.
1441 if "surface_microphysical" in title:
1442 if "amount" in title: 1442 ↛ 1444line 1442 didn't jump to line 1444 because the condition on line 1442 was never true
1443 # Compute histogram following Klingaman et al. (2017): ASoP
1444 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1445 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1446 density = False
1447 else:
1448 bins = 10.0 ** (
1449 np.arange(-10, 27, 1) / 10.0
1450 ) # Suggestion from RMED toolbox.
1451 bins = np.insert(bins, 0, 0)
1452 ax.set_yscale("log")
1453 vmin = bins[1]
1454 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1455 ax.set_xscale("log")
1456 elif "lightning" in title:
1457 bins = [0, 1, 2, 3, 4, 5]
1458 else:
1459 bins = np.linspace(vmin, vmax, 51)
1460 logging.debug(
1461 "Plotting histogram with %s bins %s - %s.",
1462 np.size(bins),
1463 np.min(bins),
1464 np.max(bins),
1465 )
1467 # Reshape cube data into a single array to allow for a single histogram.
1468 # Otherwise we plot xdim histograms stacked.
1469 cube_data_1d = (cube.data).flatten()
1471 label = None
1472 color = "black"
1473 if model_colors_map: 1473 ↛ 1474line 1473 didn't jump to line 1474 because the condition on line 1473 was never true
1474 label = cube.attributes.get("model_name")
1475 color = model_colors_map[label]
1476 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1478 # Compute area under curve.
1479 if "surface_microphysical" in title and "amount" in title: 1479 ↛ 1480line 1479 didn't jump to line 1480 because the condition on line 1479 was never true
1480 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1481 x = x * bin_mean / x.sum()
1482 x = x[1:]
1483 y = y[1:]
1485 ax.plot(
1486 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1487 )
1489 # Add some labels and tweak the style.
1490 ax.set_title(title, fontsize=16)
1491 ax.set_xlabel(
1492 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1493 )
1494 ax.set_ylabel("Normalised probability density", fontsize=14)
1495 if "surface_microphysical" in title and "amount" in title: 1495 ↛ 1496line 1495 didn't jump to line 1496 because the condition on line 1495 was never true
1496 ax.set_ylabel(
1497 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1498 )
1499 ax.set_xlim(vmin, vmax)
1500 ax.tick_params(axis="both", labelsize=12)
1502 # Overlay grid-lines onto histogram plot.
1503 ax.grid(linestyle="--", color="grey", linewidth=1)
1504 if model_colors_map: 1504 ↛ 1505line 1504 didn't jump to line 1505 because the condition on line 1504 was never true
1505 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1507 # Save plot.
1508 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1509 logging.info("Saved histogram plot to %s", filename)
1510 plt.close(fig)
1513def _plot_and_save_postage_stamp_histogram_series(
1514 cube: iris.cube.Cube,
1515 filename: str,
1516 title: str,
1517 stamp_coordinate: str,
1518 vmin: float,
1519 vmax: float,
1520 **kwargs,
1521):
1522 """Plot and save postage (ensemble members) stamps for a histogram series.
1524 Parameters
1525 ----------
1526 cube: Cube
1527 2 dimensional Cube of the data to plot as histogram.
1528 filename: str
1529 Filename of the plot to write.
1530 title: str
1531 Plot title.
1532 stamp_coordinate: str
1533 Coordinate that becomes different plots.
1534 vmin: float
1535 minimum for pdf x-axis
1536 vmax: float
1537 maximum for pdf x-axis
1538 """
1539 # Use the smallest square grid that will fit the members.
1540 nmember = len(cube.coord(stamp_coordinate).points)
1541 grid_rows = int(math.sqrt(nmember))
1542 grid_size = math.ceil(nmember / grid_rows)
1544 fig = plt.figure(
1545 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1546 )
1547 # Make a subplot for each member.
1548 for member, subplot in zip(
1549 cube.slices_over(stamp_coordinate),
1550 range(1, grid_size * grid_rows + 1),
1551 strict=False,
1552 ):
1553 # Implicit interface is much easier here, due to needing to have the
1554 # cartopy GeoAxes generated.
1555 plt.subplot(grid_rows, grid_size, subplot)
1556 # Reshape cube data into a single array to allow for a single histogram.
1557 # Otherwise we plot xdim histograms stacked.
1558 member_data_1d = (member.data).flatten()
1559 plt.hist(member_data_1d, density=True, stacked=True)
1560 axes = plt.gca()
1561 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1562 axes.set_title(f"{mtitle}")
1563 axes.set_xlim(vmin, vmax)
1565 # Overall figure title.
1566 fig.suptitle(title, fontsize=16)
1568 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1569 logging.info("Saved histogram postage stamp plot to %s", filename)
1570 plt.close(fig)
1573def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1574 cube: iris.cube.Cube,
1575 filename: str,
1576 title: str,
1577 stamp_coordinate: str,
1578 vmin: float,
1579 vmax: float,
1580 **kwargs,
1581):
1582 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1583 ax.set_title(title, fontsize=16)
1584 ax.set_xlim(vmin, vmax)
1585 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1586 ax.set_ylabel("normalised probability density", fontsize=14)
1587 # Loop over all slices along the stamp_coordinate
1588 for member in cube.slices_over(stamp_coordinate):
1589 # Flatten the member data to 1D
1590 member_data_1d = member.data.flatten()
1591 # Plot the histogram using plt.hist
1592 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1593 plt.hist(
1594 member_data_1d,
1595 density=True,
1596 stacked=True,
1597 label=f"{mtitle}",
1598 )
1600 # Add a legend
1601 ax.legend(fontsize=16)
1603 # Save the figure to a file
1604 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1605 logging.info("Saved histogram postage stamp plot to %s", filename)
1607 # Close the figure
1608 plt.close(fig)
1611def _plot_and_save_scattermap_plot(
1612 cube: iris.cube.Cube, filename: str, title: str, projection=None, **kwargs
1613):
1614 """Plot and save a geographical scatter plot.
1616 Parameters
1617 ----------
1618 cube: Cube
1619 1 dimensional Cube of the data points with auxiliary latitude and
1620 longitude coordinates,
1621 filename: str
1622 Filename of the plot to write.
1623 title: str
1624 Plot title.
1625 projection: str
1626 Mapping projection to be used by cartopy.
1627 """
1628 # Setup plot details, size, resolution, etc.
1629 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1630 if projection is not None:
1631 # Apart from the default, the only projection we currently support is
1632 # a stereographic projection over the North Pole.
1633 if projection == "NP_Stereo":
1634 axes = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0.0))
1635 else:
1636 raise ValueError(f"Unknown projection: {projection}")
1637 else:
1638 axes = plt.axes(projection=ccrs.PlateCarree())
1640 # Scatter plot of the field. The marker size is chosen to give
1641 # symbols that decrease in size as the number of observations
1642 # increases, although the fraction of the figure covered by
1643 # symbols increases roughly as N^(1/2), disregarding overlaps,
1644 # and has been selected for the default figure size of (10, 10).
1645 # Should this be changed, the marker size should be adjusted in
1646 # proportion to the area of the figure.
1647 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
1648 klon = None
1649 klat = None
1650 for kc in range(len(cube.aux_coords)):
1651 if cube.aux_coords[kc].standard_name == "latitude":
1652 klat = kc
1653 elif cube.aux_coords[kc].standard_name == "longitude":
1654 klon = kc
1655 scatter_map = iplt.scatter(
1656 cube.aux_coords[klon],
1657 cube.aux_coords[klat],
1658 c=cube.data[:],
1659 s=mrk_size,
1660 cmap="jet",
1661 edgecolors="k",
1662 )
1664 # Add coastlines and borderlines.
1665 try:
1666 axes.coastlines(resolution="10m")
1667 axes.add_feature(cfeature.BORDERS)
1668 except AttributeError:
1669 pass
1671 # Add title.
1672 axes.set_title(title, fontsize=16)
1674 # Add colour bar.
1675 cbar = fig.colorbar(scatter_map)
1676 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=20)
1678 # Save plot.
1679 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1680 logging.info("Saved geographical scatter plot to %s", filename)
1681 plt.close(fig)
1684def _spatial_plot(
1685 method: Literal["contourf", "pcolormesh"],
1686 cube: iris.cube.Cube,
1687 filename: str | None,
1688 sequence_coordinate: str,
1689 stamp_coordinate: str,
1690 overlay_cube: iris.cube.Cube | None = None,
1691 contour_cube: iris.cube.Cube | None = None,
1692 **kwargs,
1693):
1694 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1696 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1697 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1698 is present then postage stamp plots will be produced.
1700 If an overlay_cube and/or contour_cube are specified, multiple variables can
1701 be overplotted on the same figure.
1703 Parameters
1704 ----------
1705 method: "contourf" | "pcolormesh"
1706 The plotting method to use.
1707 cube: Cube
1708 Iris cube of the data to plot. It should have two spatial dimensions,
1709 such as lat and lon, and may also have a another two dimension to be
1710 plotted sequentially and/or as postage stamp plots.
1711 filename: str | None
1712 Name of the plot to write, used as a prefix for plot sequences. If None
1713 uses the recipe name.
1714 sequence_coordinate: str
1715 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1716 This coordinate must exist in the cube.
1717 stamp_coordinate: str
1718 Coordinate about which to plot postage stamp plots. Defaults to
1719 ``"realization"``.
1720 overlay_cube: Cube | None, optional
1721 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1722 contour_cube: Cube | None, optional
1723 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1725 Raises
1726 ------
1727 ValueError
1728 If the cube doesn't have the right dimensions.
1729 TypeError
1730 If the cube isn't a single cube.
1731 """
1732 # Ensure we've got a single cube.
1733 cube = check_single_cube(cube)
1735 # Set title based on recipe metadata or use cube name
1736 recipe_title = get_recipe_metadata().get("title", cube.name())
1738 # Check if there is a valid stamp coordinate in cube dimensions.
1739 if stamp_coordinate == "realization": 1739 ↛ 1744line 1739 didn't jump to line 1744 because the condition on line 1739 was always true
1740 stamp_coordinate = check_stamp_coordinate(cube)
1742 # Make postage stamp plots if stamp_coordinate exists and has more than a
1743 # single point.
1744 plotting_func = _plot_and_save_spatial_plot
1745 try:
1746 if cube.coord(stamp_coordinate).shape[0] > 1:
1747 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1748 except iris.exceptions.CoordinateNotFoundError:
1749 pass
1751 # Produce a geographical scatter plot if the data have a
1752 # dimension called observation or model_obs_error
1753 if any( 1753 ↛ 1757line 1753 didn't jump to line 1757 because the condition on line 1753 was never true
1754 crd.var_name == "station" or crd.var_name == "model_obs_error"
1755 for crd in cube.coords()
1756 ):
1757 plotting_func = _plot_and_save_scattermap_plot
1759 # Must have a sequence coordinate.
1760 try:
1761 cube.coord(sequence_coordinate)
1762 except iris.exceptions.CoordinateNotFoundError as err:
1763 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1765 # Create a plot for each value of the sequence coordinate.
1766 plot_index = []
1767 nplot = np.size(cube.coord(sequence_coordinate).points)
1769 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1770 # Set plot titles and filename
1771 seq_coord = cube_slice.coord(sequence_coordinate)
1772 plot_title, plot_filename = _set_title_and_filename(
1773 seq_coord, nplot, recipe_title, filename
1774 )
1776 # Extract sequence slice for overlay_cube and contour_cube if required.
1777 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1778 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1780 # Do the actual plotting.
1781 plotting_func(
1782 cube_slice,
1783 filename=plot_filename,
1784 stamp_coordinate=stamp_coordinate,
1785 title=plot_title,
1786 method=method,
1787 overlay_cube=overlay_slice,
1788 contour_cube=contour_slice,
1789 **kwargs,
1790 )
1791 plot_index.append(plot_filename)
1793 # Add list of plots to plot metadata.
1794 complete_plot_index = _append_to_plot_index(plot_index)
1796 # Make a page to display the plots.
1797 _make_plot_html_page(complete_plot_index)
1800####################
1801# Public functions #
1802####################
1805def spatial_contour_plot(
1806 cube: iris.cube.Cube,
1807 filename: str = None,
1808 sequence_coordinate: str = "time",
1809 stamp_coordinate: str = "realization",
1810 **kwargs,
1811) -> iris.cube.Cube:
1812 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1814 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1815 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1816 is present then postage stamp plots will be produced.
1818 Parameters
1819 ----------
1820 cube: Cube
1821 Iris cube of the data to plot. It should have two spatial dimensions,
1822 such as lat and lon, and may also have a another two dimension to be
1823 plotted sequentially and/or as postage stamp plots.
1824 filename: str, optional
1825 Name of the plot to write, used as a prefix for plot sequences. Defaults
1826 to the recipe name.
1827 sequence_coordinate: str, optional
1828 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1829 This coordinate must exist in the cube.
1830 stamp_coordinate: str, optional
1831 Coordinate about which to plot postage stamp plots. Defaults to
1832 ``"realization"``.
1834 Returns
1835 -------
1836 Cube
1837 The original cube (so further operations can be applied).
1839 Raises
1840 ------
1841 ValueError
1842 If the cube doesn't have the right dimensions.
1843 TypeError
1844 If the cube isn't a single cube.
1845 """
1846 _spatial_plot(
1847 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1848 )
1849 return cube
1852def spatial_pcolormesh_plot(
1853 cube: iris.cube.Cube,
1854 filename: str = None,
1855 sequence_coordinate: str = "time",
1856 stamp_coordinate: str = "realization",
1857 **kwargs,
1858) -> iris.cube.Cube:
1859 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1861 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1862 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1863 is present then postage stamp plots will be produced.
1865 This function is significantly faster than ``spatial_contour_plot``,
1866 especially at high resolutions, and should be preferred unless contiguous
1867 contour areas are important.
1869 Parameters
1870 ----------
1871 cube: Cube
1872 Iris cube of the data to plot. It should have two spatial dimensions,
1873 such as lat and lon, and may also have a another two dimension to be
1874 plotted sequentially and/or as postage stamp plots.
1875 filename: str, optional
1876 Name of the plot to write, used as a prefix for plot sequences. Defaults
1877 to the recipe name.
1878 sequence_coordinate: str, optional
1879 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1880 This coordinate must exist in the cube.
1881 stamp_coordinate: str, optional
1882 Coordinate about which to plot postage stamp plots. Defaults to
1883 ``"realization"``.
1885 Returns
1886 -------
1887 Cube
1888 The original cube (so further operations can be applied).
1890 Raises
1891 ------
1892 ValueError
1893 If the cube doesn't have the right dimensions.
1894 TypeError
1895 If the cube isn't a single cube.
1896 """
1897 _spatial_plot(
1898 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1899 )
1900 return cube
1903def spatial_multi_pcolormesh_plot(
1904 cube: iris.cube.Cube,
1905 overlay_cube: iris.cube.Cube | None = None,
1906 contour_cube: iris.cube.Cube | None = None,
1907 filename: str = None,
1908 sequence_coordinate: str = "time",
1909 stamp_coordinate: str = "realization",
1910 **kwargs,
1911) -> iris.cube.Cube:
1912 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1914 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1915 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1916 is present then postage stamp plots will be produced.
1918 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1920 If specified, contours of a contour_cube can be overplotted on top of those.
1922 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1924 This function is significantly faster than ``spatial_contour_plot``,
1925 especially at high resolutions, and should be preferred unless contiguous
1926 contour areas are important.
1928 Parameters
1929 ----------
1930 cube: Cube
1931 Iris cube of the data to plot. It should have two spatial dimensions,
1932 such as lat and lon, and may also have a another two dimension to be
1933 plotted sequentially and/or as postage stamp plots.
1934 overlay_cube: Cube, optional
1935 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1936 such as lat and lon, and may also have a another two dimension to be
1937 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.
1938 If not provided, output plot generated without overlay cube.
1939 contour_cube: Cube, optional
1940 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,
1941 such as lat and lon, and may also have a another two dimension to be
1942 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
1943 filename: str, optional
1944 Name of the plot to write, used as a prefix for plot sequences. Defaults
1945 to the recipe name.
1946 sequence_coordinate: str, optional
1947 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1948 This coordinate must exist in the cube.
1949 stamp_coordinate: str, optional
1950 Coordinate about which to plot postage stamp plots. Defaults to
1951 ``"realization"``.
1953 Returns
1954 -------
1955 Cube
1956 The original cube (so further operations can be applied).
1958 Raises
1959 ------
1960 ValueError
1961 If the cube doesn't have the right dimensions.
1962 TypeError
1963 If the cube isn't a single cube.
1964 """
1965 _spatial_plot(
1966 "pcolormesh",
1967 cube,
1968 filename,
1969 sequence_coordinate,
1970 stamp_coordinate,
1971 overlay_cube=overlay_cube,
1972 contour_cube=contour_cube,
1973 )
1974 return cube, overlay_cube, contour_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: 2854 ↛ 2867line 2854 didn't jump to line 2867 because the condition on line 2854 was always true
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): 2878 ↛ 2879line 2878 didn't jump to line 2879 because the condition on line 2878 was never true
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)