Coverage for src/CSET/operators/plot.py: 81%
1011 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-20 13:56 +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.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(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 ↛ 758line 755 didn't jump to line 758 because the condition on line 755 was always true
756 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
757 # Tick labels for rainfall rates from Nimrod radar data.
758 if "rainfall rate composite" in cube.name(): 758 ↛ 759line 758 didn't jump to line 759 because the condition on line 758 was never true
759 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
760 # Tick labels for rain accumulations from Nimrod radar data.
761 if "rain accumulation" in cube.name(): 761 ↛ 762line 761 didn't jump to line 762 because the condition on line 761 was never true
762 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
763 if "wts accumulation" in cube.name(): 763 ↛ 764line 763 didn't jump to line 764 because the condition on line 763 was never true
764 tick_levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]
765 cbar.minorticks_off()
766 cbar.set_ticks(tick_levels)
767 cbar.set_ticklabels([f"{level:.3g}" for level in tick_levels])
768 cbar.set_label(label=f"{cube.name()}", size=14)
769 # Tick labels for model rainfall data.
770 if "surface_microphysical" in cube.name(): 770 ↛ 773line 770 didn't jump to line 773 because the condition on line 770 was always true
771 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
772 # Tick labels for Nimrod weights data.
773 logging.debug("Set colorbar ticks and labels.")
775 # Save plot.
776 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
777 logging.info("Saved spatial plot to %s", filename)
778 plt.close(fig)
781def _plot_and_save_postage_stamp_spatial_plot(
782 cube: iris.cube.Cube,
783 filename: str,
784 stamp_coordinate: str,
785 title: str,
786 method: Literal["contourf", "pcolormesh"],
787 overlay_cube: iris.cube.Cube | None = None,
788 contour_cube: iris.cube.Cube | None = None,
789 **kwargs,
790):
791 """Plot postage stamp spatial plots from an ensemble.
793 Parameters
794 ----------
795 cube: Cube
796 Iris cube of data to be plotted. It must have the stamp coordinate.
797 filename: str
798 Filename of the plot to write.
799 stamp_coordinate: str
800 Coordinate that becomes different plots.
801 method: "contourf" | "pcolormesh"
802 The plotting method to use.
803 overlay_cube: Cube, optional
804 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
805 contour_cube: Cube, optional
806 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
808 Raises
809 ------
810 ValueError
811 If the cube doesn't have the right dimensions.
812 """
813 # Use the smallest square grid that will fit the members.
814 nmember = len(cube.coord(stamp_coordinate).points)
815 grid_rows = int(math.sqrt(nmember))
816 grid_size = math.ceil(nmember / grid_rows)
818 fig = plt.figure(
819 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
820 )
822 # Specify the color bar
823 cmap, levels, norm = colorbar_map_levels(cube)
824 # If overplotting, set required colorbars
825 if overlay_cube: 825 ↛ 826line 825 didn't jump to line 826 because the condition on line 825 was never true
826 over_cmap, over_levels, over_norm = colorbar_map_levels(overlay_cube)
827 if contour_cube: 827 ↛ 828line 827 didn't jump to line 828 because the condition on line 827 was never true
828 cntr_cmap, cntr_levels, cntr_norm = colorbar_map_levels(contour_cube)
830 # Make a subplot for each member.
831 for member, subplot in zip(
832 cube.slices_over(stamp_coordinate),
833 range(1, grid_size * grid_rows + 1),
834 strict=False,
835 ):
836 # Setup subplot map projection, extent and coastlines and borderlines.
837 axes = _setup_spatial_map(
838 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
839 )
840 if method == "contourf":
841 # Filled contour plot of the field.
842 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
843 elif method == "pcolormesh":
844 if levels is not None:
845 vmin = min(levels)
846 vmax = max(levels)
847 else:
848 raise TypeError("Unknown vmin and vmax range.")
849 vmin, vmax = None, None
850 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
851 # if levels are defined.
852 if norm is not None: 852 ↛ 853line 852 didn't jump to line 853 because the condition on line 852 was never true
853 vmin = None
854 vmax = None
855 # pcolormesh plot of the field.
856 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
857 else:
858 raise ValueError(f"Unknown plotting method: {method}")
860 # Overplot overlay field, if required
861 if overlay_cube: 861 ↛ 862line 861 didn't jump to line 862 because the condition on line 861 was never true
862 try:
863 over_vmin = min(over_levels)
864 over_vmax = max(over_levels)
865 except TypeError:
866 over_vmin, over_vmax = None, None
867 if over_norm is not None:
868 over_vmin = None
869 over_vmax = None
870 iplt.pcolormesh(
871 overlay_cube[member.coord(stamp_coordinate).points[0]],
872 cmap=over_cmap,
873 norm=over_norm,
874 alpha=0.6,
875 vmin=over_vmin,
876 vmax=over_vmax,
877 )
878 # Overplot contour field, if required
879 if contour_cube: 879 ↛ 880line 879 didn't jump to line 880 because the condition on line 879 was never true
880 iplt.contour(
881 contour_cube[member.coord(stamp_coordinate).points[0]],
882 colors="darkgray",
883 levels=cntr_levels,
884 norm=cntr_norm,
885 alpha=0.6,
886 linestyles="--",
887 linewidths=1,
888 )
889 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
890 axes.set_title(f"{mtitle}")
892 # Put the shared colorbar in its own axes.
893 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
894 colorbar = fig.colorbar(
895 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
896 )
897 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
899 # Overall figure title.
900 fig.suptitle(title, fontsize=16)
902 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
903 logging.info("Saved contour postage stamp plot to %s", filename)
904 plt.close(fig)
907def _plot_and_save_line_series(
908 cubes: iris.cube.CubeList,
909 coords: list[iris.coords.Coord],
910 ensemble_coord: str,
911 filename: str,
912 title: str,
913 **kwargs,
914):
915 """Plot and save a 1D line series.
917 Parameters
918 ----------
919 cubes: Cube or CubeList
920 Cube or CubeList containing the cubes to plot on the y-axis.
921 coords: list[Coord]
922 Coordinates to plot on the x-axis, one per cube.
923 ensemble_coord: str
924 Ensemble coordinate in the cube.
925 filename: str
926 Filename of the plot to write.
927 title: str
928 Plot title.
929 """
930 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
932 model_colors_map = get_model_colors_map(cubes)
934 # Store min/max ranges.
935 y_levels = []
937 # Check match-up across sequence coords gives consistent sizes
938 validate_cubes_coords(cubes, coords)
940 for cube, coord in zip(cubes, coords, strict=True):
941 label = None
942 color = "black"
943 if model_colors_map:
944 label = cube.attributes.get("model_name")
945 color = model_colors_map.get(label)
946 for cube_slice in cube.slices_over(ensemble_coord):
947 # Label with (control) if part of an ensemble or not otherwise.
948 if cube_slice.coord(ensemble_coord).points == [0]:
949 iplt.plot(
950 coord,
951 cube_slice,
952 color=color,
953 marker="o",
954 ls="-",
955 lw=3,
956 label=f"{label} (control)"
957 if len(cube.coord(ensemble_coord).points) > 1
958 else label,
959 )
960 # Label with (perturbed) if part of an ensemble and not the control.
961 else:
962 iplt.plot(
963 coord,
964 cube_slice,
965 color=color,
966 ls="-",
967 lw=1.5,
968 alpha=0.75,
969 label=f"{label} (member)",
970 )
972 # Calculate the global min/max if multiple cubes are given.
973 _, levels, _ = colorbar_map_levels(cube, axis="y")
974 if levels is not None: 974 ↛ 975line 974 didn't jump to line 975 because the condition on line 974 was never true
975 y_levels.append(min(levels))
976 y_levels.append(max(levels))
978 # Get the current axes.
979 ax = plt.gca()
981 # Add some labels and tweak the style.
982 # check if cubes[0] works for single cube if not CubeList
983 if coords[0].name() == "time":
984 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
985 else:
986 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
987 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
988 ax.set_title(title, fontsize=16)
990 ax.ticklabel_format(axis="y", useOffset=False)
991 ax.tick_params(axis="x", labelrotation=15)
992 ax.tick_params(axis="both", labelsize=12)
994 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
995 if y_levels: 995 ↛ 996line 995 didn't jump to line 996 because the condition on line 995 was never true
996 ax.set_ylim(min(y_levels), max(y_levels))
997 # Add zero line.
998 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
999 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
1000 logging.debug(
1001 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
1002 )
1003 else:
1004 ax.autoscale()
1006 # Add gridlines
1007 ax.grid(linestyle="--", color="grey", linewidth=1)
1008 # Ientify unique labels for legend
1009 handles = list(
1010 {
1011 label: handle
1012 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1013 }.values()
1014 )
1015 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1017 # Save plot.
1018 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1019 logging.info("Saved line plot to %s", filename)
1020 plt.close(fig)
1023def _plot_and_save_line_power_spectrum_series(
1024 cubes: iris.cube.Cube | iris.cube.CubeList,
1025 coords: list[iris.coords.Coord],
1026 ensemble_coord: str,
1027 filename: str,
1028 title: str,
1029 series_coordinate: str,
1030 **kwargs,
1031):
1032 """Plot and save a 1D line series.
1034 Parameters
1035 ----------
1036 cubes: Cube or CubeList
1037 Cube or CubeList containing the cubes to plot on the y-axis.
1038 coords: list[Coord]
1039 Coordinates to plot on the x-axis, one per cube.
1040 ensemble_coord: str
1041 Ensemble coordinate in the cube.
1042 filename: str
1043 Filename of the plot to write.
1044 title: str
1045 Plot title.
1046 series_coordinate: str
1047 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
1048 """
1049 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1050 model_colors_map = get_model_colors_map(cubes)
1051 ax = plt.gca()
1053 # Store min/max ranges.
1054 y_levels = []
1056 line_marker = None
1057 line_width = 1
1059 for cube in iter_maybe(cubes):
1060 # next 2 lines replace chunk of code.
1061 xcoord = _select_series_coord(cube, series_coordinate)
1062 xname = xcoord.points
1064 yfield = cube.data # power spectrum
1065 label = None
1066 color = "black"
1067 if model_colors_map: 1067 ↛ 1070line 1067 didn't jump to line 1070 because the condition on line 1067 was always true
1068 label = cube.attributes.get("model_name")
1069 color = model_colors_map.get(label)
1070 for cube_slice in cube.slices_over(ensemble_coord):
1071 # Label with (control) if part of an ensemble or not otherwise.
1072 if cube_slice.coord(ensemble_coord).points == [0]: 1072 ↛ 1086line 1072 didn't jump to line 1086 because the condition on line 1072 was always true
1073 ax.plot(
1074 xname,
1075 yfield,
1076 color=color,
1077 marker=line_marker,
1078 ls="-",
1079 lw=line_width,
1080 label=f"{label} (control)"
1081 if len(cube.coord(ensemble_coord).points) > 1
1082 else label,
1083 )
1084 # Label with (perturbed) if part of an ensemble and not the control.
1085 else:
1086 ax.plot(
1087 xname,
1088 yfield,
1089 color=color,
1090 ls="-",
1091 lw=1.5,
1092 alpha=0.75,
1093 label=f"{label} (member)",
1094 )
1096 # Calculate the global min/max if multiple cubes are given.
1097 _, levels, _ = colorbar_map_levels(cube, axis="y")
1098 if levels is not None: 1098 ↛ 1099line 1098 didn't jump to line 1099 because the condition on line 1098 was never true
1099 y_levels.append(min(levels))
1100 y_levels.append(max(levels))
1102 # Add some labels and tweak the style.
1104 title = f"{title}"
1105 ax.set_title(title, fontsize=16)
1107 # Set appropriate x-axis label based on coordinate
1108 if series_coordinate == "wavelength" or ( 1108 ↛ 1111line 1108 didn't jump to line 1111 because the condition on line 1108 was never true
1109 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
1110 ):
1111 ax.set_xlabel("Wavelength (km)", fontsize=14)
1112 elif series_coordinate == "physical_wavenumber" or ( 1112 ↛ 1115line 1112 didn't jump to line 1115 because the condition on line 1112 was never true
1113 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
1114 ):
1115 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1116 else: # frequency or check units
1117 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1": 1117 ↛ 1118line 1117 didn't jump to line 1118 because the condition on line 1117 was never true
1118 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
1119 else:
1120 ax.set_xlabel("Wavenumber", fontsize=14)
1122 ax.set_ylabel("Power Spectral Density", fontsize=14)
1123 ax.tick_params(axis="both", labelsize=12)
1125 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
1127 # Set log-log scale
1128 ax.set_xscale("log")
1129 ax.set_yscale("log")
1131 # Add gridlines
1132 ax.grid(linestyle="--", color="grey", linewidth=1)
1133 # Ientify unique labels for legend
1134 handles = list(
1135 {
1136 label: handle
1137 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1138 }.values()
1139 )
1140 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1142 # Save plot.
1143 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1144 logging.info("Saved line plot to %s", filename)
1145 plt.close(fig)
1148def _plot_and_save_vertical_line_series(
1149 cubes: iris.cube.CubeList,
1150 coords: list[iris.coords.Coord],
1151 ensemble_coord: str,
1152 filename: str,
1153 series_coordinate: str,
1154 title: str,
1155 vmin: float,
1156 vmax: float,
1157 **kwargs,
1158):
1159 """Plot and save a 1D line series in vertical.
1161 Parameters
1162 ----------
1163 cubes: CubeList
1164 1 dimensional Cube or CubeList of the data to plot on x-axis.
1165 coord: list[Coord]
1166 Coordinates to plot on the y-axis, one per cube.
1167 ensemble_coord: str
1168 Ensemble coordinate in the cube.
1169 filename: str
1170 Filename of the plot to write.
1171 series_coordinate: str
1172 Coordinate to use as vertical axis.
1173 title: str
1174 Plot title.
1175 vmin: float
1176 Minimum value for the x-axis.
1177 vmax: float
1178 Maximum value for the x-axis.
1179 """
1180 # plot the vertical pressure axis using log scale
1181 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1183 model_colors_map = get_model_colors_map(cubes)
1185 # Check match-up across sequence coords gives consistent sizes
1186 validate_cubes_coords(cubes, coords)
1188 for cube, coord in zip(cubes, coords, strict=True):
1189 label = None
1190 color = "black"
1191 if model_colors_map: 1191 ↛ 1192line 1191 didn't jump to line 1192 because the condition on line 1191 was never true
1192 label = cube.attributes.get("model_name")
1193 color = model_colors_map.get(label)
1195 for cube_slice in cube.slices_over(ensemble_coord):
1196 # If ensemble data given plot control member with (control)
1197 # unless single forecast.
1198 if cube_slice.coord(ensemble_coord).points == [0]:
1199 iplt.plot(
1200 cube_slice,
1201 coord,
1202 color=color,
1203 marker="o",
1204 ls="-",
1205 lw=3,
1206 label=f"{label} (control)"
1207 if len(cube.coord(ensemble_coord).points) > 1
1208 else label,
1209 )
1210 # If ensemble data given plot perturbed members with (perturbed).
1211 else:
1212 iplt.plot(
1213 cube_slice,
1214 coord,
1215 color=color,
1216 ls="-",
1217 lw=1.5,
1218 alpha=0.75,
1219 label=f"{label} (member)",
1220 )
1222 # Get the current axis
1223 ax = plt.gca()
1225 # Special handling for pressure level data.
1226 if series_coordinate == "pressure": 1226 ↛ 1248line 1226 didn't jump to line 1248 because the condition on line 1226 was always true
1227 # Invert y-axis and set to log scale.
1228 ax.invert_yaxis()
1229 ax.set_yscale("log")
1231 # Define y-ticks and labels for pressure log axis.
1232 y_tick_labels = [
1233 "1000",
1234 "850",
1235 "700",
1236 "500",
1237 "300",
1238 "200",
1239 "100",
1240 ]
1241 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1243 # Set y-axis limits and ticks.
1244 ax.set_ylim(1100, 100)
1246 # Test if series_coordinate is model level data. The UM data uses
1247 # model_level_number and lfric uses full_levels as coordinate.
1248 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1249 # Define y-ticks and labels for vertical axis.
1250 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1251 y_tick_labels = [str(int(i)) for i in y_ticks]
1252 ax.set_ylim(min(y_ticks), max(y_ticks))
1254 ax.set_yticks(y_ticks)
1255 ax.set_yticklabels(y_tick_labels)
1257 # Set x-axis limits.
1258 ax.set_xlim(vmin, vmax)
1259 # Mark y=0 if present in plot.
1260 if vmin < 0.0 and vmax > 0.0: 1260 ↛ 1261line 1260 didn't jump to line 1261 because the condition on line 1260 was never true
1261 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1263 # Add some labels and tweak the style.
1264 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1265 ax.set_xlabel(
1266 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1267 )
1268 ax.set_title(title, fontsize=16)
1269 ax.ticklabel_format(axis="x")
1270 ax.tick_params(axis="y")
1271 ax.tick_params(axis="both", labelsize=12)
1273 # Add gridlines
1274 ax.grid(linestyle="--", color="grey", linewidth=1)
1275 # Ientify unique labels for legend
1276 handles = list(
1277 {
1278 label: handle
1279 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1280 }.values()
1281 )
1282 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1284 # Save plot.
1285 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1286 logging.info("Saved line plot to %s", filename)
1287 plt.close(fig)
1290def _plot_and_save_scatter_plot(
1291 cube_x: iris.cube.Cube | iris.cube.CubeList,
1292 cube_y: iris.cube.Cube | iris.cube.CubeList,
1293 filename: str,
1294 title: str,
1295 one_to_one: bool,
1296 model_names: list[str] = None,
1297 **kwargs,
1298):
1299 """Plot and save a 2D scatter plot.
1301 Parameters
1302 ----------
1303 cube_x: Cube | CubeList
1304 1 dimensional Cube or CubeList of the data to plot on x-axis.
1305 cube_y: Cube | CubeList
1306 1 dimensional Cube or CubeList of the data to plot on y-axis.
1307 filename: str
1308 Filename of the plot to write.
1309 title: str
1310 Plot title.
1311 one_to_one: bool
1312 Whether a 1:1 line is plotted.
1313 """
1314 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1315 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1316 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1317 # over the pairs simultaneously.
1319 # Ensure cube_x and cube_y are iterable
1320 cube_x_iterable = iter_maybe(cube_x)
1321 cube_y_iterable = iter_maybe(cube_y)
1323 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1324 iplt.scatter(cube_x_iter, cube_y_iter)
1325 if one_to_one is True:
1326 plt.plot(
1327 [
1328 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1329 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1330 ],
1331 [
1332 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1333 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1334 ],
1335 "k",
1336 linestyle="--",
1337 )
1338 ax = plt.gca()
1340 # Add some labels and tweak the style.
1341 if model_names is None:
1342 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1343 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1344 else:
1345 # Add the model names, these should be order of base (x) and other (y).
1346 ax.set_xlabel(
1347 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1348 )
1349 ax.set_ylabel(
1350 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1351 )
1352 ax.set_title(title, fontsize=16)
1353 ax.ticklabel_format(axis="y", useOffset=False)
1354 ax.tick_params(axis="x", labelrotation=15)
1355 ax.tick_params(axis="both", labelsize=12)
1356 ax.autoscale()
1358 # Save plot.
1359 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1360 logging.info("Saved scatter plot to %s", filename)
1361 plt.close(fig)
1364def _plot_and_save_vector_plot(
1365 cube_u: iris.cube.Cube,
1366 cube_v: iris.cube.Cube,
1367 filename: str,
1368 title: str,
1369 method: Literal["contourf", "pcolormesh"],
1370 **kwargs,
1371):
1372 """Plot and save a 2D vector plot.
1374 Parameters
1375 ----------
1376 cube_u: Cube
1377 2 dimensional Cube of u component of the data.
1378 cube_v: Cube
1379 2 dimensional Cube of v component of the data.
1380 filename: str
1381 Filename of the plot to write.
1382 title: str
1383 Plot title.
1384 """
1385 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1386 # Create a cube containing the magnitude of the vector field.
1387 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1388 cube_vec_mag.rename(f"{cube_u.long_name}_{cube_v.long_name}_magnitude")
1389 if "eastward_wind" in cube_u.long_name and "northward_wind" in cube_v.long_name:
1390 cube_vec_mag.rename(
1391 "wind_speed" + cube_u.long_name.replace("eastward_wind", "")
1392 )
1394 # Specify the color bar
1395 cmap, levels, norm = colorbar_map_levels(cube_vec_mag)
1397 # Setup plot map projection, extent and coastlines and borderlines.
1398 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1400 if method == "contourf":
1401 # Filled contour plot of the field.
1402 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1403 elif method == "pcolormesh":
1404 try:
1405 vmin = min(levels)
1406 vmax = max(levels)
1407 except TypeError:
1408 vmin, vmax = None, None
1409 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1410 # if levels are defined.
1411 if norm is not None:
1412 vmin = None
1413 vmax = None
1414 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1415 else:
1416 raise ValueError(f"Unknown plotting method: {method}")
1418 # Check to see if transect, and if so, adjust y axis.
1419 if is_transect(cube_vec_mag):
1420 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1421 axes.invert_yaxis()
1422 axes.set_yscale("log")
1423 axes.set_ylim(1100, 100)
1424 # If both model_level_number and level_height exists, iplt can construct
1425 # plot as a function of height above orography (NOT sea level).
1426 elif {"model_level_number", "level_height"}.issubset(
1427 {coord.name() for coord in cube_vec_mag.coords()}
1428 ):
1429 axes.set_yscale("log")
1431 axes.set_title(
1432 f"{title}\n"
1433 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1434 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1435 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1436 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1437 fontsize=16,
1438 )
1440 else:
1441 # Add title.
1442 axes.set_title(title, fontsize=16)
1444 # Add watermark with min/max/mean. Currently not user togglable.
1445 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1446 axes.annotate(
1447 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}",
1448 xy=(0.05, -0.05),
1449 xycoords="axes fraction",
1450 xytext=(-5, 5),
1451 textcoords="offset points",
1452 ha="right",
1453 va="bottom",
1454 size=11,
1455 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1456 )
1458 # Add colour bar.
1459 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1460 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1461 # add ticks and tick_labels for every levels if less than 20 levels exist
1462 if levels is not None and len(levels) < 20:
1463 cbar.set_ticks(levels)
1464 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1466 # 30 barbs along the longest axis of the plot, or a barb per point for data
1467 # with less than 30 points.
1468 step = max(max(cube_u.shape) // 30, 1)
1469 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1471 # Save plot.
1472 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1473 logging.info("Saved vector plot to %s", filename)
1474 plt.close(fig)
1477def _plot_and_save_histogram_series(
1478 cubes: iris.cube.Cube | iris.cube.CubeList,
1479 filename: str,
1480 title: str,
1481 vmin: float,
1482 vmax: float,
1483 **kwargs,
1484):
1485 """Plot and save a histogram series.
1487 Parameters
1488 ----------
1489 cubes: Cube or CubeList
1490 2 dimensional Cube or CubeList of the data to plot as histogram.
1491 filename: str
1492 Filename of the plot to write.
1493 title: str
1494 Plot title.
1495 vmin: float
1496 minimum for colorbar
1497 vmax: float
1498 maximum for colorbar
1499 """
1500 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1501 ax = plt.gca()
1503 model_colors_map = get_model_colors_map(cubes)
1505 # Set default that histograms will produce probability density function
1506 # at each bin (integral over range sums to 1).
1507 density = True
1509 for cube in iter_maybe(cubes):
1510 # Easier to check title (where var name originates)
1511 # than seeing if long names exist etc.
1512 # Exception case, where distribution better fits log scales/bins.
1513 if (
1514 ("surface_microphysical" in title)
1515 or ("rain accumulation" in title)
1516 or ("Rainfall rate Composite" in title)
1517 or ("Nimrod_5min" in title)
1518 ):
1519 if "amount" in title: 1519 ↛ 1521line 1519 didn't jump to line 1521 because the condition on line 1519 was never true
1520 # Compute histogram following Klingaman et al. (2017): ASoP
1521 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1522 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1523 density = False
1524 else:
1525 bins = 10.0 ** (
1526 np.arange(-10, 27, 1) / 10.0
1527 ) # Suggestion from RMED toolbox.
1528 bins = np.insert(bins, 0, 0)
1529 ax.set_yscale("log")
1530 vmin = bins[1]
1531 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1532 ax.set_xscale("log")
1533 elif "lightning" in title:
1534 bins = [0, 1, 2, 3, 4, 5]
1535 else:
1536 bins = np.linspace(vmin, vmax, 51)
1537 logging.debug(
1538 "Plotting histogram with %s bins %s - %s.",
1539 np.size(bins),
1540 np.min(bins),
1541 np.max(bins),
1542 )
1544 # Reshape cube data into a single array to allow for a single histogram.
1545 # Otherwise we plot xdim histograms stacked.
1546 cube_data_1d = (cube.data).flatten()
1548 label = None
1549 color = "black"
1550 if model_colors_map:
1551 label = cube.attributes.get("model_name")
1552 color = model_colors_map[label]
1553 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1555 # Compute area under curve.
1556 if ( 1556 ↛ 1562line 1556 didn't jump to line 1562 because the condition on line 1556 was never true
1557 ("surface_microphysical" in title and "amount" in title)
1558 or ("rain_accumulation" in title)
1559 or ("Rainfall rate Composite" in title)
1560 or ("Nimrod_5min" in title)
1561 ):
1562 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1563 x = x * bin_mean / x.sum()
1564 x = x[1:]
1565 y = y[1:]
1567 ax.plot(
1568 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1569 )
1571 # Add some labels and tweak the style.
1572 ax.set_title(title, fontsize=16)
1573 ax.set_xlabel(
1574 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1575 )
1576 ax.set_ylabel("Normalised probability density", fontsize=14)
1577 if ( 1577 ↛ 1582line 1577 didn't jump to line 1582 because the condition on line 1577 was never true
1578 ("surface_microphysical" in title and "amount" in title)
1579 or ("rain accumulation" in title)
1580 or ("Nimrod_5min" in title)
1581 ):
1582 ax.set_ylabel(
1583 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1584 )
1585 ax.set_xlim(vmin, vmax)
1586 ax.tick_params(axis="both", labelsize=12)
1588 # Overlay grid-lines onto histogram plot.
1589 ax.grid(linestyle="--", color="grey", linewidth=1)
1590 if model_colors_map:
1591 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1593 # Save plot.
1594 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1595 logging.info("Saved histogram plot to %s", filename)
1596 plt.close(fig)
1599def _plot_and_save_postage_stamp_histogram_series(
1600 cube: iris.cube.Cube,
1601 filename: str,
1602 title: str,
1603 stamp_coordinate: str,
1604 vmin: float,
1605 vmax: float,
1606 **kwargs,
1607):
1608 """Plot and save postage (ensemble members) stamps for a histogram series.
1610 Parameters
1611 ----------
1612 cube: Cube
1613 2 dimensional Cube of the data to plot as histogram.
1614 filename: str
1615 Filename of the plot to write.
1616 title: str
1617 Plot title.
1618 stamp_coordinate: str
1619 Coordinate that becomes different plots.
1620 vmin: float
1621 minimum for pdf x-axis
1622 vmax: float
1623 maximum for pdf x-axis
1624 """
1625 # Use the smallest square grid that will fit the members.
1626 nmember = len(cube.coord(stamp_coordinate).points)
1627 grid_rows = int(math.sqrt(nmember))
1628 grid_size = math.ceil(nmember / grid_rows)
1630 fig = plt.figure(
1631 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1632 )
1633 # Make a subplot for each member.
1634 for member, subplot in zip(
1635 cube.slices_over(stamp_coordinate),
1636 range(1, grid_size * grid_rows + 1),
1637 strict=False,
1638 ):
1639 # Implicit interface is much easier here, due to needing to have the
1640 # cartopy GeoAxes generated.
1641 plt.subplot(grid_rows, grid_size, subplot)
1642 # Reshape cube data into a single array to allow for a single histogram.
1643 # Otherwise we plot xdim histograms stacked.
1644 member_data_1d = (member.data).flatten()
1645 plt.hist(member_data_1d, density=True, stacked=True)
1646 axes = plt.gca()
1647 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1648 axes.set_title(f"{mtitle}")
1649 axes.set_xlim(vmin, vmax)
1651 # Overall figure title.
1652 fig.suptitle(title, fontsize=16)
1654 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1655 logging.info("Saved histogram postage stamp plot to %s", filename)
1656 plt.close(fig)
1659def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1660 cube: iris.cube.Cube,
1661 filename: str,
1662 title: str,
1663 stamp_coordinate: str,
1664 vmin: float,
1665 vmax: float,
1666 **kwargs,
1667):
1668 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1669 ax.set_title(title, fontsize=16)
1670 ax.set_xlim(vmin, vmax)
1671 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1672 ax.set_ylabel("normalised probability density", fontsize=14)
1673 # Loop over all slices along the stamp_coordinate
1674 for member in cube.slices_over(stamp_coordinate):
1675 # Flatten the member data to 1D
1676 member_data_1d = member.data.flatten()
1677 # Plot the histogram using plt.hist
1678 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1679 plt.hist(
1680 member_data_1d,
1681 density=True,
1682 stacked=True,
1683 label=f"{mtitle}",
1684 )
1686 # Add a legend
1687 ax.legend(fontsize=16)
1689 # Save the figure to a file
1690 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1691 logging.info("Saved histogram postage stamp plot to %s", filename)
1693 # Close the figure
1694 plt.close(fig)
1697def _spatial_plot(
1698 method: Literal["contourf", "pcolormesh", "scatter"],
1699 cube: iris.cube.Cube,
1700 filename: str | None,
1701 sequence_coordinate: str,
1702 stamp_coordinate: str,
1703 overlay_cube: iris.cube.Cube | None = None,
1704 contour_cube: iris.cube.Cube | None = None,
1705 point_cube: iris.cube.Cube | None = None,
1706 **kwargs,
1707):
1708 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1710 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1711 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1712 is present then postage stamp plots will be produced.
1714 If any optional overlay_cube, contour_cube or point_cube are specified, multiple data layers can
1715 be overplotted on the same figure.
1717 Parameters
1718 ----------
1719 method: "contourf" | "pcolormesh" | "scatter"
1720 The plotting method to use.
1721 Select choice of "contourf" or "pcolormesh" for gridded data.
1722 Use "scatter" for point-based data.
1723 cube: Cube
1724 Iris cube of the data to plot. It should have two spatial dimensions,
1725 such as lat and lon, and may also have a another two dimension to be
1726 plotted sequentially and/or as postage stamp plots.
1727 filename: str | None
1728 Name of the plot to write, used as a prefix for plot sequences. If None
1729 uses the recipe name.
1730 sequence_coordinate: str
1731 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1732 This coordinate must exist in the cube.
1733 stamp_coordinate: str
1734 Coordinate about which to plot postage stamp plots. Defaults to
1735 ``"realization"``.
1736 overlay_cube: Cube | None, optional
1737 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1738 contour_cube: Cube | None, optional
1739 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1740 point_cube: Cube | None, optional
1741 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
1743 Raises
1744 ------
1745 ValueError
1746 If the cube doesn't have the right dimensions.
1747 TypeError
1748 If the cube isn't a single cube.
1749 """
1750 # Ensure we've got a single cube.
1751 cube = check_single_cube(cube)
1753 # Set title based on recipe metadata or use cube name
1754 recipe_title = get_recipe_metadata().get("title", cube.name())
1756 # Check if there is a valid stamp coordinate in cube dimensions.
1757 if stamp_coordinate == "realization": 1757 ↛ 1762line 1757 didn't jump to line 1762 because the condition on line 1757 was always true
1758 stamp_coordinate = check_stamp_coordinate(cube)
1760 # Make postage stamp plots if stamp_coordinate exists and has more than a
1761 # single point.
1762 plotting_func = _plot_and_save_spatial_plot
1763 try:
1764 if cube.coord(stamp_coordinate).shape[0] > 1:
1765 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1766 except iris.exceptions.CoordinateNotFoundError:
1767 pass
1769 # Produce a geographical scatter plot if the data have a
1770 # dimension called observation or model_obs_error
1771 if any( 1771 ↛ 1775line 1771 didn't jump to line 1775 because the condition on line 1771 was never true
1772 crd.var_name == "station" or crd.var_name == "model_obs_error"
1773 for crd in cube.coords()
1774 ):
1775 plotting_func = _plot_and_save_spatial_plot
1776 method = "scatter"
1778 # Must have a sequence coordinate.
1779 try:
1780 cube.coord(sequence_coordinate)
1781 except iris.exceptions.CoordinateNotFoundError as err:
1782 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1784 # Create a plot for each value of the sequence coordinate.
1785 plot_index = []
1786 nplot = np.size(cube.coord(sequence_coordinate).points)
1788 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1789 # Set plot titles and filename
1790 seq_coord = cube_slice.coord(sequence_coordinate)
1791 plot_title, plot_filename = _set_title_and_filename(
1792 seq_coord, nplot, recipe_title, filename
1793 )
1795 # Extract sequence slice for overlay_cube, contour_cube and point_cube if required.
1796 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1797 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1798 point_slice = slice_over_maybe(point_cube, sequence_coordinate, iseq)
1800 # Do the actual plotting.
1801 plotting_func(
1802 cube_slice,
1803 filename=plot_filename,
1804 stamp_coordinate=stamp_coordinate,
1805 title=plot_title,
1806 method=method,
1807 overlay_cube=overlay_slice,
1808 contour_cube=contour_slice,
1809 point_cube=point_slice,
1810 **kwargs,
1811 )
1812 plot_index.append(plot_filename)
1814 # Add list of plots to plot metadata.
1815 complete_plot_index = _append_to_plot_index(plot_index)
1817 # Make a page to display the plots.
1818 _make_plot_html_page(complete_plot_index)
1821####################
1822# Public functions #
1823####################
1826def spatial_contour_plot(
1827 cube: iris.cube.Cube,
1828 filename: str = None,
1829 sequence_coordinate: str = "time",
1830 stamp_coordinate: str = "realization",
1831 **kwargs,
1832) -> iris.cube.Cube:
1833 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1835 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1836 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1837 is present then postage stamp plots will be produced.
1839 Parameters
1840 ----------
1841 cube: Cube
1842 Iris cube of the data to plot. It should have two spatial dimensions,
1843 such as lat and lon, and may also have a another two dimension to be
1844 plotted sequentially and/or as postage stamp plots.
1845 filename: str, optional
1846 Name of the plot to write, used as a prefix for plot sequences. Defaults
1847 to the recipe name.
1848 sequence_coordinate: str, optional
1849 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1850 This coordinate must exist in the cube.
1851 stamp_coordinate: str, optional
1852 Coordinate about which to plot postage stamp plots. Defaults to
1853 ``"realization"``.
1855 Returns
1856 -------
1857 Cube
1858 The original cube (so further operations can be applied).
1860 Raises
1861 ------
1862 ValueError
1863 If the cube doesn't have the right dimensions.
1864 TypeError
1865 If the cube isn't a single cube.
1866 """
1867 _spatial_plot(
1868 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1869 )
1870 return cube
1873def spatial_pcolormesh_plot(
1874 cube: iris.cube.Cube,
1875 filename: str = None,
1876 sequence_coordinate: str = "time",
1877 stamp_coordinate: str = "realization",
1878 **kwargs,
1879) -> iris.cube.Cube:
1880 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1882 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1883 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1884 is present then postage stamp plots will be produced.
1886 This function is significantly faster than ``spatial_contour_plot``,
1887 especially at high resolutions, and should be preferred unless contiguous
1888 contour areas are important.
1890 Parameters
1891 ----------
1892 cube: Cube
1893 Iris cube of the data to plot. It should have two spatial dimensions,
1894 such as lat and lon, and may also have a another two dimension to be
1895 plotted sequentially and/or as postage stamp plots.
1896 filename: str, optional
1897 Name of the plot to write, used as a prefix for plot sequences. Defaults
1898 to the recipe name.
1899 sequence_coordinate: str, optional
1900 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1901 This coordinate must exist in the cube.
1902 stamp_coordinate: str, optional
1903 Coordinate about which to plot postage stamp plots. Defaults to
1904 ``"realization"``.
1906 Returns
1907 -------
1908 Cube
1909 The original cube (so further operations can be applied).
1911 Raises
1912 ------
1913 ValueError
1914 If the cube doesn't have the right dimensions.
1915 TypeError
1916 If the cube isn't a single cube.
1917 """
1918 _spatial_plot(
1919 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1920 )
1921 return cube
1924def spatial_multi_pcolormesh_plot(
1925 cube: iris.cube.Cube,
1926 overlay_cube: iris.cube.Cube | None = None,
1927 contour_cube: iris.cube.Cube | None = None,
1928 point_cube: iris.cube.Cube | None = None,
1929 filename: str = None,
1930 sequence_coordinate: str = "time",
1931 stamp_coordinate: str = "realization",
1932 **kwargs,
1933) -> iris.cube.Cube:
1934 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
1936 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
1937 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1938 is present then postage stamp plots will be produced.
1940 If specified, a masked overlay_cube can be overplotted on top of the base cube.
1942 If specified, contours of a contour_cube can be overplotted on top of those.
1944 If specified, a spatial scatter map of point_cube can be overplotted.
1946 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
1948 This function is significantly faster than ``spatial_contour_plot``,
1949 especially at high resolutions, and should be preferred unless contiguous
1950 contour areas are important.
1952 Parameters
1953 ----------
1954 cube: Cube
1955 Iris cube of the data to plot. It should have two spatial dimensions,
1956 such as lat and lon, and may also have two additional dimensions to be
1957 plotted sequentially and/or as postage stamp plots.
1958 overlay_cube: Cube, optional
1959 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
1960 such as lat and lon, and may also have two additional dimensions to be
1961 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.
1962 If not provided, output plot generated without overlay cube.
1963 contour_cube: Cube, optional
1964 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,
1965 such as lat and lon, and may also have two additional dimensions to be
1966 plotted sequentially and/or as postage stamp plots. If not provided, output plot generated without contours.
1967 point_cube: Cube, optional
1968 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
1969 spatial dimensions, such as lat and lon, but these can describe a 1-D cube (e.g. list of
1970 observation stations with lat/lon coordinates) and may also have two additional dimensions to be plotted sequentially and/or as
1971 postage stamp plots. If not provided, output plot generated without point-based layer.
1972 filename: str, optional
1973 Name of the plot to write, used as a prefix for plot sequences. Defaults
1974 to the recipe name.
1975 sequence_coordinate: str, optional
1976 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1977 This coordinate must exist in the cube.
1978 stamp_coordinate: str, optional
1979 Coordinate about which to plot postage stamp plots. Defaults to
1980 ``"realization"``.
1982 Returns
1983 -------
1984 Cube
1985 The original cube (so further operations can be applied).
1987 Raises
1988 ------
1989 ValueError
1990 If the cube doesn't have the right dimensions.
1991 TypeError
1992 If the cube isn't a single cube.
1993 """
1994 _spatial_plot(
1995 "pcolormesh",
1996 cube,
1997 filename,
1998 sequence_coordinate,
1999 stamp_coordinate,
2000 overlay_cube=overlay_cube,
2001 contour_cube=contour_cube,
2002 point_cube=point_cube,
2003 )
2004 return cube, overlay_cube, contour_cube, point_cube
2007# TODO: Expand function to handle ensemble data.
2008# line_coordinate: str, optional
2009# Coordinate about which to plot multiple lines. Defaults to
2010# ``"realization"``.
2011def plot_line_series(
2012 cube: iris.cube.Cube | iris.cube.CubeList,
2013 filename: str = None,
2014 series_coordinate: str = "time",
2015 sequence_coordinate: str = "time",
2016 # add the following for ensembles
2017 stamp_coordinate: str = "realization",
2018 single_plot: bool = False,
2019 **kwargs,
2020) -> iris.cube.Cube | iris.cube.CubeList:
2021 """Plot a line plot for the specified coordinate.
2023 The Cube or CubeList must be 1D.
2025 Parameters
2026 ----------
2027 iris.cube | iris.cube.CubeList
2028 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2029 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2030 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2031 filename: str, optional
2032 Name of the plot to write, used as a prefix for plot sequences. Defaults
2033 to the recipe name.
2034 series_coordinate: str, optional
2035 Coordinate about which to make a series. Defaults to ``"time"``. This
2036 coordinate must exist in the cube.
2038 Returns
2039 -------
2040 iris.cube.Cube | iris.cube.CubeList
2041 The original Cube or CubeList (so further operations can be applied).
2043 Raises
2044 ------
2045 ValueError
2046 If the cubes don't have the right dimensions.
2047 TypeError
2048 If the cube isn't a Cube or CubeList.
2049 """
2050 # Ensure we have a name for the plot file.
2051 recipe_title = get_recipe_metadata().get("title", iter_maybe(cube)[0].name())
2053 num_models = get_num_models(cube)
2055 validate_cube_shape(cube, num_models)
2057 # Iterate over all cubes and extract coordinate to plot.
2058 cubes = iter_maybe(cube)
2059 coords = []
2060 for cube in cubes:
2061 try:
2062 coords.append(cube.coord(series_coordinate))
2063 except iris.exceptions.CoordinateNotFoundError as err:
2064 raise ValueError(
2065 f"Cube must have a {series_coordinate} coordinate."
2066 ) from err
2067 if cube.coords("realization"): 2067 ↛ 2071line 2067 didn't jump to line 2071 because the condition on line 2067 was always true
2068 if cube.ndim > 3: 2068 ↛ 2069line 2068 didn't jump to line 2069 because the condition on line 2068 was never true
2069 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2070 else:
2071 raise ValueError("Cube must have a realization coordinate.")
2073 plot_index = []
2075 # Check if this is a spectral plot by looking for spectral coordinates
2076 is_spectral_plot = series_coordinate in [
2077 "frequency",
2078 "physical_wavenumber",
2079 "wavelength",
2080 ]
2082 if is_spectral_plot:
2083 # If series coordinate is frequency, physical_wavenumber or wavelength, for example power spectra with series
2084 # coordinate frequency/wavenumber.
2085 # If several power spectra are plotted with time as sequence_coordinate for the
2086 # time slider option.
2088 # Internal plotting function.
2089 plotting_func = _plot_and_save_line_power_spectrum_series
2091 for cube in cubes:
2092 try:
2093 cube.coord(sequence_coordinate)
2094 except iris.exceptions.CoordinateNotFoundError as err:
2095 raise ValueError(
2096 f"Cube must have a {sequence_coordinate} coordinate."
2097 ) from err
2099 if num_models == 1: 2099 ↛ 2113line 2099 didn't jump to line 2113 because the condition on line 2099 was always true
2100 # check for ensembles
2101 if ( 2101 ↛ 2105line 2101 didn't jump to line 2105 because the condition on line 2101 was never true
2102 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2103 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2104 ):
2105 if single_plot:
2106 # Plot spectra, mean and ensemble spread on 1 plot
2107 plotting_func = _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
2108 else:
2109 # Plot postage stamps
2110 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
2111 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2112 else:
2113 all_points = sorted(
2114 set(
2115 itertools.chain.from_iterable(
2116 cb.coord(sequence_coordinate).points for cb in cubes
2117 )
2118 )
2119 )
2120 all_slices = list(
2121 itertools.chain.from_iterable(
2122 cb.slices_over(sequence_coordinate) for cb in cubes
2123 )
2124 )
2125 # Matched slices (matched by seq coord point; it may happen that
2126 # evaluated models do not cover the same seq coord range, hence matching
2127 # necessary)
2128 cube_iterables = [
2129 iris.cube.CubeList(
2130 s
2131 for s in all_slices
2132 if s.coord(sequence_coordinate).points[0] == point
2133 )
2134 for point in all_points
2135 ]
2137 nplot = np.size(cube.coord(sequence_coordinate).points)
2139 # Create a plot for each value of the sequence coordinate. Allowing for
2140 # multiple cubes in a CubeList to be plotted in the same plot for similar
2141 # sequence values. Passing a CubeList into the internal plotting function
2142 # for similar values of the sequence coordinate. cube_slice can be an
2143 # iris.cube.Cube or an iris.cube.CubeList.
2145 for cube_slice in cube_iterables:
2146 # Normalize cube_slice to a list of cubes
2147 if isinstance(cube_slice, iris.cube.CubeList): 2147 ↛ 2148line 2147 didn't jump to line 2148 because the condition on line 2147 was never true
2148 cubes = list(cube_slice)
2149 elif isinstance(cube_slice, iris.cube.Cube): 2149 ↛ 2152line 2149 didn't jump to line 2152 because the condition on line 2149 was always true
2150 cubes = [cube_slice]
2151 else:
2152 raise TypeError(f"Expected Cube or CubeList, got {type(cube_slice)}")
2154 # Use sequence value so multiple sequences can merge.
2155 seq_coord = cube_slice[0].coord(sequence_coordinate)
2156 plot_title, plot_filename = _set_title_and_filename(
2157 seq_coord, nplot, recipe_title, filename
2158 )
2160 # Format the coordinate value in a unit appropriate way.
2161 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.points[0])}]"
2163 # Use sequence (e.g. time) bounds if plotting single non-sequence outputs
2164 if nplot == 1 and seq_coord.has_bounds: 2164 ↛ 2169line 2164 didn't jump to line 2169 because the condition on line 2164 was always true
2165 if np.size(seq_coord.bounds) > 1: 2165 ↛ 2166line 2165 didn't jump to line 2166 because the condition on line 2165 was never true
2166 title = f"{recipe_title}\n [{seq_coord.units.title(seq_coord.bounds[0][0])} to {seq_coord.units.title(seq_coord.bounds[0][1])}]"
2168 # Do the actual plotting.
2169 plotting_func(
2170 cube_slice,
2171 coords,
2172 stamp_coordinate,
2173 plot_filename,
2174 title,
2175 series_coordinate,
2176 )
2178 plot_index.append(plot_filename)
2179 else:
2180 # Format the title and filename using plotted series coordinate
2181 nplot = 1
2182 seq_coord = coords[0]
2183 plot_title, plot_filename = _set_title_and_filename(
2184 seq_coord, nplot, recipe_title, filename
2185 )
2186 # Do the actual plotting for all other series coordinate options.
2187 _plot_and_save_line_series(
2188 cubes, coords, stamp_coordinate, plot_filename, plot_title
2189 )
2191 plot_index.append(plot_filename)
2193 # append plot to list of plots
2194 complete_plot_index = _append_to_plot_index(plot_index)
2196 # Make a page to display the plots.
2197 _make_plot_html_page(complete_plot_index)
2199 return cube
2202def plot_vertical_line_series(
2203 cubes: iris.cube.Cube | iris.cube.CubeList,
2204 filename: str = None,
2205 series_coordinate: str = "model_level_number",
2206 sequence_coordinate: str = "time",
2207 # line_coordinate: str = "realization",
2208 **kwargs,
2209) -> iris.cube.Cube | iris.cube.CubeList:
2210 """Plot a line plot against a type of vertical coordinate.
2212 The Cube or CubeList must be 1D.
2214 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2215 then a sequence of plots will be produced.
2217 Parameters
2218 ----------
2219 iris.cube | iris.cube.CubeList
2220 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2221 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2222 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2223 filename: str, optional
2224 Name of the plot to write, used as a prefix for plot sequences. Defaults
2225 to the recipe name.
2226 series_coordinate: str, optional
2227 Coordinate to plot on the y-axis. Can be ``pressure`` or
2228 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2229 for LFRic. Defaults to ``model_level_number``.
2230 This coordinate must exist in the cube.
2231 sequence_coordinate: str, optional
2232 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2233 This coordinate must exist in the cube.
2235 Returns
2236 -------
2237 iris.cube.Cube | iris.cube.CubeList
2238 The original Cube or CubeList (so further operations can be applied).
2239 Plotted data.
2241 Raises
2242 ------
2243 ValueError
2244 If the cubes doesn't have the right dimensions.
2245 TypeError
2246 If the cube isn't a Cube or CubeList.
2247 """
2248 # Ensure we have a name for the plot file.
2249 recipe_title = get_recipe_metadata().get("title", iter_maybe(cubes)[0].name())
2251 cubes = iter_maybe(cubes)
2252 # Initialise empty list to hold all data from all cubes in a CubeList
2253 all_data = []
2255 # Store min/max ranges for x range.
2256 x_levels = []
2258 num_models = get_num_models(cubes)
2260 validate_cube_shape(cubes, num_models)
2262 # Iterate over all cubes in cube or CubeList and plot.
2263 coords = []
2264 for cube in cubes:
2265 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2266 try:
2267 coords.append(cube.coord(series_coordinate))
2268 except iris.exceptions.CoordinateNotFoundError as err:
2269 raise ValueError(
2270 f"Cube must have a {series_coordinate} coordinate."
2271 ) from err
2273 try:
2274 if cube.ndim > 1 or not cube.coords("realization"): 2274 ↛ 2282line 2274 didn't jump to line 2282 because the condition on line 2274 was always true
2275 cube.coord(sequence_coordinate)
2276 except iris.exceptions.CoordinateNotFoundError as err:
2277 raise ValueError(
2278 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2279 ) from err
2281 # Get minimum and maximum from levels information.
2282 _, levels, _ = colorbar_map_levels(cube, axis="x")
2283 if levels is not None: 2283 ↛ 2287line 2283 didn't jump to line 2287 because the condition on line 2283 was always true
2284 x_levels.append(min(levels))
2285 x_levels.append(max(levels))
2286 else:
2287 all_data.append(cube.data)
2289 if len(x_levels) == 0: 2289 ↛ 2291line 2289 didn't jump to line 2291 because the condition on line 2289 was never true
2290 # Combine all data into a single NumPy array
2291 combined_data = np.concatenate(all_data)
2293 # Set the lower and upper limit for the x-axis to ensure all plots have
2294 # same range. This needs to read the whole cube over the range of the
2295 # sequence and if applicable postage stamp coordinate.
2296 vmin = np.floor(combined_data.min())
2297 vmax = np.ceil(combined_data.max())
2298 else:
2299 vmin = min(x_levels)
2300 vmax = max(x_levels)
2302 # Matching the slices (matching by seq coord point; it may happen that
2303 # evaluated models do not cover the same seq coord range, hence matching
2304 # necessary)
2305 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2307 # Create a plot for each value of the sequence coordinate.
2308 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2309 # similar sequence values. Passing a CubeList into the internal plotting function
2310 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2311 # or an iris.cube.CubeList.
2312 plot_index = []
2313 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2314 for cubes_slice in cube_iterables:
2315 # Format the coordinate value in a unit appropriate way.
2316 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2317 plot_title, plot_filename = _set_title_and_filename(
2318 seq_coord, nplot, recipe_title, filename
2319 )
2321 # Do the actual plotting.
2322 _plot_and_save_vertical_line_series(
2323 cubes_slice,
2324 coords,
2325 "realization",
2326 plot_filename,
2327 series_coordinate,
2328 title=plot_title,
2329 vmin=vmin,
2330 vmax=vmax,
2331 )
2332 plot_index.append(plot_filename)
2334 # Add list of plots to plot metadata.
2335 complete_plot_index = _append_to_plot_index(plot_index)
2337 # Make a page to display the plots.
2338 _make_plot_html_page(complete_plot_index)
2340 return cubes
2343def qq_plot(
2344 cubes: iris.cube.CubeList,
2345 coordinates: list[str],
2346 percentiles: list[float],
2347 model_names: list[str],
2348 filename: str = None,
2349 one_to_one: bool = True,
2350 **kwargs,
2351) -> iris.cube.CubeList:
2352 """Plot a Quantile-Quantile plot between two models for common time points.
2354 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2355 collapsed within the operator over all specified coordinates such as
2356 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2357 ensemble members to ensure a 1D cube (array).
2359 Parameters
2360 ----------
2361 cubes: iris.cube.CubeList
2362 Two cubes of the same variable with different models.
2363 coordinate: list[str]
2364 The list of coordinates to collapse over. This list should be
2365 every coordinate within the cube to result in a 1D cube around
2366 the percentile coordinate.
2367 percent: list[float]
2368 A list of percentiles to appear in the plot.
2369 model_names: list[str]
2370 A list of model names to appear on the axis of the plot.
2371 filename: str, optional
2372 Filename of the plot to write.
2373 one_to_one: bool, optional
2374 If True a 1:1 line is plotted; if False it is not. Default is True.
2376 Raises
2377 ------
2378 ValueError
2379 When the cubes are not compatible.
2381 Notes
2382 -----
2383 The quantile-quantile plot is a variant on the scatter plot representing
2384 two datasets by their quantiles (percentiles) for common time points.
2385 This plot does not use a theoretical distribution to compare against, but
2386 compares percentiles of two datasets. This plot does
2387 not use all raw data points, but plots the selected percentiles (quantiles) of
2388 each variable instead for the two datasets, thereby normalising the data for a
2389 direct comparison between the selected percentiles of the two dataset distributions.
2391 Quantile-quantile plots are valuable for comparing against
2392 observations and other models. Identical percentiles between the variables
2393 will lie on the one-to-one line implying the values correspond well to each
2394 other. Where there is a deviation from the one-to-one line a range of
2395 possibilities exist depending on how and where the data is shifted (e.g.,
2396 Wilks 2011 [Wilks2011]_).
2398 For distributions above the one-to-one line the distribution is left-skewed;
2399 below is right-skewed. A distinct break implies a bimodal distribution, and
2400 closer values/values further apart at the tails imply poor representation of
2401 the extremes.
2403 References
2404 ----------
2405 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2406 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2407 """
2408 # Check cubes using same functionality as the difference operator.
2409 if len(cubes) != 2:
2410 raise ValueError("cubes should contain exactly 2 cubes.")
2411 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2412 other: Cube = cubes.extract_cube(
2413 iris.Constraint(
2414 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2415 )
2416 )
2418 # Get spatial coord names.
2419 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2420 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2422 # Ensure cubes to compare are on common differencing grid.
2423 # This is triggered if either
2424 # i) latitude and longitude shapes are not the same. Note grid points
2425 # are not compared directly as these can differ through rounding
2426 # errors.
2427 # ii) or variables are known to often sit on different grid staggering
2428 # in different models (e.g. cell center vs cell edge), as is the case
2429 # for UM and LFRic comparisons.
2430 # In future greater choice of regridding method might be applied depending
2431 # on variable type. Linear regridding can in general be appropriate for smooth
2432 # variables. Care should be taken with interpretation of differences
2433 # given this dependency on regridding.
2434 if (
2435 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2436 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2437 ) or (
2438 base.long_name
2439 in [
2440 "eastward_wind_at_10m",
2441 "northward_wind_at_10m",
2442 "northward_wind_at_cell_centres",
2443 "eastward_wind_at_cell_centres",
2444 "zonal_wind_at_pressure_levels",
2445 "meridional_wind_at_pressure_levels",
2446 "potential_vorticity_at_pressure_levels",
2447 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2448 ]
2449 ):
2450 logging.debug(
2451 "Linear regridding base cube to other grid to compute differences"
2452 )
2453 base = regrid_onto_cube(base, other, method="Linear")
2455 # Extract just common time points.
2456 base, other = _extract_common_time_points(base, other)
2458 # Equalise attributes so we can merge.
2459 fully_equalise_attributes([base, other])
2460 logging.debug("Base: %s\nOther: %s", base, other)
2462 # Collapse cubes.
2463 base = collapse(
2464 base,
2465 coordinate=coordinates,
2466 method="PERCENTILE",
2467 additional_percent=percentiles,
2468 )
2469 other = collapse(
2470 other,
2471 coordinate=coordinates,
2472 method="PERCENTILE",
2473 additional_percent=percentiles,
2474 )
2476 # Ensure we have a name for the plot file.
2477 recipe_title = get_recipe_metadata().get("title", "QQ_plot")
2478 title = f"{recipe_title}"
2480 if filename is None:
2481 filename = slugify(recipe_title)
2483 # Add file extension.
2484 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2486 # Do the actual plotting on a scatter plot
2487 _plot_and_save_scatter_plot(
2488 base, other, plot_filename, title, one_to_one, model_names
2489 )
2491 # Add list of plots to plot metadata.
2492 plot_index = _append_to_plot_index([plot_filename])
2494 # Make a page to display the plots.
2495 _make_plot_html_page(plot_index)
2497 return iris.cube.CubeList([base, other])
2500def hinton_plot(change, signif, xaxis_labels, yaxis_labels, magnitude=None):
2501 """
2502 Plot a Hinton style triangle/scorecard plot.
2504 This plot type can be useful for summarising high level information, such as comparing
2505 how 'skillful' two models are when verified against observations for a variety of metrics,
2506 as a function of lead-time. A few parameters of the plot style are fixed in function rather
2507 than customisable by the user as input arguments; many have been designed to automatically
2508 scale the plot depending on the number of x and y components.
2510 Parameters
2511 ----------
2512 change: np.ndarray
2513 A 2d numpy array containing the values (scaled to 1 to -1) that determine the triangle
2514 size/direction.
2515 signif: np.ndarray
2516 A 2d numpy array containing 0s and 1s to determine if triangle is significant or not.
2517 xaxis_labels: list
2518 List of labels for the xaxis (must match the second dimension length of signif and change,
2519 along with magnitude if not None).
2520 yaxis_labels: list
2521 List of labels for the yaxis (must match the first dimension length of signif and change,
2522 along with magnitude if not None).
2523 magnitude: np.ndarray | None
2524 Optional 2D array, matching the shape of change, signif, which contains numerical values
2525 the user wishes to display under each respective triangle.
2527 Returns
2528 -------
2529 matplotlib axes object to either display or do further modifications to.
2530 """
2531 # Setup colors of triangles
2532 color_pos = "#7CAE00"
2533 color_neg = "#7B68EE"
2535 # Setup cell/text size ratios
2536 figsize = None
2537 cell_size_in = 0.35
2538 text_row_ratio = 0.25
2540 # Ensure arrays, and change to bool for sig.
2541 change = np.asarray(change)
2542 signif = np.asarray(signif).astype(bool)
2543 if magnitude is not None: 2543 ↛ 2544line 2543 didn't jump to line 2544 because the condition on line 2543 was never true
2544 magnitude = np.asarray(magnitude)
2546 # Get the number of x and y elements
2547 ny, nx = change.shape
2549 # Build non-uniform y coordinates
2550 tri_height = 1.0
2551 txt_height = text_row_ratio
2553 tri_y = []
2554 txt_y = []
2555 y_edges = [0.0]
2557 y = 0.0
2558 for _j in range(ny):
2559 tri_y.append(y + tri_height / 2)
2560 y += tri_height
2561 y_edges.append(y)
2563 if magnitude is not None: 2563 ↛ 2564line 2563 didn't jump to line 2564 because the condition on line 2563 was never true
2564 txt_y.append(y + txt_height / 2)
2565 y += txt_height
2566 y_edges.append(y)
2568 total_height = y
2570 # Dynamic figure size
2571 if figsize is None: 2571 ↛ 2576line 2571 didn't jump to line 2576 because the condition on line 2571 was always true
2572 width = nx * cell_size_in
2573 height = total_height * cell_size_in + 2
2574 figsize = (width, height)
2576 fig, ax = plt.subplots(figsize=figsize)
2578 # Setup axes and grid.
2579 ax.set_aspect("equal", adjustable="box")
2580 ax.set_xlim(-0.5, nx - 0.5)
2581 ax.set_ylim(0, total_height)
2583 ax.set_xticks(np.arange(nx))
2584 ax.set_xticklabels(xaxis_labels, rotation=90)
2586 ax.set_yticks(tri_y)
2587 ax.set_yticklabels(yaxis_labels)
2589 ax.set_xticks(np.arange(-0.5, nx, 1), minor=True)
2590 ax.set_yticks(y_edges, minor=True)
2592 ax.set_axisbelow(True)
2593 ax.grid(which="minor", linestyle=":", linewidth=0.3, color="0.7")
2594 ax.grid(False, which="major")
2595 ax.tick_params(which="minor", length=0)
2597 ax.invert_yaxis()
2599 # Compute marker scaling (fixed overlap)
2600 fig.canvas.draw()
2602 bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
2603 width_in, height_in = bbox.width, bbox.height
2605 cell_w = (width_in * fig.dpi) / nx
2606 cell_h = (height_in * fig.dpi) / total_height
2607 cell_pixels = min(cell_w, cell_h)
2609 max_marker_size = (0.6 * cell_pixels) ** 2
2611 text_fontsize = cell_pixels * 0.15
2613 # Plot triangles + text
2614 for j in range(ny):
2615 for i in range(nx):
2616 val = change[j, i]
2617 if np.isnan(val): 2617 ↛ 2618line 2617 didn't jump to line 2618 because the condition on line 2617 was never true
2618 continue
2620 if abs(val) < 0.01: 2620 ↛ 2621line 2620 didn't jump to line 2621 because the condition on line 2620 was never true
2621 continue
2623 sig = signif[j, i]
2624 size = max_marker_size * abs(val)
2626 # Triangle style
2627 if val >= 0:
2628 marker = "^"
2629 color = color_pos
2630 else:
2631 marker = "v"
2632 color = color_neg
2634 if sig:
2635 edgecolor = "black"
2636 linewidth = 0.6
2637 else:
2638 edgecolor = "none"
2639 linewidth = 0.0
2641 # Triangle
2642 ax.scatter(
2643 i,
2644 tri_y[j],
2645 s=size,
2646 marker=marker,
2647 c=color,
2648 edgecolors=edgecolor,
2649 linewidths=linewidth,
2650 zorder=3,
2651 clip_on=True, # ensures no rendering bleed
2652 )
2654 # Text row
2655 if magnitude is not None: 2655 ↛ 2656line 2655 didn't jump to line 2656 because the condition on line 2655 was never true
2656 mag_val = magnitude[j, i]
2658 if not np.isnan(mag_val):
2659 ax.text(
2660 i,
2661 txt_y[j],
2662 f"{mag_val:.1f}",
2663 ha="center",
2664 va="center",
2665 fontsize=text_fontsize,
2666 color="black",
2667 zorder=4,
2668 )
2670 plt.tight_layout()
2671 return fig, ax
2674def scatter_plot(
2675 cube_x: iris.cube.Cube | iris.cube.CubeList,
2676 cube_y: iris.cube.Cube | iris.cube.CubeList,
2677 filename: str = None,
2678 one_to_one: bool = True,
2679 **kwargs,
2680) -> iris.cube.CubeList:
2681 """Plot a scatter plot between two variables.
2683 Both cubes must be 1D.
2685 Parameters
2686 ----------
2687 cube_x: Cube | CubeList
2688 1 dimensional Cube of the data to plot on y-axis.
2689 cube_y: Cube | CubeList
2690 1 dimensional Cube of the data to plot on x-axis.
2691 filename: str, optional
2692 Filename of the plot to write.
2693 one_to_one: bool, optional
2694 If True a 1:1 line is plotted; if False it is not. Default is True.
2696 Returns
2697 -------
2698 cubes: CubeList
2699 CubeList of the original x and y cubes for further processing.
2701 Raises
2702 ------
2703 ValueError
2704 If the cube doesn't have the right dimensions and cubes not the same
2705 size.
2706 TypeError
2707 If the cube isn't a single cube.
2709 Notes
2710 -----
2711 Scatter plots are used for determining if there is a relationship between
2712 two variables. Positive relations have a slope going from bottom left to top
2713 right; Negative relations have a slope going from top left to bottom right.
2714 """
2715 # Iterate over all cubes in cube or CubeList and plot.
2716 for cube_iter in iter_maybe(cube_x):
2717 # Check cubes are correct shape.
2718 cube_iter = check_single_cube(cube_iter)
2719 if cube_iter.ndim > 1:
2720 raise ValueError("cube_x must be 1D.")
2722 # Iterate over all cubes in cube or CubeList and plot.
2723 for cube_iter in iter_maybe(cube_y):
2724 # Check cubes are correct shape.
2725 cube_iter = check_single_cube(cube_iter)
2726 if cube_iter.ndim > 1:
2727 raise ValueError("cube_y must be 1D.")
2729 # Ensure we have a name for the plot file.
2730 recipe_title = get_recipe_metadata().get("title", "Scatter_plot")
2731 title = f"{recipe_title}"
2733 if filename is None:
2734 filename = slugify(recipe_title)
2736 # Add file extension.
2737 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2739 # Do the actual plotting.
2740 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2742 # Add list of plots to plot metadata.
2743 plot_index = _append_to_plot_index([plot_filename])
2745 # Make a page to display the plots.
2746 _make_plot_html_page(plot_index)
2748 return iris.cube.CubeList([cube_x, cube_y])
2751def vector_plot(
2752 cube_u: iris.cube.Cube,
2753 cube_v: iris.cube.Cube,
2754 filename: str = None,
2755 sequence_coordinate: str = "time",
2756 **kwargs,
2757) -> iris.cube.CubeList:
2758 """Plot a vector plot based on the input u and v components."""
2759 recipe_title = get_recipe_metadata().get("title", "Vector_plot")
2761 # Cubes must have a matching sequence coordinate.
2762 try:
2763 # Check that the u and v cubes have the same sequence coordinate.
2764 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2764 ↛ anywhereline 2764 didn't jump anywhere: it always raised an exception.
2765 raise ValueError("Coordinates do not match.")
2766 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2767 raise ValueError(
2768 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2769 ) from err
2771 # Create a plot for each value of the sequence coordinate.
2772 plot_index = []
2773 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2774 for cube_u_slice, cube_v_slice in zip(
2775 cube_u.slices_over(sequence_coordinate),
2776 cube_v.slices_over(sequence_coordinate),
2777 strict=True,
2778 ):
2779 # Format the coordinate value in a unit appropriate way.
2780 seq_coord = cube_u_slice.coord(sequence_coordinate)
2781 plot_title, plot_filename = _set_title_and_filename(
2782 seq_coord, nplot, recipe_title, filename
2783 )
2785 # Do the actual plotting.
2786 _plot_and_save_vector_plot(
2787 cube_u_slice,
2788 cube_v_slice,
2789 filename=plot_filename,
2790 title=plot_title,
2791 method="pcolormesh",
2792 )
2793 plot_index.append(plot_filename)
2795 # Add list of plots to plot metadata.
2796 complete_plot_index = _append_to_plot_index(plot_index)
2798 # Make a page to display the plots.
2799 _make_plot_html_page(complete_plot_index)
2801 return iris.cube.CubeList([cube_u, cube_v])
2804def plot_histogram_series(
2805 cubes: iris.cube.Cube | iris.cube.CubeList,
2806 filename: str = None,
2807 sequence_coordinate: str = "time",
2808 stamp_coordinate: str = "realization",
2809 single_plot: bool = False,
2810 **kwargs,
2811) -> iris.cube.Cube | iris.cube.CubeList:
2812 """Plot a histogram plot for each vertical level provided.
2814 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2815 is present then a sequence of plots will be produced using the time slider
2816 functionality to scroll through histograms against time. If a
2817 stamp_coordinate is present then postage stamp plots will be produced. If
2818 stamp_coordinate and single_plot is True, all postage stamp plots will be
2819 plotted in a single plot instead of separate postage stamp plots.
2821 Parameters
2822 ----------
2823 cubes: Cube | iris.cube.CubeList
2824 Iris cube or CubeList of the data to plot. It should have a single dimension other
2825 than the stamp coordinate.
2826 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2827 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2828 filename: str, optional
2829 Name of the plot to write, used as a prefix for plot sequences. Defaults
2830 to the recipe name.
2831 sequence_coordinate: str, optional
2832 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2833 This coordinate must exist in the cube and will be used for the time
2834 slider.
2835 stamp_coordinate: str, optional
2836 Coordinate about which to plot postage stamp plots. Defaults to
2837 ``"realization"``.
2838 single_plot: bool, optional
2839 If True, all postage stamp plots will be plotted in a single plot. If
2840 False, each postage stamp plot will be plotted separately. Is only valid
2841 if stamp_coordinate exists and has more than a single point.
2843 Returns
2844 -------
2845 iris.cube.Cube | iris.cube.CubeList
2846 The original Cube or CubeList (so further operations can be applied).
2847 Plotted data.
2849 Raises
2850 ------
2851 ValueError
2852 If the cube doesn't have the right dimensions.
2853 TypeError
2854 If the cube isn't a Cube or CubeList.
2855 """
2856 recipe_title = get_recipe_metadata().get("title", "Histogram")
2858 cubes = iter_maybe(cubes)
2859 # Ensure we have a name for the plot file.
2860 if filename is None:
2861 filename = slugify(recipe_title)
2863 # Internal plotting function.
2864 plotting_func = _plot_and_save_histogram_series
2866 num_models = get_num_models(cubes)
2868 validate_cube_shape(cubes, num_models)
2870 # If several histograms are plotted, check sequence_coordinate
2871 check_sequence_coordinate(cubes, sequence_coordinate)
2873 # Get axis minimum and maximum from levels information.
2874 # If no levels set, derive minima and maxima from data in CubeList.
2875 vmin, vmax = _set_axis_range(cubes)
2877 # Make postage stamp plots if stamp_coordinate exists and has more than a
2878 # single point. If single_plot is True:
2879 # -- all postage stamp plots will be plotted in a single plot instead of
2880 # separate postage stamp plots.
2881 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2882 # produced per single model only
2883 if num_models == 1:
2884 if ( 2884 ↛ 2888line 2884 didn't jump to line 2888 because the condition on line 2884 was never true
2885 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2886 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2887 ):
2888 if single_plot:
2889 plotting_func = (
2890 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2891 )
2892 else:
2893 plotting_func = _plot_and_save_postage_stamp_histogram_series
2894 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2895 else:
2896 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2898 plot_index = []
2899 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2900 # Create a plot for each value of the sequence coordinate. Allowing for
2901 # multiple cubes in a CubeList to be plotted in the same plot for similar
2902 # sequence values. Passing a CubeList into the internal plotting function
2903 # for similar values of the sequence coordinate. cube_slice can be an
2904 # iris.cube.Cube or an iris.cube.CubeList.
2905 for cube_slice in cube_iterables:
2906 single_cube = cube_slice
2907 if isinstance(cube_slice, iris.cube.CubeList):
2908 single_cube = cube_slice[0]
2910 # Ensure valid stamp coordinate in cube dimensions
2911 if stamp_coordinate == "realization": 2911 ↛ 2914line 2911 didn't jump to line 2914 because the condition on line 2911 was always true
2912 stamp_coordinate = check_stamp_coordinate(single_cube)
2913 # Set plot titles and filename, based on sequence coordinate
2914 seq_coord = single_cube.coord(sequence_coordinate)
2915 # Use time coordinate in title and filename if single histogram output.
2916 if sequence_coordinate == "realization" and nplot == 1: 2916 ↛ 2917line 2916 didn't jump to line 2917 because the condition on line 2916 was never true
2917 seq_coord = single_cube.coord("time")
2918 plot_title, plot_filename = _set_title_and_filename(
2919 seq_coord, nplot, recipe_title, filename
2920 )
2922 # Do the actual plotting.
2923 plotting_func(
2924 cube_slice,
2925 filename=plot_filename,
2926 stamp_coordinate=stamp_coordinate,
2927 title=plot_title,
2928 vmin=vmin,
2929 vmax=vmax,
2930 )
2931 plot_index.append(plot_filename)
2933 # Add list of plots to plot metadata.
2934 complete_plot_index = _append_to_plot_index(plot_index)
2936 # Make a page to display the plots.
2937 _make_plot_html_page(complete_plot_index)
2939 return cubes
2942def _plot_and_save_postage_stamp_power_spectrum_series(
2943 cubes: iris.cube.Cube,
2944 coords: list[iris.coords.Coord],
2945 stamp_coordinate: str,
2946 filename: str,
2947 title: str,
2948 series_coordinate: str = None,
2949 **kwargs,
2950):
2951 """Plot and save postage (ensemble members) stamps for a power spectrum series.
2953 Parameters
2954 ----------
2955 cubes: Cube or CubeList
2956 Cube or Cubelist of the power spectrum data.
2957 coords: list[Coord]
2958 Coordinates to plot on the x-axis, one per cube.
2959 stamp_coordinate: str
2960 Coordinate that becomes different plots.
2961 filename: str
2962 Filename of the plot to write.
2963 title: str
2964 Plot title.
2965 series_coordinate: str, optional
2966 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
2968 """
2969 # Use the smallest square grid that will fit the members.
2970 grid_size = int(math.ceil(math.sqrt(len(cubes.coord(stamp_coordinate).points))))
2972 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
2973 model_colors_map = get_model_colors_map(cubes)
2974 # ax = plt.gca()
2975 # Make a subplot for each member.
2976 for member, subplot in zip(
2977 cubes.slices_over(stamp_coordinate), range(1, grid_size**2 + 1), strict=False
2978 ):
2979 ax = plt.subplot(grid_size, grid_size, subplot)
2981 # Store min/max ranges.
2982 y_levels = []
2984 line_marker = None
2985 line_width = 1
2987 for cube in iter_maybe(member):
2988 xcoord = _select_series_coord(cube, series_coordinate)
2989 xname = xcoord.points
2991 yfield = cube.data # power spectrum
2992 label = None
2993 color = "black"
2994 if model_colors_map: 2994 ↛ 2995line 2994 didn't jump to line 2995 because the condition on line 2994 was never true
2995 label = cube.attributes.get("model_name")
2996 color = model_colors_map.get(label)
2998 if member.coord(stamp_coordinate).points == [0]:
2999 ax.plot(
3000 xname,
3001 yfield,
3002 color=color,
3003 marker=line_marker,
3004 ls="-",
3005 lw=line_width,
3006 label=f"{label} (control)"
3007 if len(cube.coord(stamp_coordinate).points) > 1
3008 else label,
3009 )
3010 # Label with member if part of an ensemble and not the control.
3011 else:
3012 ax.plot(
3013 xname,
3014 yfield,
3015 color=color,
3016 ls="-",
3017 lw=1.5,
3018 alpha=0.75,
3019 label=f"{label} (member)",
3020 )
3022 # Calculate the global min/max if multiple cubes are given.
3023 _, levels, _ = colorbar_map_levels(cube, axis="y")
3024 if levels is not None: 3024 ↛ 3025line 3024 didn't jump to line 3025 because the condition on line 3024 was never true
3025 y_levels.append(min(levels))
3026 y_levels.append(max(levels))
3028 # Add some labels and tweak the style.
3029 title = f"{title}"
3030 ax.set_title(title, fontsize=16)
3032 # Set appropriate x-axis label based on coordinate
3033 if series_coordinate == "wavelength" or ( 3033 ↛ 3036line 3033 didn't jump to line 3036 because the condition on line 3033 was never true
3034 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3035 ):
3036 ax.set_xlabel("Wavelength (km)", fontsize=14)
3037 elif series_coordinate == "physical_wavenumber" or ( 3037 ↛ 3042line 3037 didn't jump to line 3042 because the condition on line 3037 was always true
3038 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3039 ):
3040 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3041 else: # frequency or check units
3042 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3043 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3044 else:
3045 ax.set_xlabel("Wavenumber", fontsize=14)
3047 ax.set_ylabel("Power Spectral Density", fontsize=14)
3048 ax.tick_params(axis="both", labelsize=12)
3050 # Set log-log scale
3051 ax.set_xscale("log")
3052 ax.set_yscale("log")
3054 # Add gridlines
3055 ax.grid(linestyle="--", color="grey", linewidth=1)
3056 # Ientify unique labels for legend
3057 handles = list(
3058 {
3059 label: handle
3060 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3061 }.values()
3062 )
3063 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3065 ax = plt.gca()
3066 ax.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
3068 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3069 logging.info("Saved histogram postage stamp plot to %s", filename)
3070 plt.close(fig)
3073def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
3074 cubes: iris.cube.Cube,
3075 coords: list[iris.coords.Coord],
3076 stamp_coordinate: str,
3077 filename: str,
3078 title: str,
3079 series_coordinate: str = None,
3080 **kwargs,
3081):
3082 """Plot and save power spectra for ensemble members in single plot.
3084 Parameters
3085 ----------
3086 cubes: Cube or CubeList
3087 Cube or Cubelist of the power spectrum data.
3088 coords: list[Coord]
3089 Coordinates to plot on the x-axis, one per cube.
3090 stamp_coordinate: str
3091 Coordinate that becomes different plots.
3092 filename: str
3093 Filename of the plot to write.
3094 title: str
3095 Plot title.
3096 series_coordinate: str, optional
3097 Coordinate being plotted on x-axis. In case of spectra frequency, physical_wavenumber, or wavelength.
3099 """
3100 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
3101 model_colors_map = get_model_colors_map(cubes)
3103 line_marker = None
3104 line_width = 1
3106 # Compute ensemble statistics to show spread
3107 mean_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MEAN)
3108 min_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MIN)
3109 max_cube = cubes.collapsed(stamp_coordinate, iris.analysis.MAX)
3111 xcoord_global = mean_cube.coord(series_coordinate)
3112 x_global = xcoord_global.points
3114 for i, member in enumerate(cubes.slices_over(stamp_coordinate)):
3115 xcoord = _select_series_coord(member, series_coordinate)
3116 xname = xcoord.points
3118 yfield = member.data # power spectrum
3119 color = "black"
3120 if model_colors_map: 3120 ↛ 3124line 3120 didn't jump to line 3124 because the condition on line 3120 was always true
3121 label = member.attributes.get("model_name") if i == 0 else None
3122 color = model_colors_map.get(label)
3124 if member.coord(stamp_coordinate).points == [0]:
3125 ax.plot(
3126 xname,
3127 yfield,
3128 color=color,
3129 marker=line_marker,
3130 ls="-",
3131 lw=line_width,
3132 label=f"{label} (control)"
3133 if len(member.coord(stamp_coordinate).points) > 1
3134 else label,
3135 )
3136 # Label with member number if part of an ensemble and not the control.
3137 else:
3138 ax.plot(
3139 xname,
3140 yfield,
3141 color=color,
3142 ls="-",
3143 lw=1.5,
3144 alpha=0.75,
3145 label=label,
3146 )
3148 # Set appropriate x-axis label based on coordinate
3149 if series_coordinate == "wavelength" or ( 3149 ↛ 3152line 3149 didn't jump to line 3152 because the condition on line 3149 was never true
3150 hasattr(xcoord, "long_name") and xcoord.long_name == "wavelength"
3151 ):
3152 ax.set_xlabel("Wavelength (km)", fontsize=14)
3153 elif series_coordinate == "physical_wavenumber" or ( 3153 ↛ 3158line 3153 didn't jump to line 3158 because the condition on line 3153 was always true
3154 hasattr(xcoord, "long_name") and xcoord.long_name == "physical_wavenumber"
3155 ):
3156 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3157 else: # frequency or check units
3158 if hasattr(xcoord, "units") and str(xcoord.units) == "km-1":
3159 ax.set_xlabel("Wavenumber (km⁻¹)", fontsize=14)
3160 else:
3161 ax.set_xlabel("Wavenumber", fontsize=14)
3163 # Add ensemble spread shading
3164 ax.fill_between(
3165 x_global,
3166 min_cube.data,
3167 max_cube.data,
3168 color="grey",
3169 alpha=0.3,
3170 label="Ensemble spread",
3171 )
3173 # Add ensemble mean line
3174 ax.plot(x_global, mean_cube.data, color="black", lw=1, label="Ensemble mean")
3176 ax.set_ylabel("Power Spectral Density", fontsize=14)
3177 ax.tick_params(axis="both", labelsize=12)
3179 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
3180 # Set log-log scale
3181 ax.set_xscale("log")
3182 ax.set_yscale("log")
3184 # Add gridlines
3185 ax.grid(linestyle="--", color="grey", linewidth=1)
3186 # Identify unique labels for legend
3187 handles = list(
3188 {
3189 label: handle
3190 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
3191 }.values()
3192 )
3193 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
3195 # Figure title.
3196 ax.set_title(title, fontsize=16)
3198 # Save the figure to a file
3199 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
3201 # Close the figure
3202 plt.close(fig)