Merge pull request #6801 from stephenworsley/v3.14.x.mergeback

V3.14.x.mergeback

Author: bjlittle

docs/src/userguide/plotting_examples/masking_brazil_plot.py

@@ -5,7 +5,7 @@
 
 import iris
 import iris.quickplot as qplt
-from iris.util import mask_cube_from_shapefile
+from iris.util import mask_cube_from_shape
 
 country_shp_reader = shpreader.Reader(
     shpreader.natural_earth(
@@ -19,7 +19,7 @@
 ][0]
 
 cube = iris.load_cube(iris.sample_data_path("air_temp.pp"))
-brazil_cube = mask_cube_from_shapefile(cube, brazil_shp)
+brazil_cube = mask_cube_from_shape(cube=cube, shape=brazil_shp)
 
 qplt.pcolormesh(brazil_cube)
 plt.show()

docs/src/userguide/plotting_examples/masking_stereographic_plot.py

@@ -0,0 +1,53 @@
+"""Masking data with a stereographic projection and plotted with quickplot."""
+
+import cartopy.crs as ccrs
+import cartopy.io.shapereader as shpreader
+import matplotlib.pyplot as plt
+
+import iris
+import iris.quickplot as qplt
+from iris.util import mask_cube_from_shape
+
+# Define WGS84 coordinate reference system
+wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(ellipse="WGS84"))
+
+country_shp_reader = shpreader.Reader(
+    shpreader.natural_earth(
+        resolution="10m", category="cultural", name="admin_0_countries"
+    )
+)
+uk_shp = [
+    country.geometry
+    for country in country_shp_reader.records()
+    if "United Kingdom" in country.attributes["NAME_LONG"]
+][0]
+
+cube = iris.load_cube(iris.sample_data_path("toa_brightness_stereographic.nc"))
+uk_cube = mask_cube_from_shape(cube=cube, shape=uk_shp, shape_crs=wgs84)
+
+plt.figure(figsize=(12, 5))
+# Plot #1: original data
+ax = plt.subplot(131)
+qplt.pcolormesh(cube, vmin=210, vmax=330)
+plt.gca().coastlines()
+
+# Plot #2: UK geometry
+ax = plt.subplot(132, title="Mask Geometry", projection=ccrs.Orthographic(-5, 45))
+ax.set_extent([-12, 5, 49, 61])
+ax.add_geometries(
+    [
+        uk_shp,
+    ],
+    crs=wgs84,
+    edgecolor="None",
+    facecolor="orange",
+)
+plt.gca().coastlines()
+
+# Plot #3 masked data
+ax = plt.subplot(133)
+qplt.pcolormesh(uk_cube, vmin=210, vmax=330)
+plt.gca().coastlines()
+
+plt.tight_layout()
+plt.show()

docs/src/userguide/subsetting_a_cube.rst

@@ -332,9 +332,55 @@ on bounds can be done in the following way::
 The above example constrains to cells where either the upper or lower bound occur
 after 1st January 2008.
 
+.. _cube_masking:
+
 Cube Masking
 --------------
 
+Masking a cube allows you to hide unwanted data points without changing the 
+shape or size of the cube.  This can be achieved by two methods:
+
+1. Masking a cube using a boolean mask array via :func:`iris.util.mask_cube`.
+2. Masking a cube using a shapefile via :func:`iris.util.mask_cube_from_shape`.
+
+.. _masking-from-boolean:
+
+Masking a cube using a boolean mask array
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The :func:`iris.util.mask_cube` function allows you to mask unwanted data points
+in a cube using a boolean mask array. The mask array must have the same shape as
+the data array in the cube, with ``True`` values indicating points to be masked.
+
+For example, the mask could be based on a threshold value. In the following
+example we mask all points in a cube where the air potential temperature is
+greater than 290 K.
+
+    >>> filename = iris.sample_data_path('uk_hires.pp')
+    >>> cube_temp = iris.load_cube(filename, 'air_potential_temperature')
+    >>> print(cube_temp.summary(shorten=True))
+    air_potential_temperature / (K)     (time: 3; model_level_number: 7; grid_latitude: 204; grid_longitude: 187)
+    >>> type(cube_temp.data)
+    <class 'numpy.ndarray'>
+
+Note that this example cube has 4 dimensions: time and model level number in
+addition to grid latitude, and grid longitude. The data array associated with
+the cube is a regular :py:class:`numpy.ndarray`.
+
+We can build a boolean mask array by applying a condition to the cube's data
+array:
+
+    >>> mask = cube_temp.data > 290
+    >>> cube_masked = iris.util.mask_cube(cube_temp, mask)
+    >>> print(cube_masked.summary(shorten=True))
+    air_potential_temperature / (K)     (time: 3; model_level_number: 7; grid_latitude: 204; grid_longitude: 187)
+    >>> type(cube_masked.data)
+    <class 'numpy.ma.MaskedArray'>
+
+The masked cube will have the same shape and coordinates as the original cube,
+but the data array now includes an associated boolean mask, and the cube's
+`data` property is now a :py:class:`numpy.ma.MaskedArray`.
+
 .. _masking-from-shapefile:
 
 Masking from a shapefile
@@ -346,26 +392,77 @@ Often we want to perform some kind of analysis over a complex geographical featu
 - over a continent, country, or list of countries
 - over a river watershed or lake basin
 - over states or administrative regions of a country
+- extract data along the trajectory of a storm track
+- extract data at specific points of interest such as cities or weather stations
+
+These geographical features can often be described by `ESRI Shapefiles`_. 
+Shapefiles are a file format first developed for GIS software in the 1990s, and
+`Natural Earth`_ maintain a large freely usable database of shapefiles of many
+geographical and political divisions, accessible via `cartopy`_. Users may also
+provide or create their own custom shapefiles for `cartopy`_ to load, or or any
+other source that can be interpreted as a `shapely.Geometry`_  object, such as
+shapes encoded in a geoJSON or KML file.
+
+The :func:`iris.util.mask_cube_from_shape` function facilitates cube masking
+from shapefiles. Once a shape is loaded as a `shapely.Geometry`, and passed to
+the function, any data outside the bounds of the shape geometry is masked.
+
+.. important::
+    For best masking results, both the cube **and** masking shape (geometry)
+    should have a coordinate reference system (CRS) defined. Note that the CRS of
+    the masking geometry must be provided explicitly to :func:`iris.util.mask_cube_from_shape`
+    (via the ``shape_crs`` keyword argument), whereas the :class:`iris.cube.Cube`
+    CRS is read from the cube itself. 
+    
+    The cube **must** have a :attr:`iris.coords.Coord.coord_system` defined
+    otherwise an error will be raised.
 
-These geographical features can often be described by `ESRI Shapefiles`_. Shapefiles are a file format first developed for GIS software in the 1990s, and `Natural Earth`_ maintain a large freely usable database of shapefiles of many geographical and political divisions,
-accessible via `cartopy`_. Users may also provide their own custom shapefiles for `cartopy`_ to load, or their own underlying geometry in the same format as a shapefile geometry.
+.. note::
+    Because shape vectors are inherently Cartesian in nature, they contain no
+    inherent understanding of the spherical geometry underpinning geographic
+    coordinate systems. For this reason, **shapefiles or shape vectors that
+    cross the antimeridian or poles are not supported by this function** to
+    avoid unexpected masking behaviour.  
+    
+    For shapes that do cross these boundaries, this function expects the user
+    to undertake fixes upstream of Iris, using tools like `GDAL`_ or
+    `antimeridian`_ to ensure correct geometry wrapping.
+
+As an introductory example, we load a shapefile of country borders for Brazil
+from `Natural Earth`_ via the `Cartopy_shapereader`_. The `.geometry` attribute
+of the records in the reader contain the `Shapely`_ polygon we're interested in.
+They contain the coordinates that define the polygon being masked under the
+WGS84 coordinate system. We pass this to the :class:`iris.util.mask_cube_from_shape`
+function and this returns a copy of the cube with a :py:class:`numpy.masked_array`
+as the data payload, where the data outside the shape is hidden by the masked
+array.
 
-These shapefiles can be used to mask an iris cube, so that any data outside the bounds of the shapefile is hidden from further analysis or plotting.
+.. plot:: userguide/plotting_examples/masking_brazil_plot.py
+   :include-source:
 
-First, we load the correct shapefile from NaturalEarth via the `Cartopy_shapereader`_ instructions. Here we get one for Brazil.
-The `.geometry` attribute of the records in the reader contain the `Shapely`_ polygon we're interested in. They contain the coordinates that define the polygon (or set of lines) being masked
-and once we have those we just need to provide them to the :class:`iris.util.mask_cube_from_shapefile` function. 
-This returns a copy of the cube with a :class:`numpy.masked_array` as the data payload, where the data outside the shape is hidden by the masked array. We can see this in the following example.
+We can see that the dimensions of the cube haven't changed - the plot still has
+a global extent. But only the data over Brazil is plotted - the rest has been
+masked out.
 
+.. important::
+    Because we do not explicitly pass a CRS for the shape geometry to 
+    :func:`iris.util.mask_cube_from_shape`, the function assumes the geometry
+    has the same CRS as the cube.
 
-.. plot:: userguide/plotting_examples/masking_brazil_plot.py
-   :include-source:
+However, a :class:`iris.cube.Cube` and `Shapely`_ geometry do not need to have
+the same CRS, as long as both have a CRS defined.  Where the CRS of the 
+:class:`iris.cube.Cube` and geometry differ, :func:`iris.util.mask_cube_from_shape`
+will reproject the geometry (via `GDAL`_) onto the cube's CRS prior to masking.
+The masked cube will be returned in the same CRS as the input cube.
 
-We can see that the dimensions of the cube haven't changed - the plot is still global. But only the data over Brazil is plotted - the rest has been masked out. 
+In the following example, we load a cube containing satellite derived temperature
+data in a stereographic projection (with projected coordinates with units of 
+metres), and mask it to only show data over the United Kingdom, based on a
+shapefile of the UK boundary defined in WGS84 lat-lon coordinates.
+
+.. plot:: userguide/plotting_examples/masking_stereographic_plot.py
+   :include-source:
 
-.. note::
-    While Iris will try to dynamically adjust the shapefile to mask cubes of different projections, it can struggle with rotated pole projections and cubes with Meridians not at 0°
-    Converting your Cube's coordinate system may help if you get a fully masked cube as the output from this function unexpectedly.
 
 
 Cube Iteration
@@ -481,8 +578,11 @@ Similarly, Iris cubes have indexing capability::
 	print(cube[1, ::-2])
 
 
+.. _antimeridian: https://www.gadom.ski/antimeridian/latest/
 .. _Cartopy_shapereader: https://cartopy.readthedocs.io/stable/tutorials/using_the_shapereader.html#id1
-.. _Natural Earth: https://www.naturalearthdata.com/
 .. _ESRI Shapefiles: https://support.esri.com/en-us/technical-paper/esri-shapefile-technical-description-279
+.. _GDAL: https://gdal.org/en/stable/programs/ogr2ogr.html
+.. _Natural Earth: https://www.naturalearthdata.com/
+.. _shapely.Geometry: https://shapely.readthedocs.io/en/stable/geometry.html
 
 

docs/src/whatsnew/3.14.rst

@@ -1,13 +1,13 @@
 .. include:: ../common_links.inc
 
-v3.14 (31 Oct 2025 [release candidate])
-***************************************
+v3.14 (14 Nov 2025)
+*******************
 
 This document explains the changes made to Iris for this release
 (:doc:`View all changes <index>`.)
 
 
-.. dropdown:: 3.14 Release Highlights
+.. dropdown:: v3.14 Release Highlights
    :color: primary
    :icon: info
    :animate: fade-in
@@ -166,8 +166,9 @@ This document explains the changes made to Iris for this release
 #. `@rcomer`_ updated all Cartopy references to point to the new location at
    https://cartopy.readthedocs.io (:pull:`6636`)
 
-#. `@hsteptoe`_ added additional worked examples to the :func:`iris.util.mask_cube_from_shape` 
-   documentation, to demonstrate how to use the function with different types of shapefiles. 
+#. `@hsteptoe`_ added additional worked examples to :ref:`cube_masking` in the user guide, 
+   and :func:`iris.util.mask_cube_from_shape` documentation, to demonstrate how to use the
+   function with different types of shapefiles. 
    (:pull:`6129`)
 
 

lib/iris/util.py

@@ -2262,7 +2262,8 @@ def mask_cube_from_shapefile(
     Parameters
     ----------
     cube : :class:`~iris.cube.Cube` object
-        The ``Cube`` object to masked. Must be singular, rather than a ``CubeList``.
+        The :class:`~iris.cube.Cube` object to masked. Must be singular,
+         rather than a :class:`~iris.cube.CubeList`.
     shape : Shapely.Geometry object
         A single `shape` of the area to remain unmasked on the ``cube``.
         If it a line object of some kind then minimum_weight will be ignored,
@@ -2284,7 +2285,7 @@ def mask_cube_from_shapefile(
     :func:`~iris.util.mask_cube`
         Mask any cells in the cube’s data array.
     :func:`~iris.util.mask_cube_from_shape`
-        Mask any cells in the cube’s data array.
+        Mask all points in a cube that do not intersect a shape object.
 
     Notes
     -----
@@ -2296,7 +2297,9 @@ def mask_cube_from_shapefile(
 
     Warnings
     --------
-    This function requires additional dependencies: ``rasterio`` and ``affine``.
+    This function requires additional dependencies:
+    `rasterio <https://rasterio.readthedocs.io/en/stable/>`_
+    and `affine <https://affine.readthedocs.io/en/latest/>`_.
     """
     message = (
         "iris.util.mask_cube_from_shapefile has been deprecated, and will be removed in a "
@@ -2349,13 +2352,14 @@ def mask_cube_from_shape(
     Parameters
     ----------
     cube : :class:`~iris.cube.Cube` object
-        The ``Cube`` object to masked. Must be singular, rather than a ``CubeList``.
+        The :class:`~iris.cube.Cube` object to masked. Must be singular,
+        rather than a :class:`~iris.cube.CubeList`.
     shape : shapely.Geometry object
         A single ``shape`` of the area to remain unmasked on the ``cube``.
         If it a line object of some kind then minimum_weight will be ignored,
         because you cannot compare the area of a 1D line and 2D Cell.
     shape_crs : cartopy.crs.CRS, default=None
-        The coordinate reference system of the shape object.
+        The coordinate reference system of the ``shape`` object.
     in_place : bool, default=False
         Whether to mask the ``cube`` in-place or return a newly masked ``cube``.
         Defaults to ``False``.
@@ -2440,22 +2444,24 @@ def mask_cube_from_shape(
     Notes
     -----
     Iris does not handle the shape loading so it is agnostic to the source type of the shape.
-    The shape can be loaded from an Esri shapefile, created using the ``shapely`` library, or
-    any other source that can be interpreted as a ``shapely.Geometry`` object, such as shapes
-    encoded in a geoJSON or KML file.
+    The shape can be loaded from an Esri shapefile, created using the
+    `shapely <https://shapely.readthedocs.io/en/stable/>`_ library, or any other source that
+    can be interpreted as a `shapely.Geometry <https://shapely.readthedocs.io/en/stable/geometry.html>`_
+    object, such as shapes encoded in a geoJSON or KML file.
 
     Warnings
     --------
-    For best masking results, both the cube _and_ masking geometry should have a
+    For best masking results, both the cube **and** masking geometry should have a
     coordinate reference system (CRS) defined. Note that CRS of the masking geometry
     must be provided explicitly to this function (via ``shape_crs``), whereas the
     cube CRS is read from the cube itself. The cube **must** have a coord_system defined.
 
     Masking results will be most consistent when the cube and masking geometry have the same CRS.
 
-    If a CRS is _not_ provided for the the masking geometry, the CRS of the cube is assumed.
+    If a CRS is **not** provided for the the masking geometry, the CRS of the cube is assumed.
 
-    This function requires additional dependencies: ``rasterio`` and ``affine``.
+    This function requires additional dependencies: `rasterio <https://rasterio.readthedocs.io/en/stable/>`_
+    and `affine <https://affine.readthedocs.io/en/latest/>`_.
 
     Because shape vectors are inherently Cartesian in nature, they contain no inherent
     understanding of the spherical geometry underpinning geographic coordinate systems.