Gallery: update seasonal ensemble example

docs/iris/gallery_code/meteorology/plot_lagged_ensemble.py

@@ -19,6 +19,7 @@
 """
 
 import matplotlib.pyplot as plt
+import matplotlib.ticker
 import numpy as np
 
 import iris
@@ -32,45 +33,48 @@ def realization_metadata(cube, field, fname):
     in the cube.
 
     """
-    # add an ensemble member coordinate if one doesn't already exist
+    # Add an ensemble member coordinate if one doesn't already exist.
     if not cube.coords("realization"):
-        # the ensemble member is encoded in the filename as *_???.pp where ???
-        # is the ensemble member
+        # The ensemble member is encoded in the filename as *_???.pp where ???
+        # is the ensemble member.
         realization_number = fname[-6:-3]
-
-        import iris.coords
-
         realization_coord = iris.coords.AuxCoord(
             np.int32(realization_number), "realization", units="1"
         )
         cube.add_aux_coord(realization_coord)
 
 
 def main():
-    # extract surface temperature cubes which have an ensemble member
-    # coordinate, adding appropriate lagged ensemble metadata
+    # Create a constraint to extract surface temperature cubes which have a
+    # "realization" coordinate.
+    constraint = iris.Constraint(
+        "surface_temperature", realization=lambda value: True
+    )
+    # Use this to load our ensemble.  The callback ensures all our members
+    # have the "realization" coordinate and therefore they will all be loaded.
     surface_temp = iris.load_cube(
         iris.sample_data_path("GloSea4", "ensemble_???.pp"),
-        iris.Constraint("surface_temperature", realization=lambda value: True),
+        constraint,
         callback=realization_metadata,
     )
 
     # -------------------------------------------------------------------------
     # Plot #1: Ensemble postage stamps
     # -------------------------------------------------------------------------
 
-    # for the purposes of this example, take the last time element of the cube
-    last_timestep = surface_temp[:, -1, :, :]
+    # For the purposes of this example, take the last time element of the cube.
+    # First get hold of the last time by slicing the coordinate.
+    last_time_coord = surface_temp.coord("time")[-1]
+    last_timestep = surface_temp.subset(last_time_coord)
 
-    # Make 50 evenly spaced levels which span the dataset
-    contour_levels = np.linspace(
-        np.min(last_timestep.data), np.max(last_timestep.data), 50
-    )
+    # Find the maximum and minimum across the dataset.
+    data_min = np.min(last_timestep.data)
+    data_max = np.max(last_timestep.data)
 
-    # Create a wider than normal figure to support our many plots
+    # Create a wider than normal figure to support our many plots.
     plt.figure(figsize=(12, 6), dpi=100)
 
-    # Also manually adjust the spacings which are used when creating subplots
+    # Also manually adjust the spacings which are used when creating subplots.
     plt.gcf().subplots_adjust(
         hspace=0.05,
         wspace=0.05,
@@ -80,83 +84,68 @@ def main():
         right=0.925,
     )
 
-    # iterate over all possible latitude longitude slices
+    # Iterate over all possible latitude longitude slices.
     for cube in last_timestep.slices(["latitude", "longitude"]):
 
-        # get the ensemble member number from the ensemble coordinate
+        # Get the ensemble member number from the ensemble coordinate.
         ens_member = cube.coord("realization").points[0]
 
-        # plot the data in a 4x4 grid, with each plot's position in the grid
-        # being determined by ensemble member number the special case for the
-        # 13th ensemble member is to have the plot at the bottom right
+        # Plot the data in a 4x4 grid, with each plot's position in the grid
+        # being determined by ensemble member number.  The special case for the
+        # 13th ensemble member is to have the plot at the bottom right.
         if ens_member == 13:
             plt.subplot(4, 4, 16)
         else:
             plt.subplot(4, 4, ens_member + 1)
 
-        cf = iplt.contourf(cube, contour_levels)
+        # Plot with 50 evenly spaced contour levels (49 intervals).
+        cf = iplt.contourf(cube, 49, vmin=data_min, vmax=data_max)
 
-        # add coastlines
+        # Add coastlines.
         plt.gca().coastlines()
 
-    # make an axes to put the shared colorbar in
+    # Make an axes to put the shared colorbar in.
     colorbar_axes = plt.gcf().add_axes([0.35, 0.1, 0.3, 0.05])
     colorbar = plt.colorbar(cf, colorbar_axes, orientation="horizontal")
-    colorbar.set_label("%s" % last_timestep.units)
-
-    # limit the colorbar to 8 tick marks
-    import matplotlib.ticker
+    colorbar.set_label(last_timestep.units)
 
+    # Limit the colorbar to 8 tick marks.
     colorbar.locator = matplotlib.ticker.MaxNLocator(8)
     colorbar.update_ticks()
 
-    # get the time for the entire plot
-    time_coord = last_timestep.coord("time")
-    time = time_coord.units.num2date(time_coord.bounds[0, 0])
+    # Get the time for the entire plot.
+    time = last_time_coord.units.num2date(last_time_coord.bounds[0, 0])
 
-    # set a global title for the postage stamps with the date formated by
-    # "monthname year"
-    plt.suptitle(
-        "Surface temperature ensemble forecasts for %s"
-        % (time.strftime("%B %Y"),)
-    )
+    # Set a global title for the postage stamps with the date formated by
+    # "monthname year".
+    time_string = time.strftime("%B %Y")
+    plt.suptitle(f"Surface temperature ensemble forecasts for {time_string}")
 
     iplt.show()
 
     # -------------------------------------------------------------------------
     # Plot #2: ENSO plumes
     # -------------------------------------------------------------------------
 
-    # Nino 3.4 lies between: 170W and 120W, 5N and 5S, so define a constraint
-    # which matches this
-    nino_3_4_constraint = iris.Constraint(
-        longitude=lambda v: -170 + 360 <= v <= -120 + 360,
-        latitude=lambda v: -5 <= v <= 5,
+    # Nino 3.4 lies between: 170W and 120W, 5N and 5S, so use the intersection
+    # method to restrict to this region.
+    nino_cube = surface_temp.intersection(
+        latitude=[-5, 5], longitude=[-170, -120]
     )
 
-    nino_cube = surface_temp.extract(nino_3_4_constraint)
-
-    # Subsetting a circular longitude coordinate always results in a circular
-    # coordinate, so set the coordinate to be non-circular
-    nino_cube.coord("longitude").circular = False
-
-    # Calculate the horizontal mean for the nino region
+    # Calculate the horizontal mean for the nino region.
     mean = nino_cube.collapsed(["latitude", "longitude"], iris.analysis.MEAN)
 
-    # Calculate the ensemble mean of the horizontal mean. To do this, remove
-    # the "forecast_period" and "forecast_reference_time" coordinates which
-    # span both "relalization" and "time".
-    mean.remove_coord("forecast_reference_time")
-    mean.remove_coord("forecast_period")
+    # Calculate the ensemble mean of the horizontal mean.
     ensemble_mean = mean.collapsed("realization", iris.analysis.MEAN)
 
-    # take the ensemble mean from each ensemble member
-    mean -= ensemble_mean.data
+    # Take the ensemble mean from each ensemble member.
+    mean -= ensemble_mean
 
     plt.figure()
 
     for ensemble_member in mean.slices(["time"]):
-        # draw each ensemble member as a dashed line in black
+        # Draw each ensemble member as a dashed line in black.
         iplt.plot(ensemble_member, "--k")
 
     plt.title("Mean temperature anomaly for ENSO 3.4 region")

docs/iris/src/whatsnew/latest.rst

@@ -47,7 +47,8 @@ This document explains the changes made to Iris for this release
 📚 Documentation
 ================
 
-* N/A
+* `@rcomer`_ updated the "Seasonal ensemble model plots" Gallery example.
+  (:pull:`3933`)
 
 
 💼 Internal
@@ -56,3 +57,4 @@ This document explains the changes made to Iris for this release
 * N/A
 
 .. _@gcaria: https://github.com/gcaria
+.. _@rcomer: https://github.com/rcomer