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