Coverage for src / CSET / operators / plot.py: 76%
962 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-08 17:20 +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
36import scipy.fft as fft
37from cartopy.mpl.geoaxes import GeoAxes
38from iris.cube import Cube
39from markdown_it import MarkdownIt
40from mpl_toolkits.axes_grid1.inset_locator import inset_axes
42from CSET._common import (
43 filename_slugify,
44 get_recipe_metadata,
45 iter_maybe,
46 render_file,
47 slugify,
48)
49from CSET.operators._plot_colormaps import (
50 _colorbar_map_levels,
51 _get_model_colors_map,
52)
53from CSET.operators._utils import (
54 check_stamp_coordinate,
55 fully_equalise_attributes,
56 get_cube_yxcoordname,
57 is_transect,
58 slice_over_maybe,
59)
60from CSET.operators.collapse import collapse
61from CSET.operators.misc import _extract_common_time_points
62from CSET.operators.regrid import regrid_onto_cube
64# Use a non-interactive plotting backend.
65mpl.use("agg")
67############################
68# Private helper functions #
69############################
72def _append_to_plot_index(plot_index: list) -> list:
73 """Add plots into the plot index, returning the complete plot index."""
74 with open("meta.json", "r+t", encoding="UTF-8") as fp:
75 fcntl.flock(fp, fcntl.LOCK_EX)
76 fp.seek(0)
77 meta = json.load(fp)
78 complete_plot_index = meta.get("plots", [])
79 complete_plot_index = complete_plot_index + plot_index
80 meta["plots"] = complete_plot_index
81 if os.getenv("CYLC_TASK_CYCLE_POINT") and not bool(
82 os.getenv("DO_CASE_AGGREGATION")
83 ):
84 meta["case_date"] = os.getenv("CYLC_TASK_CYCLE_POINT", "")
85 fp.seek(0)
86 fp.truncate()
87 json.dump(meta, fp, indent=2)
88 return complete_plot_index
91def _check_single_cube(cube: iris.cube.Cube | iris.cube.CubeList) -> iris.cube.Cube:
92 """Ensure a single cube is given.
94 If a CubeList of length one is given that the contained cube is returned,
95 otherwise an error is raised.
97 Parameters
98 ----------
99 cube: Cube | CubeList
100 The cube to check.
102 Returns
103 -------
104 cube: Cube
105 The checked cube.
107 Raises
108 ------
109 TypeError
110 If the input cube is not a Cube or CubeList of a single Cube.
111 """
112 if isinstance(cube, iris.cube.Cube):
113 return cube
114 if isinstance(cube, iris.cube.CubeList):
115 if len(cube) == 1:
116 return cube[0]
117 raise TypeError("Must have a single cube", cube)
120def _make_plot_html_page(plots: list):
121 """Create a HTML page to display a plot image."""
122 # Debug check that plots actually contains some strings.
123 assert isinstance(plots[0], str)
125 # Load HTML template file.
126 operator_files = importlib.resources.files()
127 template_file = operator_files.joinpath("_plot_page_template.html")
129 # Get some metadata.
130 meta = get_recipe_metadata()
131 title = meta.get("title", "Untitled")
132 description = MarkdownIt().render(meta.get("description", "*No description.*"))
134 # Prepare template variables.
135 variables = {
136 "title": title,
137 "description": description,
138 "initial_plot": plots[0],
139 "plots": plots,
140 "title_slug": slugify(title),
141 }
143 # Render template.
144 html = render_file(template_file, **variables)
146 # Save completed HTML.
147 with open("index.html", "wt", encoding="UTF-8") as fp:
148 fp.write(html)
151def _setup_spatial_map(
152 cube: iris.cube.Cube,
153 figure,
154 cmap,
155 grid_size: tuple[int, int] | None = None,
156 subplot: int | None = None,
157):
158 """Define map projections, extent and add coastlines and borderlines for spatial plots.
160 For spatial map plots, a relevant map projection for rotated or non-rotated inputs
161 is specified, and map extent defined based on the input data.
163 Parameters
164 ----------
165 cube: Cube
166 2 dimensional (lat and lon) Cube of the data to plot.
167 figure:
168 Matplotlib Figure object holding all plot elements.
169 cmap:
170 Matplotlib colormap.
171 grid_size: (int, int), optional
172 Size of grid (rows, cols) for subplots if multiple spatial subplots in figure.
173 subplot: int, optional
174 Subplot index if multiple spatial subplots in figure.
176 Returns
177 -------
178 axes:
179 Matplotlib GeoAxes definition.
180 """
181 # Identify min/max plot bounds.
182 try:
183 lat_axis, lon_axis = get_cube_yxcoordname(cube)
184 x1 = np.nanmin(cube.coord(lon_axis).points)
185 x2 = np.nanmax(cube.coord(lon_axis).points)
186 y1 = np.nanmin(cube.coord(lat_axis).points)
187 y2 = np.nanmax(cube.coord(lat_axis).points)
189 # Adjust bounds within +/- 180.0 if x dimension extends beyond half-globe.
190 if np.abs(x2 - x1) > 180.0:
191 x1 = x1 - 180.0
192 x2 = x2 - 180.0
193 logging.debug("Adjusting plot bounds to fit global extent.")
195 # Consider map projection orientation.
196 # Adapting orientation enables plotting across international dateline.
197 # Users can adapt the default central_longitude if alternative projections views.
198 if x2 > 180.0 or x1 < -180.0:
199 central_longitude = 180.0
200 else:
201 central_longitude = 0.0
203 # Define spatial map projection.
204 coord_system = cube.coord(lat_axis).coord_system
205 if isinstance(coord_system, iris.coord_systems.RotatedGeogCS):
206 # Define rotated pole map projection for rotated pole inputs.
207 projection = ccrs.RotatedPole(
208 pole_longitude=coord_system.grid_north_pole_longitude,
209 pole_latitude=coord_system.grid_north_pole_latitude,
210 central_rotated_longitude=central_longitude,
211 )
212 crs = projection
213 elif isinstance(coord_system, iris.coord_systems.TransverseMercator): 213 ↛ 215line 213 didn't jump to line 215 because the condition on line 213 was never true
214 # Define Transverse Mercator projection for TM inputs.
215 projection = ccrs.TransverseMercator(
216 central_longitude=coord_system.longitude_of_central_meridian,
217 central_latitude=coord_system.latitude_of_projection_origin,
218 false_easting=coord_system.false_easting,
219 false_northing=coord_system.false_northing,
220 scale_factor=coord_system.scale_factor_at_central_meridian,
221 )
222 crs = projection
223 else:
224 # Assume polar projection for regional grids encompassing N. Pole
225 if y1 > 20.0 and y2 > 80.0: 225 ↛ 226line 225 didn't jump to line 226 because the condition on line 225 was never true
226 projection = ccrs.NorthPolarStereo(central_longitude=0.0)
227 crs = ccrs.PlateCarree()
228 elif y1 < -80.0 and y2 < -20.0: 228 ↛ 229line 228 didn't jump to line 229 because the condition on line 228 was never true
229 projection = ccrs.SouthPolarStereo(central_longitude=0.0)
230 crs = ccrs.PlateCarree()
231 # Define regular map projection for non-rotated pole inputs.
232 # Alternatives might include e.g. for global model outputs:
233 # projection=ccrs.Robinson(central_longitude=X.y, globe=None)
234 # See also https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html.
235 else:
236 projection = ccrs.PlateCarree(central_longitude=central_longitude)
237 crs = ccrs.PlateCarree()
239 # Define axes for plot (or subplot) with required map projection.
240 if subplot is not None:
241 axes = figure.add_subplot(
242 grid_size[0], grid_size[1], subplot, projection=projection
243 )
244 else:
245 axes = figure.add_subplot(projection=projection)
247 # Add coastlines and borderlines if cube contains x and y map coordinates.
248 # Avoid adding lines for 2D masked data or specific fixed ancillary spatial plots.
249 if (cube.ndim > 1 and iris.util.is_masked(cube.data)) or any( 249 ↛ 252line 249 didn't jump to line 252 because the condition on line 249 was never true
250 name in cube.name() for name in ["land_", "orography", "altitude"]
251 ):
252 pass
253 else:
254 if cmap.name in ["viridis", "Greys"]:
255 coastcol = "magenta"
256 else:
257 coastcol = "black"
258 logging.debug("Plotting coastlines and borderlines in colour %s.", coastcol)
259 axes.coastlines(resolution="10m", color=coastcol)
260 axes.add_feature(cfeature.BORDERS, edgecolor=coastcol)
262 # Add gridlines.
263 gl = axes.gridlines(
264 alpha=0.3,
265 draw_labels=True,
266 dms=False,
267 x_inline=False,
268 y_inline=False,
269 )
270 gl.top_labels = False
271 gl.right_labels = False
272 if subplot:
273 gl.bottom_labels = False
274 gl.left_labels = False
275 if subplot % grid_size[1] == 1:
276 gl.left_labels = True
277 if subplot > ((grid_size[0] - 1) * grid_size[1]): 277 ↛ 282line 277 didn't jump to line 282 because the condition on line 277 was always true
278 gl.bottom_labels = True
280 # If is lat/lon spatial map, fix extent to keep plot tight.
281 # Specifying crs within set_extent helps ensure only data region is shown.
282 if isinstance(coord_system, iris.coord_systems.GeogCS):
283 axes.set_extent([x1, x2, y1, y2], crs=crs)
285 except ValueError:
286 # Skip if not both x and y map coordinates.
287 axes = figure.gca()
288 pass
290 return axes
293def _get_plot_resolution() -> int:
294 """Get resolution of rasterised plots in pixels per inch."""
295 return get_recipe_metadata().get("plot_resolution", 100)
298def _get_start_end_strings(seq_coord: iris.coords.Coord, use_bounds: bool):
299 """Return title and filename based on start and end points or bounds."""
300 if use_bounds and seq_coord.has_bounds():
301 vals = seq_coord.bounds.flatten()
302 else:
303 vals = seq_coord.points
304 start = seq_coord.units.title(vals[0])
305 end = seq_coord.units.title(vals[-1])
307 if start == end:
308 sequence_title = f"\n [{start}]"
309 sequence_fname = f"_{filename_slugify(start)}"
310 else:
311 sequence_title = f"\n [{start} to {end}]"
312 sequence_fname = f"_{filename_slugify(start)}_{filename_slugify(end)}"
314 # Do not include time if coord set to zero.
315 if (
316 seq_coord.units == "hours since 0001-01-01 00:00:00"
317 and vals[0] == 0
318 and vals[-1] == 0
319 ):
320 sequence_title = ""
321 sequence_fname = ""
323 return sequence_title, sequence_fname
326def _set_title_and_filename(
327 seq_coord: iris.coords.Coord,
328 nplot: int,
329 recipe_title: str,
330 filename: str,
331):
332 """Set plot title and filename based on cube coordinate.
334 Parameters
335 ----------
336 sequence_coordinate: iris.coords.Coord
337 Coordinate about which to make a plot sequence.
338 nplot: int
339 Number of output plots to generate - controls title/naming.
340 recipe_title: str
341 Default plot title, potentially to update.
342 filename: str
343 Input plot filename, potentially to update.
345 Returns
346 -------
347 plot_title: str
348 Output formatted plot title string, based on plotted data.
349 plot_filename: str
350 Output formatted plot filename string.
351 """
352 ndim = seq_coord.ndim
353 npoints = np.size(seq_coord.points)
354 sequence_title = ""
355 sequence_fname = ""
357 # Case 1: Multiple dimension sequence input - list number of aggregated cases
358 # (e.g. aggregation histogram plots)
359 if ndim > 1:
360 ncase = np.shape(seq_coord)[0]
361 sequence_title = f"\n [{ncase} cases]"
362 sequence_fname = f"_{ncase}cases"
364 # Case 2: Single dimension input
365 else:
366 # Single sequence point
367 if npoints == 1:
368 if nplot > 1:
369 # Default labels for sequence inputs
370 sequence_value = seq_coord.units.title(seq_coord.points[0])
371 sequence_value = sequence_value.replace(" unknown", "")
372 sequence_title = f"\n [{sequence_value}]"
373 sequence_fname = f"_{filename_slugify(sequence_value)}"
374 else:
375 # Aggregated attribute available where input collapsed over aggregation
376 try:
377 ncase = seq_coord.attributes["number_reference_times"]
378 sequence_title = f"\n [{ncase} cases]"
379 sequence_fname = f"_{ncase}cases"
380 except KeyError:
381 sequence_title, sequence_fname = _get_start_end_strings(
382 seq_coord, use_bounds=seq_coord.has_bounds()
383 )
384 # Multiple sequence (e.g. time) points
385 else:
386 sequence_title, sequence_fname = _get_start_end_strings(
387 seq_coord, use_bounds=False
388 )
390 # Set plot title and filename
391 plot_title = f"{recipe_title}{sequence_title}"
393 # Set plot filename, defaulting to user input if provided.
394 if filename is None:
395 filename = slugify(recipe_title)
396 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
397 else:
398 if nplot > 1:
399 plot_filename = f"{filename.rsplit('.', 1)[0]}{sequence_fname}.png"
400 else:
401 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
403 return plot_title, plot_filename
406def _set_postage_stamp_title(stamp_coord: iris.coords.Coord) -> str:
407 """Control postage stamp plot output titles based on stamp coordinate."""
408 if stamp_coord.name() == "realization":
409 mtitle = "Member"
410 else:
411 mtitle = stamp_coord.name().capitalize()
413 if stamp_coord.name() == "time":
414 mtitle = f"{stamp_coord.units.title(stamp_coord.points[0])}"
415 else:
416 mtitle = f"{mtitle} #{stamp_coord.points[0]}"
418 return mtitle
421def _plot_and_save_spatial_plot(
422 cube: iris.cube.Cube,
423 filename: str,
424 title: str,
425 method: Literal["contourf", "pcolormesh"],
426 overlay_cube: iris.cube.Cube | None = None,
427 contour_cube: iris.cube.Cube | None = None,
428 scatter_cube: iris.cube.Cube | None = None,
429 **kwargs,
430):
431 """Plot and save a spatial plot.
433 Parameters
434 ----------
435 cube: Cube
436 2 dimensional (lat and lon) Cube of the data to plot.
437 filename: str
438 Filename of the plot to write.
439 title: str
440 Plot title.
441 method: "contourf" | "pcolormesh" | "scatter"
442 The plotting method to use.
443 overlay_cube: Cube, optional
444 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
445 contour_cube: Cube, optional
446 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
447 scatter_cube: Cube, optional
448 Optional 1 (station list) or 2 dimensional (lat and lon) Cube of data to overplot as map of scatter points over base cube
449 """
450 # Setup plot details, size, resolution, etc.
451 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
453 # Specify the color bar
454 cmap, levels, norm = _colorbar_map_levels(cube)
456 # If overplotting, set required colorbars
457 if overlay_cube:
458 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
459 if contour_cube:
460 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
462 # Setup plot map projection, extent and coastlines and borderlines.
463 axes = _setup_spatial_map(cube, fig, cmap)
465 # Plot the field.
466 try:
467 vmin = min(levels)
468 vmax = max(levels)
469 except TypeError:
470 vmin, vmax = None, None
471 # Ensure to use norm and not vmin/vmax if levels are defined.
472 if norm is not None:
473 vmin = None
474 vmax = None
475 logging.debug("Plotting using defined levels.")
476 if method == "contourf":
477 # Filled contour plot of the field.
478 plot = iplt.contourf(cube, cmap=cmap, levels=levels, norm=norm)
479 elif method == "pcolormesh":
480 plot = iplt.pcolormesh(cube, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
481 elif method == "scatter": 481 ↛ 489line 481 didn't jump to line 489 because the condition on line 481 was never true
482 # Scatter plot of the field. The marker size is chosen to give
483 # symbols that decrease in size as the number of observations
484 # increases, although the fraction of the figure covered by
485 # symbols increases roughly as N^(1/2), disregarding overlaps,
486 # and has been selected for the default figure size of (10, 10).
487 # Should this be changed, the marker size should be adjusted in
488 # proportion to the area of the figure.
489 mrk_size = int(np.sqrt(2500000.0 / len(cube.data)))
490 lat_axis, lon_axis = get_cube_yxcoordname(cube)
491 plot = iplt.scatter(
492 cube.coord(lon_axis),
493 cube.coord(lat_axis),
494 c=cube.data[:],
495 s=mrk_size,
496 cmap=cmap,
497 edgecolors="k",
498 norm=norm,
499 vmin=vmin,
500 vmax=vmax,
501 )
502 else:
503 raise ValueError(f"Unknown plotting method: {method}")
505 # Overplot overlay field, if required
506 if overlay_cube:
507 try:
508 over_vmin = min(over_levels)
509 over_vmax = max(over_levels)
510 except TypeError:
511 over_vmin, over_vmax = None, None
512 if over_norm is not None: 512 ↛ 513line 512 didn't jump to line 513 because the condition on line 512 was never true
513 over_vmin = None
514 over_vmax = None
515 overlay = iplt.pcolormesh(
516 overlay_cube,
517 cmap=over_cmap,
518 norm=over_norm,
519 alpha=0.8,
520 vmin=over_vmin,
521 vmax=over_vmax,
522 )
523 # Overplot contour field, if required, with contour labelling.
524 if contour_cube:
525 contour = iplt.contour(
526 contour_cube,
527 colors="darkgray",
528 levels=cntr_levels,
529 norm=cntr_norm,
530 alpha=0.5,
531 linestyles="--",
532 linewidths=1,
533 )
534 plt.clabel(contour)
535 # Overplot valid elements of scatter field, if required
536 if scatter_cube: 536 ↛ 537line 536 didn't jump to line 537 because the condition on line 536 was never true
537 mrk_size = int(np.sqrt(2500000.0 / len(scatter_cube.data)))
538 lat_axis, lon_axis = get_cube_yxcoordname(scatter_cube)
539 lon_coord = scatter_cube.coord(lon_axis)
540 lat_coord = scatter_cube.coord(lat_axis)
541 valid = ~scatter_cube.data.mask
542 valid_lon = iris.coords.AuxCoord(
543 lon_coord.points[valid],
544 standard_name=lon_coord.standard_name,
545 units=lon_coord.units,
546 coord_system=lon_coord.coord_system,
547 )
548 valid_lat = iris.coords.AuxCoord(
549 lat_coord.points[valid],
550 standard_name=lat_coord.standard_name,
551 units=lat_coord.units,
552 coord_system=lat_coord.coord_system,
553 )
554 iplt.scatter(
555 valid_lon,
556 valid_lat,
557 c=scatter_cube.data[valid],
558 s=mrk_size,
559 cmap=cmap,
560 edgecolors="k",
561 norm=norm,
562 vmin=vmin,
563 vmax=vmax,
564 )
566 # Check to see if transect, and if so, adjust y axis.
567 if is_transect(cube):
568 if "pressure" in [coord.name() for coord in cube.coords()]:
569 axes.invert_yaxis()
570 axes.set_yscale("log")
571 axes.set_ylim(1100, 100)
572 # If both model_level_number and level_height exists, iplt can construct
573 # plot as a function of height above orography (NOT sea level).
574 elif {"model_level_number", "level_height"}.issubset( 574 ↛ 579line 574 didn't jump to line 579 because the condition on line 574 was always true
575 {coord.name() for coord in cube.coords()}
576 ):
577 axes.set_yscale("log")
579 axes.set_title(
580 f"{title}\n"
581 f"Start Lat: {cube.attributes['transect_coords'].split('_')[0]}"
582 f" Start Lon: {cube.attributes['transect_coords'].split('_')[1]}"
583 f" End Lat: {cube.attributes['transect_coords'].split('_')[2]}"
584 f" End Lon: {cube.attributes['transect_coords'].split('_')[3]}",
585 fontsize=16,
586 )
588 # Inset code
589 axins = inset_axes(
590 axes,
591 width="20%",
592 height="20%",
593 loc="upper right",
594 axes_class=GeoAxes,
595 axes_kwargs={"map_projection": ccrs.PlateCarree()},
596 )
598 # Slightly transparent to reduce plot blocking.
599 axins.patch.set_alpha(0.4)
601 axins.coastlines(resolution="50m")
602 axins.add_feature(cfeature.BORDERS, linewidth=0.3)
604 SLat, SLon, ELat, ELon = (
605 float(coord) for coord in cube.attributes["transect_coords"].split("_")
606 )
608 # Draw line between them
609 axins.plot(
610 [SLon, ELon], [SLat, ELat], color="black", transform=ccrs.PlateCarree()
611 )
613 # Plot points (note: lon, lat order for Cartopy)
614 axins.plot(SLon, SLat, marker="x", color="green", transform=ccrs.PlateCarree())
615 axins.plot(ELon, ELat, marker="x", color="red", transform=ccrs.PlateCarree())
617 lon_min, lon_max = sorted([SLon, ELon])
618 lat_min, lat_max = sorted([SLat, ELat])
620 # Midpoints
621 lon_mid = (lon_min + lon_max) / 2
622 lat_mid = (lat_min + lat_max) / 2
624 # Maximum half-range
625 half_range = max(lon_max - lon_min, lat_max - lat_min) / 2
626 if half_range == 0: # points identical → provide small default 626 ↛ 630line 626 didn't jump to line 630 because the condition on line 626 was always true
627 half_range = 1
629 # Set square extent
630 axins.set_extent(
631 [
632 lon_mid - half_range,
633 lon_mid + half_range,
634 lat_mid - half_range,
635 lat_mid + half_range,
636 ],
637 crs=ccrs.PlateCarree(),
638 )
640 # Ensure square aspect
641 axins.set_aspect("equal")
643 else:
644 # Add title.
645 axes.set_title(title, fontsize=16)
647 # Adjust padding if spatial plot or transect
648 if is_transect(cube):
649 yinfopad = -0.1
650 ycbarpad = 0.1
651 else:
652 yinfopad = 0.01
653 ycbarpad = 0.042
655 # Add watermark with min/max/mean. Currently not user togglable.
656 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
657 axes.annotate(
658 f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
659 xy=(0.025, yinfopad),
660 xycoords="axes fraction",
661 xytext=(-5, 5),
662 textcoords="offset points",
663 ha="left",
664 va="bottom",
665 size=11,
666 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
667 )
669 # Add secondary colour bar for overlay_cube field if required.
670 if overlay_cube:
671 cbarB = fig.colorbar(
672 overlay, orientation="horizontal", location="bottom", pad=0.0, shrink=0.7
673 )
674 cbarB.set_label(label=f"{overlay_cube.name()} ({overlay_cube.units})", size=14)
675 # add ticks and tick_labels for every levels if less than 20 levels exist
676 if over_levels is not None and len(over_levels) < 20: 676 ↛ 677line 676 didn't jump to line 677 because the condition on line 676 was never true
677 cbarB.set_ticks(over_levels)
678 cbarB.set_ticklabels([f"{level:.2f}" for level in over_levels])
679 if "rainfall" or "snowfall" or "visibility" in overlay_cube.name():
680 cbarB.set_ticklabels([f"{level:.3g}" for level in over_levels])
681 logging.debug("Set secondary colorbar ticks and labels.")
683 # Add main colour bar.
684 cbar = fig.colorbar(
685 plot, orientation="horizontal", location="bottom", pad=ycbarpad, shrink=0.7
686 )
688 cbar.set_label(label=f"{cube.name()} ({cube.units})", size=14)
689 # add ticks and tick_labels for every levels if less than 20 levels exist
690 if levels is not None and len(levels) < 20:
691 cbar.set_ticks(levels)
692 cbar.set_ticklabels([f"{level:.2f}" for level in levels])
693 if "rainfall" or "snowfall" or "visibility" in cube.name(): 693 ↛ 695line 693 didn't jump to line 695 because the condition on line 693 was always true
694 cbar.set_ticklabels([f"{level:.3g}" for level in levels])
695 logging.debug("Set colorbar ticks and labels.")
697 # Save plot.
698 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
699 logging.info("Saved spatial plot to %s", filename)
700 plt.close(fig)
703def _plot_and_save_postage_stamp_spatial_plot(
704 cube: iris.cube.Cube,
705 filename: str,
706 stamp_coordinate: str,
707 title: str,
708 method: Literal["contourf", "pcolormesh"],
709 overlay_cube: iris.cube.Cube | None = None,
710 contour_cube: iris.cube.Cube | None = None,
711 **kwargs,
712):
713 """Plot postage stamp spatial plots from an ensemble.
715 Parameters
716 ----------
717 cube: Cube
718 Iris cube of data to be plotted. It must have the stamp coordinate.
719 filename: str
720 Filename of the plot to write.
721 stamp_coordinate: str
722 Coordinate that becomes different plots.
723 method: "contourf" | "pcolormesh"
724 The plotting method to use.
725 overlay_cube: Cube, optional
726 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
727 contour_cube: Cube, optional
728 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
730 Raises
731 ------
732 ValueError
733 If the cube doesn't have the right dimensions.
734 """
735 # Use the smallest square grid that will fit the members.
736 nmember = len(cube.coord(stamp_coordinate).points)
737 grid_rows = int(math.sqrt(nmember))
738 grid_size = math.ceil(nmember / grid_rows)
740 fig = plt.figure(
741 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
742 )
744 # Specify the color bar
745 cmap, levels, norm = _colorbar_map_levels(cube)
746 # If overplotting, set required colorbars
747 if overlay_cube: 747 ↛ 748line 747 didn't jump to line 748 because the condition on line 747 was never true
748 over_cmap, over_levels, over_norm = _colorbar_map_levels(overlay_cube)
749 if contour_cube: 749 ↛ 750line 749 didn't jump to line 750 because the condition on line 749 was never true
750 cntr_cmap, cntr_levels, cntr_norm = _colorbar_map_levels(contour_cube)
752 # Make a subplot for each member.
753 for member, subplot in zip(
754 cube.slices_over(stamp_coordinate),
755 range(1, grid_size * grid_rows + 1),
756 strict=False,
757 ):
758 # Setup subplot map projection, extent and coastlines and borderlines.
759 axes = _setup_spatial_map(
760 member, fig, cmap, grid_size=(grid_rows, grid_size), subplot=subplot
761 )
762 if method == "contourf":
763 # Filled contour plot of the field.
764 plot = iplt.contourf(member, cmap=cmap, levels=levels, norm=norm)
765 elif method == "pcolormesh":
766 if levels is not None:
767 vmin = min(levels)
768 vmax = max(levels)
769 else:
770 raise TypeError("Unknown vmin and vmax range.")
771 vmin, vmax = None, None
772 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
773 # if levels are defined.
774 if norm is not None: 774 ↛ 775line 774 didn't jump to line 775 because the condition on line 774 was never true
775 vmin = None
776 vmax = None
777 # pcolormesh plot of the field.
778 plot = iplt.pcolormesh(member, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
779 else:
780 raise ValueError(f"Unknown plotting method: {method}")
782 # Overplot overlay field, if required
783 if overlay_cube: 783 ↛ 784line 783 didn't jump to line 784 because the condition on line 783 was never true
784 try:
785 over_vmin = min(over_levels)
786 over_vmax = max(over_levels)
787 except TypeError:
788 over_vmin, over_vmax = None, None
789 if over_norm is not None:
790 over_vmin = None
791 over_vmax = None
792 iplt.pcolormesh(
793 overlay_cube[member.coord(stamp_coordinate).points[0]],
794 cmap=over_cmap,
795 norm=over_norm,
796 alpha=0.6,
797 vmin=over_vmin,
798 vmax=over_vmax,
799 )
800 # Overplot contour field, if required
801 if contour_cube: 801 ↛ 802line 801 didn't jump to line 802 because the condition on line 801 was never true
802 iplt.contour(
803 contour_cube[member.coord(stamp_coordinate).points[0]],
804 colors="darkgray",
805 levels=cntr_levels,
806 norm=cntr_norm,
807 alpha=0.6,
808 linestyles="--",
809 linewidths=1,
810 )
811 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
812 axes.set_title(f"{mtitle}")
814 # Put the shared colorbar in its own axes.
815 colorbar_axes = fig.add_axes([0.15, 0.05, 0.7, 0.03])
816 colorbar = fig.colorbar(
817 plot, colorbar_axes, orientation="horizontal", pad=0.042, shrink=0.7
818 )
819 colorbar.set_label(f"{cube.name()} ({cube.units})", size=14)
821 # Overall figure title.
822 fig.suptitle(title, fontsize=16)
824 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
825 logging.info("Saved contour postage stamp plot to %s", filename)
826 plt.close(fig)
829def _plot_and_save_line_series(
830 cubes: iris.cube.CubeList,
831 coords: list[iris.coords.Coord],
832 ensemble_coord: str,
833 filename: str,
834 title: str,
835 **kwargs,
836):
837 """Plot and save a 1D line series.
839 Parameters
840 ----------
841 cubes: Cube or CubeList
842 Cube or CubeList containing the cubes to plot on the y-axis.
843 coords: list[Coord]
844 Coordinates to plot on the x-axis, one per cube.
845 ensemble_coord: str
846 Ensemble coordinate in the cube.
847 filename: str
848 Filename of the plot to write.
849 title: str
850 Plot title.
851 """
852 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
854 model_colors_map = _get_model_colors_map(cubes)
856 # Store min/max ranges.
857 y_levels = []
859 # Check match-up across sequence coords gives consistent sizes
860 _validate_cubes_coords(cubes, coords)
862 for cube, coord in zip(cubes, coords, strict=True):
863 label = None
864 color = "black"
865 if model_colors_map:
866 label = cube.attributes.get("model_name")
867 color = model_colors_map.get(label)
868 for cube_slice in cube.slices_over(ensemble_coord):
869 # Label with (control) if part of an ensemble or not otherwise.
870 if cube_slice.coord(ensemble_coord).points == [0]:
871 iplt.plot(
872 coord,
873 cube_slice,
874 color=color,
875 marker="o",
876 ls="-",
877 lw=3,
878 label=f"{label} (control)"
879 if len(cube.coord(ensemble_coord).points) > 1
880 else label,
881 )
882 # Label with (perturbed) if part of an ensemble and not the control.
883 else:
884 iplt.plot(
885 coord,
886 cube_slice,
887 color=color,
888 ls="-",
889 lw=1.5,
890 alpha=0.75,
891 label=f"{label} (member)",
892 )
894 # Calculate the global min/max if multiple cubes are given.
895 _, levels, _ = _colorbar_map_levels(cube, axis="y")
896 if levels is not None: 896 ↛ 897line 896 didn't jump to line 897 because the condition on line 896 was never true
897 y_levels.append(min(levels))
898 y_levels.append(max(levels))
900 # Get the current axes.
901 ax = plt.gca()
903 # Add some labels and tweak the style.
904 # check if cubes[0] works for single cube if not CubeList
905 if coords[0].name() == "time":
906 ax.set_xlabel(f"{coords[0].name()}", fontsize=14)
907 else:
908 ax.set_xlabel(f"{coords[0].name()} / {coords[0].units}", fontsize=14)
909 ax.set_ylabel(f"{cubes[0].name()} / {cubes[0].units}", fontsize=14)
910 ax.set_title(title, fontsize=16)
912 ax.ticklabel_format(axis="y", useOffset=False)
913 ax.tick_params(axis="x", labelrotation=15)
914 ax.tick_params(axis="both", labelsize=12)
916 # Set y limits to global min and max, autoscale if colorbar doesn't exist.
917 if y_levels: 917 ↛ 918line 917 didn't jump to line 918 because the condition on line 917 was never true
918 ax.set_ylim(min(y_levels), max(y_levels))
919 # Add zero line.
920 if min(y_levels) < 0.0 and max(y_levels) > 0.0:
921 ax.axhline(y=0, xmin=0, xmax=1, ls="-", color="grey", lw=2)
922 logging.debug(
923 "Line plot with y-axis limits %s-%s", min(y_levels), max(y_levels)
924 )
925 else:
926 ax.autoscale()
928 # Add gridlines
929 ax.grid(linestyle="--", color="grey", linewidth=1)
930 # Ientify unique labels for legend
931 handles = list(
932 {
933 label: handle
934 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
935 }.values()
936 )
937 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
939 # Save plot.
940 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
941 logging.info("Saved line plot to %s", filename)
942 plt.close(fig)
945def _plot_and_save_vertical_line_series(
946 cubes: iris.cube.CubeList,
947 coords: list[iris.coords.Coord],
948 ensemble_coord: str,
949 filename: str,
950 series_coordinate: str,
951 title: str,
952 vmin: float,
953 vmax: float,
954 **kwargs,
955):
956 """Plot and save a 1D line series in vertical.
958 Parameters
959 ----------
960 cubes: CubeList
961 1 dimensional Cube or CubeList of the data to plot on x-axis.
962 coord: list[Coord]
963 Coordinates to plot on the y-axis, one per cube.
964 ensemble_coord: str
965 Ensemble coordinate in the cube.
966 filename: str
967 Filename of the plot to write.
968 series_coordinate: str
969 Coordinate to use as vertical axis.
970 title: str
971 Plot title.
972 vmin: float
973 Minimum value for the x-axis.
974 vmax: float
975 Maximum value for the x-axis.
976 """
977 # plot the vertical pressure axis using log scale
978 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
980 model_colors_map = _get_model_colors_map(cubes)
982 # Check match-up across sequence coords gives consistent sizes
983 _validate_cubes_coords(cubes, coords)
985 for cube, coord in zip(cubes, coords, strict=True):
986 label = None
987 color = "black"
988 if model_colors_map: 988 ↛ 989line 988 didn't jump to line 989 because the condition on line 988 was never true
989 label = cube.attributes.get("model_name")
990 color = model_colors_map.get(label)
992 for cube_slice in cube.slices_over(ensemble_coord):
993 # If ensemble data given plot control member with (control)
994 # unless single forecast.
995 if cube_slice.coord(ensemble_coord).points == [0]:
996 iplt.plot(
997 cube_slice,
998 coord,
999 color=color,
1000 marker="o",
1001 ls="-",
1002 lw=3,
1003 label=f"{label} (control)"
1004 if len(cube.coord(ensemble_coord).points) > 1
1005 else label,
1006 )
1007 # If ensemble data given plot perturbed members with (perturbed).
1008 else:
1009 iplt.plot(
1010 cube_slice,
1011 coord,
1012 color=color,
1013 ls="-",
1014 lw=1.5,
1015 alpha=0.75,
1016 label=f"{label} (member)",
1017 )
1019 # Get the current axis
1020 ax = plt.gca()
1022 # Special handling for pressure level data.
1023 if series_coordinate == "pressure": 1023 ↛ 1045line 1023 didn't jump to line 1045 because the condition on line 1023 was always true
1024 # Invert y-axis and set to log scale.
1025 ax.invert_yaxis()
1026 ax.set_yscale("log")
1028 # Define y-ticks and labels for pressure log axis.
1029 y_tick_labels = [
1030 "1000",
1031 "850",
1032 "700",
1033 "500",
1034 "300",
1035 "200",
1036 "100",
1037 ]
1038 y_ticks = [1000, 850, 700, 500, 300, 200, 100]
1040 # Set y-axis limits and ticks.
1041 ax.set_ylim(1100, 100)
1043 # Test if series_coordinate is model level data. The UM data uses
1044 # model_level_number and lfric uses full_levels as coordinate.
1045 elif series_coordinate in ("model_level_number", "full_levels", "half_levels"):
1046 # Define y-ticks and labels for vertical axis.
1047 y_ticks = iter_maybe(cubes)[0].coord(series_coordinate).points
1048 y_tick_labels = [str(int(i)) for i in y_ticks]
1049 ax.set_ylim(min(y_ticks), max(y_ticks))
1051 ax.set_yticks(y_ticks)
1052 ax.set_yticklabels(y_tick_labels)
1054 # Set x-axis limits.
1055 ax.set_xlim(vmin, vmax)
1056 # Mark y=0 if present in plot.
1057 if vmin < 0.0 and vmax > 0.0: 1057 ↛ 1058line 1057 didn't jump to line 1058 because the condition on line 1057 was never true
1058 ax.axvline(x=0, ymin=0, ymax=1, ls="-", color="grey", lw=2)
1060 # Add some labels and tweak the style.
1061 ax.set_ylabel(f"{coord.name()} / {coord.units}", fontsize=14)
1062 ax.set_xlabel(
1063 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1064 )
1065 ax.set_title(title, fontsize=16)
1066 ax.ticklabel_format(axis="x")
1067 ax.tick_params(axis="y")
1068 ax.tick_params(axis="both", labelsize=12)
1070 # Add gridlines
1071 ax.grid(linestyle="--", color="grey", linewidth=1)
1072 # Ientify unique labels for legend
1073 handles = list(
1074 {
1075 label: handle
1076 for (handle, label) in zip(*ax.get_legend_handles_labels(), strict=True)
1077 }.values()
1078 )
1079 ax.legend(handles=handles, loc="best", ncol=1, frameon=True, fontsize=16)
1081 # Save plot.
1082 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1083 logging.info("Saved line plot to %s", filename)
1084 plt.close(fig)
1087def _plot_and_save_scatter_plot(
1088 cube_x: iris.cube.Cube | iris.cube.CubeList,
1089 cube_y: iris.cube.Cube | iris.cube.CubeList,
1090 filename: str,
1091 title: str,
1092 one_to_one: bool,
1093 model_names: list[str] = None,
1094 **kwargs,
1095):
1096 """Plot and save a 2D scatter plot.
1098 Parameters
1099 ----------
1100 cube_x: Cube | CubeList
1101 1 dimensional Cube or CubeList of the data to plot on x-axis.
1102 cube_y: Cube | CubeList
1103 1 dimensional Cube or CubeList of the data to plot on y-axis.
1104 filename: str
1105 Filename of the plot to write.
1106 title: str
1107 Plot title.
1108 one_to_one: bool
1109 Whether a 1:1 line is plotted.
1110 """
1111 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1112 # plot the cube_x and cube_y 1D fields as a scatter plot. If they are CubeLists this ensures
1113 # to pair each cube from cube_x with the corresponding cube from cube_y, allowing to iterate
1114 # over the pairs simultaneously.
1116 # Ensure cube_x and cube_y are iterable
1117 cube_x_iterable = iter_maybe(cube_x)
1118 cube_y_iterable = iter_maybe(cube_y)
1120 for cube_x_iter, cube_y_iter in zip(cube_x_iterable, cube_y_iterable, strict=True):
1121 iplt.scatter(cube_x_iter, cube_y_iter)
1122 if one_to_one is True:
1123 plt.plot(
1124 [
1125 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1126 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1127 ],
1128 [
1129 np.nanmin([np.nanmin(cube_y.data), np.nanmin(cube_x.data)]),
1130 np.nanmax([np.nanmax(cube_y.data), np.nanmax(cube_x.data)]),
1131 ],
1132 "k",
1133 linestyle="--",
1134 )
1135 ax = plt.gca()
1137 # Add some labels and tweak the style.
1138 if model_names is None:
1139 ax.set_xlabel(f"{cube_x[0].name()} / {cube_x[0].units}", fontsize=14)
1140 ax.set_ylabel(f"{cube_y[0].name()} / {cube_y[0].units}", fontsize=14)
1141 else:
1142 # Add the model names, these should be order of base (x) and other (y).
1143 ax.set_xlabel(
1144 f"{model_names[0]}_{cube_x[0].name()} / {cube_x[0].units}", fontsize=14
1145 )
1146 ax.set_ylabel(
1147 f"{model_names[1]}_{cube_y[0].name()} / {cube_y[0].units}", fontsize=14
1148 )
1149 ax.set_title(title, fontsize=16)
1150 ax.ticklabel_format(axis="y", useOffset=False)
1151 ax.tick_params(axis="x", labelrotation=15)
1152 ax.tick_params(axis="both", labelsize=12)
1153 ax.autoscale()
1155 # Save plot.
1156 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1157 logging.info("Saved scatter plot to %s", filename)
1158 plt.close(fig)
1161def _plot_and_save_vector_plot(
1162 cube_u: iris.cube.Cube,
1163 cube_v: iris.cube.Cube,
1164 filename: str,
1165 title: str,
1166 method: Literal["contourf", "pcolormesh"],
1167 **kwargs,
1168):
1169 """Plot and save a 2D vector plot.
1171 Parameters
1172 ----------
1173 cube_u: Cube
1174 2 dimensional Cube of u component of the data.
1175 cube_v: Cube
1176 2 dimensional Cube of v component of the data.
1177 filename: str
1178 Filename of the plot to write.
1179 title: str
1180 Plot title.
1181 """
1182 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1184 # Create a cube containing the magnitude of the vector field.
1185 cube_vec_mag = (cube_u**2 + cube_v**2) ** 0.5
1186 cube_vec_mag.rename(f"{cube_u.name()}_{cube_v.name()}_magnitude")
1188 # Specify the color bar
1189 cmap, levels, norm = _colorbar_map_levels(cube_vec_mag)
1191 # Setup plot map projection, extent and coastlines and borderlines.
1192 axes = _setup_spatial_map(cube_vec_mag, fig, cmap)
1194 if method == "contourf":
1195 # Filled contour plot of the field.
1196 plot = iplt.contourf(cube_vec_mag, cmap=cmap, levels=levels, norm=norm)
1197 elif method == "pcolormesh":
1198 try:
1199 vmin = min(levels)
1200 vmax = max(levels)
1201 except TypeError:
1202 vmin, vmax = None, None
1203 # pcolormesh plot of the field and ensure to use norm and not vmin/vmax
1204 # if levels are defined.
1205 if norm is not None:
1206 vmin = None
1207 vmax = None
1208 plot = iplt.pcolormesh(cube_vec_mag, cmap=cmap, norm=norm, vmin=vmin, vmax=vmax)
1209 else:
1210 raise ValueError(f"Unknown plotting method: {method}")
1212 # Check to see if transect, and if so, adjust y axis.
1213 if is_transect(cube_vec_mag):
1214 if "pressure" in [coord.name() for coord in cube_vec_mag.coords()]:
1215 axes.invert_yaxis()
1216 axes.set_yscale("log")
1217 axes.set_ylim(1100, 100)
1218 # If both model_level_number and level_height exists, iplt can construct
1219 # plot as a function of height above orography (NOT sea level).
1220 elif {"model_level_number", "level_height"}.issubset(
1221 {coord.name() for coord in cube_vec_mag.coords()}
1222 ):
1223 axes.set_yscale("log")
1225 axes.set_title(
1226 f"{title}\n"
1227 f"Start Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[0]}"
1228 f" Start Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[1]}"
1229 f" End Lat: {cube_vec_mag.attributes['transect_coords'].split('_')[2]}"
1230 f" End Lon: {cube_vec_mag.attributes['transect_coords'].split('_')[3]}",
1231 fontsize=16,
1232 )
1234 else:
1235 # Add title.
1236 axes.set_title(title, fontsize=16)
1238 # Add watermark with min/max/mean. Currently not user togglable.
1239 # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
1240 axes.annotate(
1241 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}",
1242 xy=(0.05, -0.05),
1243 xycoords="axes fraction",
1244 xytext=(-5, 5),
1245 textcoords="offset points",
1246 ha="right",
1247 va="bottom",
1248 size=11,
1249 bbox=dict(boxstyle="round", fc="#cccccc", ec="#808080", alpha=0.9),
1250 )
1252 # Add colour bar.
1253 cbar = fig.colorbar(plot, orientation="horizontal", pad=0.042, shrink=0.7)
1254 cbar.set_label(label=f"{cube_vec_mag.name()} ({cube_vec_mag.units})", size=14)
1255 # add ticks and tick_labels for every levels if less than 20 levels exist
1256 if levels is not None and len(levels) < 20:
1257 cbar.set_ticks(levels)
1258 cbar.set_ticklabels([f"{level:.1f}" for level in levels])
1260 # 30 barbs along the longest axis of the plot, or a barb per point for data
1261 # with less than 30 points.
1262 step = max(max(cube_u.shape) // 30, 1)
1263 iplt.quiver(cube_u[::step, ::step], cube_v[::step, ::step], pivot="middle")
1265 # Save plot.
1266 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1267 logging.info("Saved vector plot to %s", filename)
1268 plt.close(fig)
1271def _plot_and_save_histogram_series(
1272 cubes: iris.cube.Cube | iris.cube.CubeList,
1273 filename: str,
1274 title: str,
1275 vmin: float,
1276 vmax: float,
1277 **kwargs,
1278):
1279 """Plot and save a histogram series.
1281 Parameters
1282 ----------
1283 cubes: Cube or CubeList
1284 2 dimensional Cube or CubeList of the data to plot as histogram.
1285 filename: str
1286 Filename of the plot to write.
1287 title: str
1288 Plot title.
1289 vmin: float
1290 minimum for colorbar
1291 vmax: float
1292 maximum for colorbar
1293 """
1294 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1295 ax = plt.gca()
1297 model_colors_map = _get_model_colors_map(cubes)
1299 # Set default that histograms will produce probability density function
1300 # at each bin (integral over range sums to 1).
1301 density = True
1303 for cube in iter_maybe(cubes):
1304 # Easier to check title (where var name originates)
1305 # than seeing if long names exist etc.
1306 # Exception case, where distribution better fits log scales/bins.
1307 if "surface_microphysical" in title:
1308 if "amount" in title: 1308 ↛ 1310line 1308 didn't jump to line 1310 because the condition on line 1308 was never true
1309 # Compute histogram following Klingaman et al. (2017): ASoP
1310 bin2 = np.exp(np.log(0.02) + 0.1 * np.linspace(0, 99, 100))
1311 bins = np.pad(bin2, (1, 0), "constant", constant_values=0)
1312 density = False
1313 else:
1314 bins = 10.0 ** (
1315 np.arange(-10, 27, 1) / 10.0
1316 ) # Suggestion from RMED toolbox.
1317 bins = np.insert(bins, 0, 0)
1318 ax.set_yscale("log")
1319 vmin = bins[1]
1320 vmax = bins[-1] # Manually set vmin/vmax to override json derived value.
1321 ax.set_xscale("log")
1322 elif "lightning" in title:
1323 bins = [0, 1, 2, 3, 4, 5]
1324 else:
1325 bins = np.linspace(vmin, vmax, 51)
1326 logging.debug(
1327 "Plotting histogram with %s bins %s - %s.",
1328 np.size(bins),
1329 np.min(bins),
1330 np.max(bins),
1331 )
1333 # Reshape cube data into a single array to allow for a single histogram.
1334 # Otherwise we plot xdim histograms stacked.
1335 cube_data_1d = (cube.data).flatten()
1337 label = None
1338 color = "black"
1339 if model_colors_map: 1339 ↛ 1340line 1339 didn't jump to line 1340 because the condition on line 1339 was never true
1340 label = cube.attributes.get("model_name")
1341 color = model_colors_map[label]
1342 x, y = np.histogram(cube_data_1d, bins=bins, density=density)
1344 # Compute area under curve.
1345 if "surface_microphysical" in title and "amount" in title: 1345 ↛ 1346line 1345 didn't jump to line 1346 because the condition on line 1345 was never true
1346 bin_mean = (bins[:-1] + bins[1:]) / 2.0
1347 x = x * bin_mean / x.sum()
1348 x = x[1:]
1349 y = y[1:]
1351 ax.plot(
1352 y[:-1], x, color=color, linewidth=3, marker="o", markersize=6, label=label
1353 )
1355 # Add some labels and tweak the style.
1356 ax.set_title(title, fontsize=16)
1357 ax.set_xlabel(
1358 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1359 )
1360 ax.set_ylabel("Normalised probability density", fontsize=14)
1361 if "surface_microphysical" in title and "amount" in title: 1361 ↛ 1362line 1361 didn't jump to line 1362 because the condition on line 1361 was never true
1362 ax.set_ylabel(
1363 f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
1364 )
1365 ax.set_xlim(vmin, vmax)
1366 ax.tick_params(axis="both", labelsize=12)
1368 # Overlay grid-lines onto histogram plot.
1369 ax.grid(linestyle="--", color="grey", linewidth=1)
1370 if model_colors_map: 1370 ↛ 1371line 1370 didn't jump to line 1371 because the condition on line 1370 was never true
1371 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1373 # Save plot.
1374 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1375 logging.info("Saved histogram plot to %s", filename)
1376 plt.close(fig)
1379def _plot_and_save_postage_stamp_histogram_series(
1380 cube: iris.cube.Cube,
1381 filename: str,
1382 title: str,
1383 stamp_coordinate: str,
1384 vmin: float,
1385 vmax: float,
1386 **kwargs,
1387):
1388 """Plot and save postage (ensemble members) stamps for a histogram series.
1390 Parameters
1391 ----------
1392 cube: Cube
1393 2 dimensional Cube of the data to plot as histogram.
1394 filename: str
1395 Filename of the plot to write.
1396 title: str
1397 Plot title.
1398 stamp_coordinate: str
1399 Coordinate that becomes different plots.
1400 vmin: float
1401 minimum for pdf x-axis
1402 vmax: float
1403 maximum for pdf x-axis
1404 """
1405 # Use the smallest square grid that will fit the members.
1406 nmember = len(cube.coord(stamp_coordinate).points)
1407 grid_rows = int(math.sqrt(nmember))
1408 grid_size = math.ceil(nmember / grid_rows)
1410 fig = plt.figure(
1411 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1412 )
1413 # Make a subplot for each member.
1414 for member, subplot in zip(
1415 cube.slices_over(stamp_coordinate),
1416 range(1, grid_size * grid_rows + 1),
1417 strict=False,
1418 ):
1419 # Implicit interface is much easier here, due to needing to have the
1420 # cartopy GeoAxes generated.
1421 plt.subplot(grid_rows, grid_size, subplot)
1422 # Reshape cube data into a single array to allow for a single histogram.
1423 # Otherwise we plot xdim histograms stacked.
1424 member_data_1d = (member.data).flatten()
1425 plt.hist(member_data_1d, density=True, stacked=True)
1426 axes = plt.gca()
1427 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1428 axes.set_title(f"{mtitle}")
1429 axes.set_xlim(vmin, vmax)
1431 # Overall figure title.
1432 fig.suptitle(title, fontsize=16)
1434 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1435 logging.info("Saved histogram postage stamp plot to %s", filename)
1436 plt.close(fig)
1439def _plot_and_save_postage_stamps_in_single_plot_histogram_series(
1440 cube: iris.cube.Cube,
1441 filename: str,
1442 title: str,
1443 stamp_coordinate: str,
1444 vmin: float,
1445 vmax: float,
1446 **kwargs,
1447):
1448 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1449 ax.set_title(title, fontsize=16)
1450 ax.set_xlim(vmin, vmax)
1451 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1452 ax.set_ylabel("normalised probability density", fontsize=14)
1453 # Loop over all slices along the stamp_coordinate
1454 for member in cube.slices_over(stamp_coordinate):
1455 # Flatten the member data to 1D
1456 member_data_1d = member.data.flatten()
1457 # Plot the histogram using plt.hist
1458 mtitle = _set_postage_stamp_title(member.coord(stamp_coordinate))
1459 plt.hist(
1460 member_data_1d,
1461 density=True,
1462 stacked=True,
1463 label=f"{mtitle}",
1464 )
1466 # Add a legend
1467 ax.legend(fontsize=16)
1469 # Save the figure to a file
1470 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1471 logging.info("Saved histogram postage stamp plot to %s", filename)
1473 # Close the figure
1474 plt.close(fig)
1477def _plot_and_save_scatter_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 scatter plot series.
1487 Parameters
1488 ----------
1489 cubes: Cube or CubeList
1490 2 dimensional Cube or CubeList of the data to plot as scatter.
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 ## TEST CUBELIST len >= 2##
1507 for cube in iter_maybe(cubes[1:]):
1508 label = None
1509 color = "black"
1510 if model_colors_map:
1511 label = cube.attributes.get("model_name")
1512 color = model_colors_map[label]
1514 iplt.scatter(cubes[0], cube, color=color, marker="o", label=label)
1516 # Add some labels and tweak the style.
1517 ax.set_title(title, fontsize=16)
1518 ax.set_xlabel(
1519 f"{iter_maybe(cubes)[0].name()} / {iter_maybe(cubes)[0].units}", fontsize=14
1520 )
1521 ax.set_ylabel(
1522 f"{iter_maybe(cubes)[1].name()} / {iter_maybe(cubes)[1].units}", fontsize=14
1523 )
1524 ax.tick_params(axis="both", labelsize=12)
1525 ax.autoscale()
1527 lims = [
1528 np.min([ax.get_xlim(), ax.get_ylim()]), # min of both axes
1529 np.max([ax.get_xlim(), ax.get_ylim()]), # max of both axes
1530 ]
1531 ax.plot(lims, lims, "k-", alpha=0.75, zorder=0)
1532 ax.set_aspect("equal")
1533 ax.set_xlim(lims)
1534 ax.set_ylim(lims)
1536 # Overlay grid-lines onto histogram plot.
1537 ax.grid(linestyle="--", color="grey", linewidth=1)
1538 if model_colors_map:
1539 ax.legend(loc="upper left", ncol=1, frameon=True, fontsize=16)
1541 # Save plot.
1542 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1543 logging.info("Saved scatter plot to %s", filename)
1544 plt.close(fig)
1547def _plot_and_save_power_spectrum_series(
1548 cubes: iris.cube.Cube | iris.cube.CubeList,
1549 filename: str,
1550 title: str,
1551 **kwargs,
1552):
1553 """Plot and save a power spectrum series.
1555 Parameters
1556 ----------
1557 cubes: Cube or CubeList
1558 2 dimensional Cube or CubeList of the data to plot as power spectrum.
1559 filename: str
1560 Filename of the plot to write.
1561 title: str
1562 Plot title.
1563 """
1564 fig = plt.figure(figsize=(10, 10), facecolor="w", edgecolor="k")
1565 ax = plt.gca()
1567 model_colors_map = _get_model_colors_map(cubes)
1569 for cube in iter_maybe(cubes):
1570 # Calculate power spectrum
1572 # Extract time coordinate and convert to datetime
1573 time_coord = cube.coord("time")
1574 time_points = time_coord.units.num2date(time_coord.points)
1576 # Choose one time point (e.g., the first one)
1577 target_time = time_points[0]
1579 # Bind target_time inside the lambda using a default argument
1580 time_constraint = iris.Constraint(
1581 time=lambda cell, target_time=target_time: cell.point == target_time
1582 )
1584 cube = cube.extract(time_constraint)
1586 if cube.ndim == 2:
1587 cube_3d = cube.data[np.newaxis, :, :]
1588 logging.debug("Adding in new axis for a 2 dimensional cube.")
1589 elif cube.ndim == 3: 1589 ↛ 1590line 1589 didn't jump to line 1590 because the condition on line 1589 was never true
1590 cube_3d = cube.data
1591 else:
1592 raise ValueError("Cube dimensions unsuitable for power spectra code")
1593 raise ValueError(
1594 f"Cube is {cube.ndim} dimensional. Cube should be 2 or 3 dimensional."
1595 )
1597 # Calculate spectra
1598 ps_array = _DCT_ps(cube_3d)
1600 ps_cube = iris.cube.Cube(
1601 ps_array,
1602 long_name="power_spectra",
1603 )
1605 ps_cube.attributes["model_name"] = cube.attributes.get("model_name")
1607 # Create a frequency/wavelength array for coordinate
1608 ps_len = ps_cube.data.shape[1]
1609 freqs = np.arange(1, ps_len + 1)
1610 freq_coord = iris.coords.DimCoord(freqs, long_name="frequency", units="1")
1612 # Convert datetime to numeric time using original units
1613 numeric_time = time_coord.units.date2num(time_points)
1614 # Create a new DimCoord with numeric time
1615 new_time_coord = iris.coords.DimCoord(
1616 numeric_time, standard_name="time", units=time_coord.units
1617 )
1619 # Add time and frequency coordinate to spectra cube.
1620 ps_cube.add_dim_coord(new_time_coord.copy(), 0)
1621 ps_cube.add_dim_coord(freq_coord.copy(), 1)
1623 # Extract data from the cube
1624 frequency = ps_cube.coord("frequency").points
1625 power_spectrum = ps_cube.data
1627 label = None
1628 color = "black"
1629 if model_colors_map: 1629 ↛ 1630line 1629 didn't jump to line 1630 because the condition on line 1629 was never true
1630 label = ps_cube.attributes.get("model_name")
1631 color = model_colors_map[label]
1632 ax.plot(frequency, power_spectrum[0], color=color, label=label)
1634 # Add some labels and tweak the style.
1635 ax.set_title(title, fontsize=16)
1636 ax.set_xlabel("Wavenumber", fontsize=14)
1637 ax.set_ylabel("Power", fontsize=14)
1638 ax.tick_params(axis="both", labelsize=12)
1640 # Set log-log scale
1641 ax.set_xscale("log")
1642 ax.set_yscale("log")
1644 # Overlay grid-lines onto power spectrum plot.
1645 ax.grid(linestyle="--", color="grey", linewidth=1)
1646 if model_colors_map: 1646 ↛ 1647line 1646 didn't jump to line 1647 because the condition on line 1646 was never true
1647 ax.legend(loc="best", ncol=1, frameon=True, fontsize=16)
1649 # Save plot.
1650 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1651 logging.info("Saved power spectrum plot to %s", filename)
1652 plt.close(fig)
1655def _plot_and_save_postage_stamp_power_spectrum_series(
1656 cube: iris.cube.Cube,
1657 filename: str,
1658 title: str,
1659 stamp_coordinate: str,
1660 **kwargs,
1661):
1662 """Plot and save postage (ensemble members) stamps for a power spectrum series.
1664 Parameters
1665 ----------
1666 cube: Cube
1667 2 dimensional Cube of the data to plot as power spectrum.
1668 filename: str
1669 Filename of the plot to write.
1670 title: str
1671 Plot title.
1672 stamp_coordinate: str
1673 Coordinate that becomes different plots.
1674 """
1675 # Use the smallest square grid that will fit the members.
1676 nmember = len(cube.coord(stamp_coordinate).points)
1677 grid_rows = int(math.sqrt(nmember))
1678 grid_size = math.ceil(nmember / grid_rows)
1680 fig = plt.figure(
1681 figsize=(10, 10 * max(grid_rows / grid_size, 0.5)), facecolor="w", edgecolor="k"
1682 )
1684 # Make a subplot for each member.
1685 for member, subplot in zip(
1686 cube.slices_over(stamp_coordinate),
1687 range(1, grid_size * grid_rows + 1),
1688 strict=False,
1689 ):
1690 # Implicit interface is much easier here, due to needing to have the
1691 # cartopy GeoAxes generated.
1692 plt.subplot(grid_rows, grid_size, subplot)
1694 frequency = member.coord("frequency").points
1696 axes = plt.gca()
1697 axes.plot(frequency, member.data)
1698 axes.set_title(f"Member #{member.coord(stamp_coordinate).points[0]}")
1700 # Overall figure title.
1701 fig.suptitle(title, fontsize=16)
1703 fig.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1704 logging.info("Saved power spectra postage stamp plot to %s", filename)
1705 plt.close(fig)
1708def _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series(
1709 cube: iris.cube.Cube,
1710 filename: str,
1711 title: str,
1712 stamp_coordinate: str,
1713 **kwargs,
1714):
1715 fig, ax = plt.subplots(figsize=(10, 10), facecolor="w", edgecolor="k")
1716 ax.set_title(title, fontsize=16)
1717 ax.set_xlabel(f"{cube.name()} / {cube.units}", fontsize=14)
1718 ax.set_ylabel("Power", fontsize=14)
1719 # Loop over all slices along the stamp_coordinate
1720 for member in cube.slices_over(stamp_coordinate):
1721 frequency = member.coord("frequency").points
1722 ax.plot(
1723 frequency,
1724 member.data,
1725 label=f"Member #{member.coord(stamp_coordinate).points[0]}",
1726 )
1728 # Add a legend
1729 ax.legend(fontsize=16)
1731 # Save the figure to a file
1732 plt.savefig(filename, bbox_inches="tight", dpi=_get_plot_resolution())
1733 logging.info("Saved power spectra plot to %s", filename)
1735 # Close the figure
1736 plt.close(fig)
1739def _spatial_plot(
1740 method: Literal["contourf", "pcolormesh"],
1741 cube: iris.cube.Cube,
1742 filename: str | None,
1743 sequence_coordinate: str,
1744 stamp_coordinate: str,
1745 overlay_cube: iris.cube.Cube | None = None,
1746 contour_cube: iris.cube.Cube | None = None,
1747 scatter_cube: iris.cube.Cube | None = None,
1748 **kwargs,
1749):
1750 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1752 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1753 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1754 is present then postage stamp plots will be produced.
1756 If an overlay_cube and/or contour_cube and/or scatter_cube are specified,
1757 multiple variables or outputs can be overplotted on the same figure.
1759 Parameters
1760 ----------
1761 method: "contourf" | "pcolormesh"
1762 The plotting method to use.
1763 cube: Cube
1764 Iris cube of the data to plot. It should have two spatial dimensions,
1765 such as lat and lon, and may also have a another two dimension to be
1766 plotted sequentially and/or as postage stamp plots.
1767 filename: str | None
1768 Name of the plot to write, used as a prefix for plot sequences. If None
1769 uses the recipe name.
1770 sequence_coordinate: str
1771 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1772 This coordinate must exist in the cube.
1773 stamp_coordinate: str
1774 Coordinate about which to plot postage stamp plots. Defaults to
1775 ``"realization"``.
1776 overlay_cube: Cube | None, optional
1777 Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube
1778 contour_cube: Cube | None, optional
1779 Optional 2 dimensional (lat and lon) Cube of data to overplot as contours over base cube
1780 scatter_cube: Cube | None, optional
1781 Optional 1 dimensional (station) or 2 dimensional (lat and lon) Cube of data to overplot as points on scatter map over base cube
1783 Raises
1784 ------
1785 ValueError
1786 If the cube doesn't have the right dimensions.
1787 TypeError
1788 If the cube isn't a single cube.
1789 """
1790 recipe_title = get_recipe_metadata().get("title", "Untitled")
1792 # Ensure we've got a single cube.
1793 cube = _check_single_cube(cube)
1795 # Check if there is a valid stamp coordinate in cube dimensions.
1796 if stamp_coordinate == "realization": 1796 ↛ 1801line 1796 didn't jump to line 1801 because the condition on line 1796 was always true
1797 stamp_coordinate = check_stamp_coordinate(cube)
1799 # Make postage stamp plots if stamp_coordinate exists and has more than a
1800 # single point.
1801 plotting_func = _plot_and_save_spatial_plot
1802 try:
1803 if cube.coord(stamp_coordinate).shape[0] > 1:
1804 plotting_func = _plot_and_save_postage_stamp_spatial_plot
1805 except iris.exceptions.CoordinateNotFoundError:
1806 pass
1808 # Produce a geographical scatter plot if the data have a
1809 # dimension called observation or model_obs_error
1810 if any( 1810 ↛ 1814line 1810 didn't jump to line 1814 because the condition on line 1810 was never true
1811 crd.var_name == "station" or crd.var_name == "model_obs_error"
1812 for crd in cube.coords()
1813 ):
1814 plotting_func = _plot_and_save_spatial_plot
1815 method = "scatter"
1817 # Must have a sequence coordinate.
1818 try:
1819 cube.coord(sequence_coordinate)
1820 except iris.exceptions.CoordinateNotFoundError as err:
1821 raise ValueError(f"Cube must have a {sequence_coordinate} coordinate.") from err
1823 # Create a plot for each value of the sequence coordinate.
1824 plot_index = []
1825 nplot = np.size(cube.coord(sequence_coordinate).points)
1827 for iseq, cube_slice in enumerate(cube.slices_over(sequence_coordinate)):
1828 # Set plot titles and filename
1829 seq_coord = cube_slice.coord(sequence_coordinate)
1830 plot_title, plot_filename = _set_title_and_filename(
1831 seq_coord, nplot, recipe_title, filename
1832 )
1834 # Extract sequence slice for overlay cubes if required.
1835 overlay_slice = slice_over_maybe(overlay_cube, sequence_coordinate, iseq)
1836 contour_slice = slice_over_maybe(contour_cube, sequence_coordinate, iseq)
1837 scatter_slice = slice_over_maybe(scatter_cube, sequence_coordinate, iseq)
1839 # Do the actual plotting.
1840 plotting_func(
1841 cube_slice,
1842 filename=plot_filename,
1843 stamp_coordinate=stamp_coordinate,
1844 title=plot_title,
1845 method=method,
1846 overlay_cube=overlay_slice,
1847 contour_cube=contour_slice,
1848 scatter_cube=scatter_slice,
1849 **kwargs,
1850 )
1851 plot_index.append(plot_filename)
1853 # Add list of plots to plot metadata.
1854 complete_plot_index = _append_to_plot_index(plot_index)
1856 # Make a page to display the plots.
1857 _make_plot_html_page(complete_plot_index)
1860def _get_num_models(cube: iris.cube.Cube | iris.cube.CubeList) -> int:
1861 """Return number of models based on cube attributes."""
1862 model_names = list(
1863 filter(
1864 lambda x: x is not None,
1865 {cb.attributes.get("model_name", None) for cb in iter_maybe(cube)},
1866 )
1867 )
1868 if not model_names:
1869 logging.debug("Missing model names. Will assume single model.")
1870 return 1
1871 else:
1872 return len(model_names)
1875def _validate_cube_shape(
1876 cube: iris.cube.Cube | iris.cube.CubeList, num_models: int
1877) -> None:
1878 """Check all cubes have a model name."""
1879 if isinstance(cube, iris.cube.CubeList) and len(cube) != num_models: 1879 ↛ 1880line 1879 didn't jump to line 1880 because the condition on line 1879 was never true
1880 raise ValueError(
1881 f"The number of model names ({num_models}) should equal the number "
1882 f"of cubes ({len(cube)})."
1883 )
1886def _validate_cubes_coords(
1887 cubes: iris.cube.CubeList, coords: list[iris.coords.Coord]
1888) -> None:
1889 """Check same number of cubes as sequence coordinate for zip functions."""
1890 if len(cubes) != len(coords): 1890 ↛ 1891line 1890 didn't jump to line 1891 because the condition on line 1890 was never true
1891 raise ValueError(
1892 f"The number of CubeList entries ({len(cubes)}) should equal the number "
1893 f"of sequence coordinates ({len(coords)})."
1894 f"Check that number of time entries in input data are consistent if "
1895 f"performing time-averaging steps prior to plotting outputs."
1896 )
1899####################
1900# Public functions #
1901####################
1904def spatial_contour_plot(
1905 cube: iris.cube.Cube,
1906 filename: str = None,
1907 sequence_coordinate: str = "time",
1908 stamp_coordinate: str = "realization",
1909 **kwargs,
1910) -> iris.cube.Cube:
1911 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1913 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1914 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1915 is present then postage stamp plots will be produced.
1917 Parameters
1918 ----------
1919 cube: Cube
1920 Iris cube of the data to plot. It should have two spatial dimensions,
1921 such as lat and lon, and may also have a another two dimension to be
1922 plotted sequentially and/or as postage stamp plots.
1923 filename: str, optional
1924 Name of the plot to write, used as a prefix for plot sequences. Defaults
1925 to the recipe name.
1926 sequence_coordinate: str, optional
1927 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1928 This coordinate must exist in the cube.
1929 stamp_coordinate: str, optional
1930 Coordinate about which to plot postage stamp plots. Defaults to
1931 ``"realization"``.
1933 Returns
1934 -------
1935 Cube
1936 The original cube (so further operations can be applied).
1938 Raises
1939 ------
1940 ValueError
1941 If the cube doesn't have the right dimensions.
1942 TypeError
1943 If the cube isn't a single cube.
1944 """
1945 _spatial_plot(
1946 "contourf", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1947 )
1948 return cube
1951def spatial_pcolormesh_plot(
1952 cube: iris.cube.Cube,
1953 filename: str = None,
1954 sequence_coordinate: str = "time",
1955 stamp_coordinate: str = "realization",
1956 **kwargs,
1957) -> iris.cube.Cube:
1958 """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube.
1960 A 2D spatial field can be plotted, but if the sequence_coordinate is present
1961 then a sequence of plots will be produced. Similarly if the stamp_coordinate
1962 is present then postage stamp plots will be produced.
1964 This function is significantly faster than ``spatial_contour_plot``,
1965 especially at high resolutions, and should be preferred unless contiguous
1966 contour areas are important.
1968 Parameters
1969 ----------
1970 cube: Cube
1971 Iris cube of the data to plot. It should have two spatial dimensions,
1972 such as lat and lon, and may also have a another two dimension to be
1973 plotted sequentially and/or as postage stamp plots.
1974 filename: str, optional
1975 Name of the plot to write, used as a prefix for plot sequences. Defaults
1976 to the recipe name.
1977 sequence_coordinate: str, optional
1978 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
1979 This coordinate must exist in the cube.
1980 stamp_coordinate: str, optional
1981 Coordinate about which to plot postage stamp plots. Defaults to
1982 ``"realization"``.
1984 Returns
1985 -------
1986 Cube
1987 The original cube (so further operations can be applied).
1989 Raises
1990 ------
1991 ValueError
1992 If the cube doesn't have the right dimensions.
1993 TypeError
1994 If the cube isn't a single cube.
1995 """
1996 _spatial_plot(
1997 "pcolormesh", cube, filename, sequence_coordinate, stamp_coordinate, **kwargs
1998 )
1999 return cube
2002def spatial_multi_pcolormesh_plot(
2003 cube: iris.cube.Cube,
2004 overlay_cube: iris.cube.Cube | None = None,
2005 contour_cube: iris.cube.Cube | None = None,
2006 scatter_cube: iris.cube.Cube | None = None,
2007 filename: str = None,
2008 sequence_coordinate: str = "time",
2009 stamp_coordinate: str = "realization",
2010 **kwargs,
2011) -> iris.cube.Cube:
2012 """Plot a set of spatial variables onto a map from a 2D, 3D, or 4D cube.
2014 A 2D basis cube spatial field can be plotted, but if the sequence_coordinate is present
2015 then a sequence of plots will be produced. Similarly if the stamp_coordinate
2016 is present then postage stamp plots will be produced.
2018 If specified, a masked overlay_cube can be overplotted on top of the base cube.
2020 If specified, contours of a contour_cube can be overplotted on top of those.
2022 If specified, a scattermap of scatter_cube can be overplotted.
2024 For single-variable equivalent of this routine, use spatial_pcolormesh_plot.
2026 This function is significantly faster than ``spatial_contour_plot``,
2027 especially at high resolutions, and should be preferred unless contiguous
2028 contour areas are important.
2030 Parameters
2031 ----------
2032 cube: Cube
2033 Iris cube of the data to plot. It should have two spatial dimensions,
2034 such as lat and lon, and may also have a another two dimension to be
2035 plotted sequentially and/or as postage stamp plots.
2036 overlay_cube: Cube
2037 Iris cube of the data to plot as an overlay on top of basis cube. It should have two spatial dimensions,
2038 such as lat and lon, and may also have a another two dimension to be
2039 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.
2040 contour_cube: Cube
2041 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,
2042 such as lat and lon, and may also have a another two dimension to be
2043 plotted sequentially and/or as postage stamp plots.
2044 scatter_cube: Cube
2045 Iris cube of the data to plot as a scatter map overlay on top of basis cube, overlay_cube and contour_cube. It should have two spatial dimensions,
2046 such as lat and lon, but these can describe a 1-D cube (e.g. list of
2047 observation stations with lat/lon coordinates) and may also have a another
2048 two dimension to be plotted sequentially and/or as postage stamp plots.
2049 filename: str, optional
2050 Name of the plot to write, used as a prefix for plot sequences. Defaults
2051 to the recipe name.
2052 sequence_coordinate: str, optional
2053 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2054 This coordinate must exist in the cube.
2055 stamp_coordinate: str, optional
2056 Coordinate about which to plot postage stamp plots. Defaults to
2057 ``"realization"``.
2059 Returns
2060 -------
2061 Cube
2062 The original cube (so further operations can be applied).
2064 Raises
2065 ------
2066 ValueError
2067 If the cube doesn't have the right dimensions.
2068 TypeError
2069 If the cube isn't a single cube.
2070 """
2071 _spatial_plot(
2072 "pcolormesh",
2073 cube,
2074 filename,
2075 sequence_coordinate,
2076 stamp_coordinate,
2077 overlay_cube=overlay_cube,
2078 contour_cube=contour_cube,
2079 scatter_cube=scatter_cube,
2080 )
2081 return cube, overlay_cube, contour_cube, scatter_cube
2084# TODO: Expand function to handle ensemble data.
2085# line_coordinate: str, optional
2086# Coordinate about which to plot multiple lines. Defaults to
2087# ``"realization"``.
2088def plot_line_series(
2089 cube: iris.cube.Cube | iris.cube.CubeList,
2090 filename: str = None,
2091 series_coordinate: str = "time",
2092 # line_coordinate: str = "realization",
2093 **kwargs,
2094) -> iris.cube.Cube | iris.cube.CubeList:
2095 """Plot a line plot for the specified coordinate.
2097 The Cube or CubeList must be 1D.
2099 Parameters
2100 ----------
2101 iris.cube | iris.cube.CubeList
2102 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2103 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2104 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2105 filename: str, optional
2106 Name of the plot to write, used as a prefix for plot sequences. Defaults
2107 to the recipe name.
2108 series_coordinate: str, optional
2109 Coordinate about which to make a series. Defaults to ``"time"``. This
2110 coordinate must exist in the cube.
2112 Returns
2113 -------
2114 iris.cube.Cube | iris.cube.CubeList
2115 The original Cube or CubeList (so further operations can be applied).
2116 plotted data.
2118 Raises
2119 ------
2120 ValueError
2121 If the cubes don't have the right dimensions.
2122 TypeError
2123 If the cube isn't a Cube or CubeList.
2124 """
2125 # Ensure we have a name for the plot file.
2126 recipe_title = get_recipe_metadata().get("title", "Untitled")
2128 num_models = _get_num_models(cube)
2130 _validate_cube_shape(cube, num_models)
2132 # Iterate over all cubes and extract coordinate to plot.
2133 cubes = iter_maybe(cube)
2134 coords = []
2135 for cube in cubes:
2136 try:
2137 coords.append(cube.coord(series_coordinate))
2138 except iris.exceptions.CoordinateNotFoundError as err:
2139 raise ValueError(
2140 f"Cube must have a {series_coordinate} coordinate."
2141 ) from err
2142 if cube.ndim > 2 or not cube.coords("realization"):
2143 raise ValueError("Cube must be 1D or 2D with a realization coordinate.")
2145 # Format the title and filename using plotted series coordinate
2146 nplot = 1
2147 seq_coord = coords[0]
2148 plot_title, plot_filename = _set_title_and_filename(
2149 seq_coord, nplot, recipe_title, filename
2150 )
2152 # Do the actual plotting.
2154 # Treat cubes with station coordinate as point observation timeseries, looping over available points
2155 if "station" in [c.name() for c in cubes[0].coords()]: 2155 ↛ 2156line 2155 didn't jump to line 2156 because the condition on line 2155 was never true
2156 for station in cubes[0].coord("station").points:
2157 station_cubes = cubes.extract(iris.Constraint(station=station))
2158 station_name = station_cubes[0].coord("Station identifier").points[0]
2159 _plot_and_save_line_series(
2160 station_cubes,
2161 coords,
2162 "realization",
2163 plot_filename.replace(".png", "_" + station_name + ".png"),
2164 f"{plot_title} {station_name}",
2165 )
2167 # Default line plotting
2168 else:
2169 _plot_and_save_line_series(
2170 cubes, coords, "realization", plot_filename, plot_title
2171 )
2173 # Add list of plots to plot metadata.
2174 plot_index = _append_to_plot_index([plot_filename])
2176 # Make a page to display the plots.
2177 _make_plot_html_page(plot_index)
2179 return cube
2182def plot_vertical_line_series(
2183 cubes: iris.cube.Cube | iris.cube.CubeList,
2184 filename: str = None,
2185 series_coordinate: str = "model_level_number",
2186 sequence_coordinate: str = "time",
2187 # line_coordinate: str = "realization",
2188 **kwargs,
2189) -> iris.cube.Cube | iris.cube.CubeList:
2190 """Plot a line plot against a type of vertical coordinate.
2192 The Cube or CubeList must be 1D.
2194 A 1D line plot with y-axis as pressure coordinate can be plotted, but if the sequence_coordinate is present
2195 then a sequence of plots will be produced.
2197 Parameters
2198 ----------
2199 iris.cube | iris.cube.CubeList
2200 Cube or CubeList of the data to plot. The individual cubes should have a single dimension.
2201 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2202 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2203 filename: str, optional
2204 Name of the plot to write, used as a prefix for plot sequences. Defaults
2205 to the recipe name.
2206 series_coordinate: str, optional
2207 Coordinate to plot on the y-axis. Can be ``pressure`` or
2208 ``model_level_number`` for UM, or ``full_levels`` or ``half_levels``
2209 for LFRic. Defaults to ``model_level_number``.
2210 This coordinate must exist in the cube.
2211 sequence_coordinate: str, optional
2212 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2213 This coordinate must exist in the cube.
2215 Returns
2216 -------
2217 iris.cube.Cube | iris.cube.CubeList
2218 The original Cube or CubeList (so further operations can be applied).
2219 Plotted data.
2221 Raises
2222 ------
2223 ValueError
2224 If the cubes doesn't have the right dimensions.
2225 TypeError
2226 If the cube isn't a Cube or CubeList.
2227 """
2228 # Ensure we have a name for the plot file.
2229 recipe_title = get_recipe_metadata().get("title", "Untitled")
2231 cubes = iter_maybe(cubes)
2232 # Initialise empty list to hold all data from all cubes in a CubeList
2233 all_data = []
2235 # Store min/max ranges for x range.
2236 x_levels = []
2238 num_models = _get_num_models(cubes)
2240 _validate_cube_shape(cubes, num_models)
2242 # Iterate over all cubes in cube or CubeList and plot.
2243 coords = []
2244 for cube in cubes:
2245 # Test if series coordinate i.e. pressure level exist for any cube with cube.ndim >=1.
2246 try:
2247 coords.append(cube.coord(series_coordinate))
2248 except iris.exceptions.CoordinateNotFoundError as err:
2249 raise ValueError(
2250 f"Cube must have a {series_coordinate} coordinate."
2251 ) from err
2253 try:
2254 if cube.ndim > 1 or not cube.coords("realization"): 2254 ↛ 2262line 2254 didn't jump to line 2262 because the condition on line 2254 was always true
2255 cube.coord(sequence_coordinate)
2256 except iris.exceptions.CoordinateNotFoundError as err:
2257 raise ValueError(
2258 f"Cube must have a {sequence_coordinate} coordinate or be 1D, or 2D with a realization coordinate."
2259 ) from err
2261 # Get minimum and maximum from levels information.
2262 _, levels, _ = _colorbar_map_levels(cube, axis="x")
2263 if levels is not None: 2263 ↛ 2267line 2263 didn't jump to line 2267 because the condition on line 2263 was always true
2264 x_levels.append(min(levels))
2265 x_levels.append(max(levels))
2266 else:
2267 all_data.append(cube.data)
2269 if len(x_levels) == 0: 2269 ↛ 2271line 2269 didn't jump to line 2271 because the condition on line 2269 was never true
2270 # Combine all data into a single NumPy array
2271 combined_data = np.concatenate(all_data)
2273 # Set the lower and upper limit for the x-axis to ensure all plots have
2274 # same range. This needs to read the whole cube over the range of the
2275 # sequence and if applicable postage stamp coordinate.
2276 vmin = np.floor(combined_data.min())
2277 vmax = np.ceil(combined_data.max())
2278 else:
2279 vmin = min(x_levels)
2280 vmax = max(x_levels)
2282 # Matching the slices (matching by seq coord point; it may happen that
2283 # evaluated models do not cover the same seq coord range, hence matching
2284 # necessary)
2285 def filter_cube_iterables(cube_iterables) -> bool:
2286 return len(cube_iterables) == len(coords)
2288 cube_iterables = filter(
2289 filter_cube_iterables,
2290 (
2291 iris.cube.CubeList(
2292 s
2293 for s in itertools.chain.from_iterable(
2294 cb.slices_over(sequence_coordinate) for cb in cubes
2295 )
2296 if s.coord(sequence_coordinate).points[0] == point
2297 )
2298 for point in sorted(
2299 set(
2300 itertools.chain.from_iterable(
2301 cb.coord(sequence_coordinate).points for cb in cubes
2302 )
2303 )
2304 )
2305 ),
2306 )
2308 # Create a plot for each value of the sequence coordinate.
2309 # Allowing for multiple cubes in a CubeList to be plotted in the same plot for
2310 # similar sequence values. Passing a CubeList into the internal plotting function
2311 # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube
2312 # or an iris.cube.CubeList.
2313 plot_index = []
2314 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2315 for cubes_slice in cube_iterables:
2316 # Format the coordinate value in a unit appropriate way.
2317 seq_coord = cubes_slice[0].coord(sequence_coordinate)
2318 plot_title, plot_filename = _set_title_and_filename(
2319 seq_coord, nplot, recipe_title, filename
2320 )
2322 # Do the actual plotting.
2323 _plot_and_save_vertical_line_series(
2324 cubes_slice,
2325 coords,
2326 "realization",
2327 plot_filename,
2328 series_coordinate,
2329 title=plot_title,
2330 vmin=vmin,
2331 vmax=vmax,
2332 )
2333 plot_index.append(plot_filename)
2335 # Add list of plots to plot metadata.
2336 complete_plot_index = _append_to_plot_index(plot_index)
2338 # Make a page to display the plots.
2339 _make_plot_html_page(complete_plot_index)
2341 return cubes
2344def qq_plot(
2345 cubes: iris.cube.CubeList,
2346 coordinates: list[str],
2347 percentiles: list[float],
2348 model_names: list[str],
2349 filename: str = None,
2350 one_to_one: bool = True,
2351 **kwargs,
2352) -> iris.cube.CubeList:
2353 """Plot a Quantile-Quantile plot between two models for common time points.
2355 The cubes will be normalised by collapsing each cube to its percentiles. Cubes are
2356 collapsed within the operator over all specified coordinates such as
2357 grid_latitude, grid_longitude, vertical levels, but also realisation representing
2358 ensemble members to ensure a 1D cube (array).
2360 Parameters
2361 ----------
2362 cubes: iris.cube.CubeList
2363 Two cubes of the same variable with different models.
2364 coordinate: list[str]
2365 The list of coordinates to collapse over. This list should be
2366 every coordinate within the cube to result in a 1D cube around
2367 the percentile coordinate.
2368 percent: list[float]
2369 A list of percentiles to appear in the plot.
2370 model_names: list[str]
2371 A list of model names to appear on the axis of the plot.
2372 filename: str, optional
2373 Filename of the plot to write.
2374 one_to_one: bool, optional
2375 If True a 1:1 line is plotted; if False it is not. Default is True.
2377 Raises
2378 ------
2379 ValueError
2380 When the cubes are not compatible.
2382 Notes
2383 -----
2384 The quantile-quantile plot is a variant on the scatter plot representing
2385 two datasets by their quantiles (percentiles) for common time points.
2386 This plot does not use a theoretical distribution to compare against, but
2387 compares percentiles of two datasets. This plot does
2388 not use all raw data points, but plots the selected percentiles (quantiles) of
2389 each variable instead for the two datasets, thereby normalising the data for a
2390 direct comparison between the selected percentiles of the two dataset distributions.
2392 Quantile-quantile plots are valuable for comparing against
2393 observations and other models. Identical percentiles between the variables
2394 will lie on the one-to-one line implying the values correspond well to each
2395 other. Where there is a deviation from the one-to-one line a range of
2396 possibilities exist depending on how and where the data is shifted (e.g.,
2397 Wilks 2011 [Wilks2011]_).
2399 For distributions above the one-to-one line the distribution is left-skewed;
2400 below is right-skewed. A distinct break implies a bimodal distribution, and
2401 closer values/values further apart at the tails imply poor representation of
2402 the extremes.
2404 References
2405 ----------
2406 .. [Wilks2011] Wilks, D.S., (2011) "Statistical Methods in the Atmospheric
2407 Sciences" Third Edition, vol. 100, Academic Press, Oxford, UK, 676 pp.
2408 """
2409 # Check cubes using same functionality as the difference operator.
2410 if len(cubes) != 2:
2411 raise ValueError("cubes should contain exactly 2 cubes.")
2412 base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1))
2413 other: Cube = cubes.extract_cube(
2414 iris.Constraint(
2415 cube_func=lambda cube: "cset_comparison_base" not in cube.attributes
2416 )
2417 )
2419 # Get spatial coord names.
2420 base_lat_name, base_lon_name = get_cube_yxcoordname(base)
2421 other_lat_name, other_lon_name = get_cube_yxcoordname(other)
2423 # Ensure cubes to compare are on common differencing grid.
2424 # This is triggered if either
2425 # i) latitude and longitude shapes are not the same. Note grid points
2426 # are not compared directly as these can differ through rounding
2427 # errors.
2428 # ii) or variables are known to often sit on different grid staggering
2429 # in different models (e.g. cell center vs cell edge), as is the case
2430 # for UM and LFRic comparisons.
2431 # In future greater choice of regridding method might be applied depending
2432 # on variable type. Linear regridding can in general be appropriate for smooth
2433 # variables. Care should be taken with interpretation of differences
2434 # given this dependency on regridding.
2435 if (
2436 base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
2437 or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
2438 ) or (
2439 base.long_name
2440 in [
2441 "eastward_wind_at_10m",
2442 "northward_wind_at_10m",
2443 "northward_wind_at_cell_centres",
2444 "eastward_wind_at_cell_centres",
2445 "zonal_wind_at_pressure_levels",
2446 "meridional_wind_at_pressure_levels",
2447 "potential_vorticity_at_pressure_levels",
2448 "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
2449 ]
2450 ):
2451 logging.debug(
2452 "Linear regridding base cube to other grid to compute differences"
2453 )
2454 base = regrid_onto_cube(base, other, method="Linear")
2456 # Extract just common time points.
2457 base, other = _extract_common_time_points(base, other)
2459 # Equalise attributes so we can merge.
2460 fully_equalise_attributes([base, other])
2461 logging.debug("Base: %s\nOther: %s", base, other)
2463 # Collapse cubes.
2464 base = collapse(
2465 base,
2466 coordinate=coordinates,
2467 method="PERCENTILE",
2468 additional_percent=percentiles,
2469 )
2470 other = collapse(
2471 other,
2472 coordinate=coordinates,
2473 method="PERCENTILE",
2474 additional_percent=percentiles,
2475 )
2477 # Ensure we have a name for the plot file.
2478 recipe_title = get_recipe_metadata().get("title", "Untitled")
2479 title = f"{recipe_title}"
2481 if filename is None:
2482 filename = slugify(recipe_title)
2484 # Add file extension.
2485 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2487 # Do the actual plotting on a scatter plot
2488 _plot_and_save_scatter_plot(
2489 base, other, plot_filename, title, one_to_one, model_names
2490 )
2492 # Add list of plots to plot metadata.
2493 plot_index = _append_to_plot_index([plot_filename])
2495 # Make a page to display the plots.
2496 _make_plot_html_page(plot_index)
2498 return iris.cube.CubeList([base, other])
2501def scatter_plot(
2502 cube_x: iris.cube.Cube | iris.cube.CubeList,
2503 cube_y: iris.cube.Cube | iris.cube.CubeList,
2504 filename: str = None,
2505 one_to_one: bool = True,
2506 **kwargs,
2507) -> iris.cube.CubeList:
2508 """Plot a scatter plot between two variables.
2510 Both cubes must be 1D.
2512 Parameters
2513 ----------
2514 cube_x: Cube | CubeList
2515 1 dimensional Cube of the data to plot on y-axis.
2516 cube_y: Cube | CubeList
2517 1 dimensional Cube of the data to plot on x-axis.
2518 filename: str, optional
2519 Filename of the plot to write.
2520 one_to_one: bool, optional
2521 If True a 1:1 line is plotted; if False it is not. Default is True.
2523 Returns
2524 -------
2525 cubes: CubeList
2526 CubeList of the original x and y cubes for further processing.
2528 Raises
2529 ------
2530 ValueError
2531 If the cube doesn't have the right dimensions and cubes not the same
2532 size.
2533 TypeError
2534 If the cube isn't a single cube.
2536 Notes
2537 -----
2538 Scatter plots are used for determining if there is a relationship between
2539 two variables. Positive relations have a slope going from bottom left to top
2540 right; Negative relations have a slope going from top left to bottom right.
2541 """
2542 # Iterate over all cubes in cube or CubeList and plot.
2543 for cube_iter in iter_maybe(cube_x):
2544 # Check cubes are correct shape.
2545 cube_iter = _check_single_cube(cube_iter)
2546 if cube_iter.ndim > 1:
2547 raise ValueError("cube_x must be 1D.")
2549 # Iterate over all cubes in cube or CubeList and plot.
2550 for cube_iter in iter_maybe(cube_y):
2551 # Check cubes are correct shape.
2552 cube_iter = _check_single_cube(cube_iter)
2553 if cube_iter.ndim > 1:
2554 raise ValueError("cube_y must be 1D.")
2556 # Ensure we have a name for the plot file.
2557 recipe_title = get_recipe_metadata().get("title", "Untitled")
2558 title = f"{recipe_title}"
2560 if filename is None:
2561 filename = slugify(recipe_title)
2563 # Add file extension.
2564 plot_filename = f"{filename.rsplit('.', 1)[0]}.png"
2566 # Do the actual plotting.
2567 _plot_and_save_scatter_plot(cube_x, cube_y, plot_filename, title, one_to_one)
2569 # Add list of plots to plot metadata.
2570 plot_index = _append_to_plot_index([plot_filename])
2572 # Make a page to display the plots.
2573 _make_plot_html_page(plot_index)
2575 return iris.cube.CubeList([cube_x, cube_y])
2578def vector_plot(
2579 cube_u: iris.cube.Cube,
2580 cube_v: iris.cube.Cube,
2581 filename: str = None,
2582 sequence_coordinate: str = "time",
2583 **kwargs,
2584) -> iris.cube.CubeList:
2585 """Plot a vector plot based on the input u and v components."""
2586 recipe_title = get_recipe_metadata().get("title", "Untitled")
2588 # Cubes must have a matching sequence coordinate.
2589 try:
2590 # Check that the u and v cubes have the same sequence coordinate.
2591 if cube_u.coord(sequence_coordinate) != cube_v.coord(sequence_coordinate): 2591 ↛ anywhereline 2591 didn't jump anywhere: it always raised an exception.
2592 raise ValueError("Coordinates do not match.")
2593 except (iris.exceptions.CoordinateNotFoundError, ValueError) as err:
2594 raise ValueError(
2595 f"Cubes should have matching {sequence_coordinate} coordinate:\n{cube_u}\n{cube_v}"
2596 ) from err
2598 # Create a plot for each value of the sequence coordinate.
2599 plot_index = []
2600 nplot = np.size(cube_u[0].coord(sequence_coordinate).points)
2601 for cube_u_slice, cube_v_slice in zip(
2602 cube_u.slices_over(sequence_coordinate),
2603 cube_v.slices_over(sequence_coordinate),
2604 strict=True,
2605 ):
2606 # Format the coordinate value in a unit appropriate way.
2607 seq_coord = cube_u_slice.coord(sequence_coordinate)
2608 plot_title, plot_filename = _set_title_and_filename(
2609 seq_coord, nplot, recipe_title, filename
2610 )
2612 # Do the actual plotting.
2613 _plot_and_save_vector_plot(
2614 cube_u_slice,
2615 cube_v_slice,
2616 filename=plot_filename,
2617 title=plot_title,
2618 method="contourf",
2619 )
2620 plot_index.append(plot_filename)
2622 # Add list of plots to plot metadata.
2623 complete_plot_index = _append_to_plot_index(plot_index)
2625 # Make a page to display the plots.
2626 _make_plot_html_page(complete_plot_index)
2628 return iris.cube.CubeList([cube_u, cube_v])
2631def _check_sequence(cubes, sequence_coordinate):
2632 # If several histograms are plotted with time as sequence_coordinate for the
2633 # time slider option.
2634 for cube in cubes:
2635 try:
2636 cube.coord(sequence_coordinate)
2637 except iris.exceptions.CoordinateNotFoundError as err:
2638 raise ValueError(
2639 f"Cube must have a {sequence_coordinate} coordinate."
2640 ) from err
2643def _set_axis_range(cubes):
2644 # Get minimum and maximum from levels information.
2645 levels = None
2646 for cube in cubes: 2646 ↛ 2662line 2646 didn't jump to line 2662 because the loop on line 2646 didn't complete
2647 # First check if user-specified "auto" range variable.
2648 # This maintains the value of levels as None, so proceed.
2649 _, levels, _ = _colorbar_map_levels(cube, axis="y")
2650 if levels is None:
2651 break
2652 # If levels is changed, recheck to use the vmin,vmax or
2653 # levels-based ranges for histogram plots.
2654 _, levels, _ = _colorbar_map_levels(cube)
2655 logging.debug("levels: %s", levels)
2656 if levels is not None: 2656 ↛ 2646line 2656 didn't jump to line 2646 because the condition on line 2656 was always true
2657 vmin = min(levels)
2658 vmax = max(levels)
2659 logging.debug("Updated vmin, vmax: %s, %s", vmin, vmax)
2660 break
2662 if levels is None:
2663 vmin = min(cb.data.min() for cb in cubes)
2664 vmax = max(cb.data.max() for cb in cubes)
2666 return vmin, vmax
2669def _find_matched_slices(cubes, sequence_coordinate):
2671 all_points = sorted(
2672 set(
2673 itertools.chain.from_iterable(
2674 cb.coord(sequence_coordinate).points for cb in cubes
2675 )
2676 )
2677 )
2678 all_slices = list(
2679 itertools.chain.from_iterable(
2680 cb.slices_over(sequence_coordinate) for cb in cubes
2681 )
2682 )
2683 # Matched slices (matched by seq coord point; it may happen that
2684 # evaluated models do not cover the same seq coord range, hence matching
2685 # necessary)
2686 cube_iterables = [
2687 iris.cube.CubeList(
2688 s for s in all_slices if s.coord(sequence_coordinate).points[0] == point
2689 )
2690 for point in all_points
2691 ]
2693 return cube_iterables
2696def plot_histogram_series(
2697 cubes: iris.cube.Cube | iris.cube.CubeList,
2698 filename: str = None,
2699 sequence_coordinate: str = "time",
2700 stamp_coordinate: str = "realization",
2701 single_plot: bool = False,
2702 **kwargs,
2703) -> iris.cube.Cube | iris.cube.CubeList:
2704 """Plot a histogram plot for each vertical level provided.
2706 A histogram plot can be plotted, but if the sequence_coordinate (i.e. time)
2707 is present then a sequence of plots will be produced using the time slider
2708 functionality to scroll through histograms against time. If a
2709 stamp_coordinate is present then postage stamp plots will be produced. If
2710 stamp_coordinate and single_plot is True, all postage stamp plots will be
2711 plotted in a single plot instead of separate postage stamp plots.
2713 Parameters
2714 ----------
2715 cubes: Cube | iris.cube.CubeList
2716 Iris cube or CubeList of the data to plot. It should have a single dimension other
2717 than the stamp coordinate.
2718 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2719 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2720 filename: str, optional
2721 Name of the plot to write, used as a prefix for plot sequences. Defaults
2722 to the recipe name.
2723 sequence_coordinate: str, optional
2724 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2725 This coordinate must exist in the cube and will be used for the time
2726 slider.
2727 stamp_coordinate: str, optional
2728 Coordinate about which to plot postage stamp plots. Defaults to
2729 ``"realization"``.
2730 single_plot: bool, optional
2731 If True, all postage stamp plots will be plotted in a single plot. If
2732 False, each postage stamp plot will be plotted separately. Is only valid
2733 if stamp_coordinate exists and has more than a single point.
2735 Returns
2736 -------
2737 iris.cube.Cube | iris.cube.CubeList
2738 The original Cube or CubeList (so further operations can be applied).
2739 Plotted data.
2741 Raises
2742 ------
2743 ValueError
2744 If the cube doesn't have the right dimensions.
2745 TypeError
2746 If the cube isn't a Cube or CubeList.
2747 """
2748 recipe_title = get_recipe_metadata().get("title", "Untitled")
2750 cubes = iter_maybe(cubes)
2752 # Internal plotting function.
2753 plotting_func = _plot_and_save_histogram_series
2755 num_models = _get_num_models(cubes)
2757 _validate_cube_shape(cubes, num_models)
2759 _check_sequence(cubes, sequence_coordinate)
2761 vmin, vmax = _set_axis_range(cubes)
2763 # Make postage stamp plots if stamp_coordinate exists and has more than a
2764 # single point. If single_plot is True:
2765 # -- all postage stamp plots will be plotted in a single plot instead of
2766 # separate postage stamp plots.
2767 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2768 # produced per single model only
2769 if num_models == 1: 2769 ↛ 2782line 2769 didn't jump to line 2782 because the condition on line 2769 was always true
2770 if ( 2770 ↛ 2774line 2770 didn't jump to line 2774 because the condition on line 2770 was never true
2771 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2772 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2773 ):
2774 if single_plot:
2775 plotting_func = (
2776 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2777 )
2778 else:
2779 plotting_func = _plot_and_save_postage_stamp_histogram_series
2780 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2781 else:
2782 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2784 plot_index = []
2785 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2786 # Create a plot for each value of the sequence coordinate. Allowing for
2787 # multiple cubes in a CubeList to be plotted in the same plot for similar
2788 # sequence values. Passing a CubeList into the internal plotting function
2789 # for similar values of the sequence coordinate. cube_slice can be an
2790 # iris.cube.Cube or an iris.cube.CubeList.
2791 for cube_slice in cube_iterables:
2792 single_cube = cube_slice
2793 if isinstance(cube_slice, iris.cube.CubeList): 2793 ↛ 2794line 2793 didn't jump to line 2794 because the condition on line 2793 was never true
2794 single_cube = cube_slice[0]
2796 # Ensure valid stamp coordinate in cube dimensions
2797 if stamp_coordinate == "realization": 2797 ↛ 2800line 2797 didn't jump to line 2800 because the condition on line 2797 was always true
2798 stamp_coordinate = check_stamp_coordinate(single_cube)
2799 # Set plot titles and filename, based on sequence coordinate
2800 seq_coord = single_cube.coord(sequence_coordinate)
2801 # Use time coordinate in title and filename if single histogram output.
2802 if sequence_coordinate == "realization" and nplot == 1: 2802 ↛ 2803line 2802 didn't jump to line 2803 because the condition on line 2802 was never true
2803 seq_coord = single_cube.coord("time")
2804 # Use station name in title and filename if model vs obs comparison
2805 if sequence_coordinate == "station": 2805 ↛ 2806line 2805 didn't jump to line 2806 because the condition on line 2805 was never true
2806 seq_coord = single_cube.coord("Station identifier")
2808 plot_title, plot_filename = _set_title_and_filename(
2809 seq_coord, nplot, recipe_title, filename
2810 )
2812 # Do the actual plotting.
2813 plotting_func(
2814 cube_slice,
2815 filename=plot_filename,
2816 stamp_coordinate=stamp_coordinate,
2817 title=plot_title,
2818 vmin=vmin,
2819 vmax=vmax,
2820 )
2821 plot_index.append(plot_filename)
2823 # Add list of plots to plot metadata.
2824 complete_plot_index = _append_to_plot_index(plot_index)
2826 # Make a page to display the plots.
2827 _make_plot_html_page(complete_plot_index)
2829 return cubes
2832def plot_scatter_series(
2833 cubes: iris.cube.Cube | iris.cube.CubeList,
2834 filename: str = None,
2835 sequence_coordinate: str = "time",
2836 stamp_coordinate: str = "realization",
2837 single_plot: bool = False,
2838 scatter: bool = False,
2839 **kwargs,
2840) -> iris.cube.Cube | iris.cube.CubeList:
2841 """Plot a scatter plot for each sequence coordinate provided.
2843 A scatter plot can be plotted, but if the sequence_coordinate (i.e. time)
2844 is present then a sequence of plots will be produced using the time slider
2845 functionality to scroll through scatter against time. If a
2846 stamp_coordinate is present then postage stamp plots will be produced. If
2847 stamp_coordinate and single_plot is True, all postage stamp plots will be
2848 plotted in a single plot instead of separate postage stamp plots.
2850 Parameters
2851 ----------
2852 cubes: Cube | iris.cube.CubeList
2853 Iris cube or CubeList of the data to plot. It should have a single dimension other
2854 than the stamp coordinate.
2855 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2856 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2857 filename: str, optional
2858 Name of the plot to write, used as a prefix for plot sequences. Defaults
2859 to the recipe name.
2860 sequence_coordinate: str, optional
2861 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2862 This coordinate must exist in the cube and will be used for the time
2863 slider.
2864 stamp_coordinate: str, optional
2865 Coordinate about which to plot postage stamp plots. Defaults to
2866 ``"realization"``.
2867 single_plot: bool, optional
2868 If True, all postage stamp plots will be plotted in a single plot. If
2869 False, each postage stamp plot will be plotted separately. Is only valid
2870 if stamp_coordinate exists and has more than a single point.
2872 Returns
2873 -------
2874 iris.cube.Cube | iris.cube.CubeList
2875 The original Cube or CubeList (so further operations can be applied).
2876 Plotted data.
2878 Raises
2879 ------
2880 ValueError
2881 If the cube doesn't have the right dimensions.
2882 TypeError
2883 If the cube isn't a Cube or CubeList.
2884 """
2885 recipe_title = get_recipe_metadata().get("title", "Untitled")
2887 cubes = iter_maybe(cubes)
2889 # Internal plotting function.
2890 plotting_func = _plot_and_save_scatter_series
2892 num_models = _get_num_models(cubes)
2894 _validate_cube_shape(cubes, num_models)
2896 _check_sequence(cubes, sequence_coordinate)
2898 vmin, vmax = _set_axis_range(cubes)
2900 # Make postage stamp plots if stamp_coordinate exists and has more than a
2901 # single point. If single_plot is True:
2902 # -- all postage stamp plots will be plotted in a single plot instead of
2903 # separate postage stamp plots.
2904 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
2905 # produced per single model only
2906 if num_models == 1:
2907 if (
2908 stamp_coordinate in [c.name() for c in cubes[0].coords()]
2909 and cubes[0].coord(stamp_coordinate).shape[0] > 1
2910 ):
2911 if single_plot:
2912 plotting_func = (
2913 _plot_and_save_postage_stamps_in_single_plot_histogram_series
2914 )
2915 else:
2916 plotting_func = _plot_and_save_postage_stamp_histogram_series
2917 cube_iterables = cubes[0].slices_over(sequence_coordinate)
2918 else:
2919 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
2921 plot_index = []
2922 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
2923 # Create a plot for each value of the sequence coordinate. Allowing for
2924 # multiple cubes in a CubeList to be plotted in the same plot for similar
2925 # sequence values. Passing a CubeList into the internal plotting function
2926 # for similar values of the sequence coordinate. cube_slice can be an
2927 # iris.cube.Cube or an iris.cube.CubeList.
2928 for cube_slice in cube_iterables:
2929 single_cube = cube_slice
2930 if isinstance(cube_slice, iris.cube.CubeList):
2931 single_cube = cube_slice[0]
2933 # Ensure valid stamp coordinate in cube dimensions
2934 if stamp_coordinate == "realization":
2935 stamp_coordinate = check_stamp_coordinate(single_cube)
2936 # Set plot titles and filename, based on sequence coordinate
2937 seq_coord = single_cube.coord(sequence_coordinate)
2938 # Use time coordinate in title and filename if single histogram output.
2939 if sequence_coordinate == "realization" and nplot == 1:
2940 seq_coord = single_cube.coord("time")
2941 # Use station name in title and filename if model vs obs comparison
2942 if sequence_coordinate == "station":
2943 seq_coord = single_cube.coord("Station identifier")
2945 plot_title, plot_filename = _set_title_and_filename(
2946 seq_coord, nplot, recipe_title, filename
2947 )
2949 # Do the actual plotting.
2950 plotting_func(
2951 cube_slice,
2952 filename=plot_filename,
2953 stamp_coordinate=stamp_coordinate,
2954 title=plot_title,
2955 vmin=vmin,
2956 vmax=vmax,
2957 )
2958 plot_index.append(plot_filename)
2960 # Add list of plots to plot metadata.
2961 complete_plot_index = _append_to_plot_index(plot_index)
2963 # Make a page to display the plots.
2964 _make_plot_html_page(complete_plot_index)
2966 return cubes
2969def plot_power_spectrum_series(
2970 cubes: iris.cube.Cube | iris.cube.CubeList,
2971 filename: str = None,
2972 sequence_coordinate: str = "time",
2973 stamp_coordinate: str = "realization",
2974 single_plot: bool = False,
2975 **kwargs,
2976) -> iris.cube.Cube | iris.cube.CubeList:
2977 """Plot a power spectrum plot for each vertical level provided.
2979 A power spectrum plot can be plotted, but if the sequence_coordinate (i.e. time)
2980 is present then a sequence of plots will be produced using the time slider
2981 functionality to scroll through power spectrum against time. If a
2982 stamp_coordinate is present then postage stamp plots will be produced. If
2983 stamp_coordinate and single_plot is True, all postage stamp plots will be
2984 plotted in a single plot instead of separate postage stamp plots.
2986 Parameters
2987 ----------
2988 cubes: Cube | iris.cube.CubeList
2989 Iris cube or CubeList of the data to plot. It should have a single dimension other
2990 than the stamp coordinate.
2991 The cubes should cover the same phenomenon i.e. all cubes contain temperature data.
2992 We do not support different data such as temperature and humidity in the same CubeList for plotting.
2993 filename: str, optional
2994 Name of the plot to write, used as a prefix for plot sequences. Defaults
2995 to the recipe name.
2996 sequence_coordinate: str, optional
2997 Coordinate about which to make a plot sequence. Defaults to ``"time"``.
2998 This coordinate must exist in the cube and will be used for the time
2999 slider.
3000 stamp_coordinate: str, optional
3001 Coordinate about which to plot postage stamp plots. Defaults to
3002 ``"realization"``.
3003 single_plot: bool, optional
3004 If True, all postage stamp plots will be plotted in a single plot. If
3005 False, each postage stamp plot will be plotted separately. Is only valid
3006 if stamp_coordinate exists and has more than a single point.
3008 Returns
3009 -------
3010 iris.cube.Cube | iris.cube.CubeList
3011 The original Cube or CubeList (so further operations can be applied).
3012 Plotted data.
3014 Raises
3015 ------
3016 ValueError
3017 If the cube doesn't have the right dimensions.
3018 TypeError
3019 If the cube isn't a Cube or CubeList.
3020 """
3021 recipe_title = get_recipe_metadata().get("title", "Untitled")
3023 cubes = iter_maybe(cubes)
3025 # Internal plotting function.
3026 plotting_func = _plot_and_save_power_spectrum_series
3028 num_models = _get_num_models(cubes)
3030 _validate_cube_shape(cubes, num_models)
3032 _check_sequence(cubes, sequence_coordinate)
3034 # Make postage stamp plots if stamp_coordinate exists and has more than a
3035 # single point. If single_plot is True:
3036 # -- all postage stamp plots will be plotted in a single plot instead of
3037 # separate postage stamp plots.
3038 # -- model names (hidden in cube attrs) are ignored, that is stamp plots are
3039 # produced per single model only
3040 if num_models == 1: 3040 ↛ 3053line 3040 didn't jump to line 3053 because the condition on line 3040 was always true
3041 if ( 3041 ↛ 3045line 3041 didn't jump to line 3045 because the condition on line 3041 was never true
3042 stamp_coordinate in [c.name() for c in cubes[0].coords()]
3043 and cubes[0].coord(stamp_coordinate).shape[0] > 1
3044 ):
3045 if single_plot:
3046 plotting_func = (
3047 _plot_and_save_postage_stamps_in_single_plot_power_spectrum_series
3048 )
3049 else:
3050 plotting_func = _plot_and_save_postage_stamp_power_spectrum_series
3051 cube_iterables = cubes[0].slices_over(sequence_coordinate)
3052 else:
3053 cube_iterables = _find_matched_slices(cubes, sequence_coordinate)
3055 plot_index = []
3056 nplot = np.size(cubes[0].coord(sequence_coordinate).points)
3057 # Create a plot for each value of the sequence coordinate. Allowing for
3058 # multiple cubes in a CubeList to be plotted in the same plot for similar
3059 # sequence values. Passing a CubeList into the internal plotting function
3060 # for similar values of the sequence coordinate. cube_slice can be an
3061 # iris.cube.Cube or an iris.cube.CubeList.
3062 for cube_slice in cube_iterables:
3063 single_cube = cube_slice
3064 if isinstance(cube_slice, iris.cube.CubeList): 3064 ↛ 3065line 3064 didn't jump to line 3065 because the condition on line 3064 was never true
3065 single_cube = cube_slice[0]
3067 # Set stamp coordinate
3068 if stamp_coordinate == "realization": 3068 ↛ 3071line 3068 didn't jump to line 3071 because the condition on line 3068 was always true
3069 stamp_coordinate = check_stamp_coordinate(single_cube)
3070 # Set plot title and filenames based on sequence values
3071 seq_coord = single_cube.coord(sequence_coordinate)
3072 plot_title, plot_filename = _set_title_and_filename(
3073 seq_coord, nplot, recipe_title, filename
3074 )
3076 # Do the actual plotting.
3077 plotting_func(
3078 cube_slice,
3079 filename=plot_filename,
3080 stamp_coordinate=stamp_coordinate,
3081 title=plot_title,
3082 )
3083 plot_index.append(plot_filename)
3085 # Add list of plots to plot metadata.
3086 complete_plot_index = _append_to_plot_index(plot_index)
3088 # Make a page to display the plots.
3089 _make_plot_html_page(complete_plot_index)
3091 return cubes
3094def _DCT_ps(y_3d):
3095 """Calculate power spectra for regional domains.
3097 Parameters
3098 ----------
3099 y_3d: 3D array
3100 3 dimensional array to calculate spectrum for.
3101 (2D field data with 3rd dimension of time)
3103 Returns: ps_array
3104 Array of power spectra values calculated for input field (for each time)
3106 Method for regional domains:
3107 Calculate power spectra over limited area domain using Discrete Cosine Transform (DCT)
3108 as described in Denis et al 2002 [Denis_etal_2002]_.
3110 References
3111 ----------
3112 .. [Denis_etal_2002] Bertrand Denis, Jean Côté and René Laprise (2002)
3113 "Spectral Decomposition of Two-Dimensional Atmospheric Fields on
3114 Limited-Area Domains Using the Discrete Cosine Transform (DCT)"
3115 Monthly Weather Review, Vol. 130, 1812-1828
3116 doi: https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2
3117 """
3118 Nt, Ny, Nx = y_3d.shape
3120 # Max coefficient
3121 Nmin = min(Nx - 1, Ny - 1)
3123 # Create alpha matrix (of wavenumbers)
3124 alpha_matrix = _create_alpha_matrix(Ny, Nx)
3126 # Prepare output array
3127 ps_array = np.zeros((Nt, Nmin))
3129 # Loop over time to get spectrum for each time.
3130 for t in range(Nt):
3131 y_2d = y_3d[t]
3133 # Apply 2D DCT to transform y_3d[t] from physical space to spectral space.
3134 # fkk is a 2D array of DCT coefficients, representing the amplitudes of
3135 # cosine basis functions at different spatial frequencies.
3137 # normalise spectrum to allow comparison between models.
3138 fkk = fft.dctn(y_2d, norm="ortho")
3140 # Normalise fkk
3141 fkk = fkk / np.sqrt(Ny * Nx)
3143 # calculate variance of spectral coefficient
3144 sigma_2 = fkk**2 / Nx / Ny
3146 # Group ellipses of alphas into the same wavenumber k/Nmin
3147 for k in range(1, Nmin + 1):
3148 alpha = k / Nmin
3149 alpha_p1 = (k + 1) / Nmin
3151 # Sum up elements matching k
3152 mask_k = np.where((alpha_matrix >= alpha) & (alpha_matrix < alpha_p1))
3153 ps_array[t, k - 1] = np.sum(sigma_2[mask_k])
3155 return ps_array
3158def _create_alpha_matrix(Ny, Nx):
3159 """Construct an array of 2D wavenumbers from 2D wavenumber pair.
3161 Parameters
3162 ----------
3163 Ny, Nx:
3164 Dimensions of the 2D field for which the power spectra is calculated. Used to
3165 create the array of 2D wavenumbers. Each Ny, Nx pair is associated with a
3166 single-scale parameter.
3168 Returns: alpha_matrix
3169 normalisation of 2D wavenumber axes, transforming the spectral domain into
3170 an elliptic coordinate system.
3172 """
3173 # Create x_indices: each row is [1, 2, ..., Nx]
3174 x_indices = np.tile(np.arange(1, Nx + 1), (Ny, 1))
3176 # Create y_indices: each column is [1, 2, ..., Ny]
3177 y_indices = np.tile(np.arange(1, Ny + 1).reshape(Ny, 1), (1, Nx))
3179 # Compute alpha_matrix
3180 alpha_matrix = np.sqrt((x_indices**2) / Nx**2 + (y_indices**2) / Ny**2)
3182 return alpha_matrix