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