Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/src/userguide/plotting_examples/masking_brazil_plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -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()
53 changes: 53 additions & 0 deletions docs/src/userguide/plotting_examples/masking_stereographic_plot.py
Original file line number Diff line number Diff line change
@@ -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()
128 changes: 114 additions & 14 deletions docs/src/userguide/subsetting_a_cube.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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


11 changes: 6 additions & 5 deletions docs/src/whatsnew/3.14.rst
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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`)


Expand Down
28 changes: 17 additions & 11 deletions lib/iris/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
-----
Expand All @@ -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 "
Expand Down Expand Up @@ -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``.
Expand Down Expand Up @@ -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.
Expand Down
Loading