Shapefile masking

docs/src/userguide/plotting_examples/masking_brazil_plot.py

@@ -0,0 +1,24 @@
+"""Global cube masked to Brazil and plotted with quickplot."""
+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_shapefile
+
+country_shp_reader = shpreader.Reader(
+    shpreader.natural_earth(
+        resolution="110m", category="cultural", name="admin_0_countries"
+    )
+)
+brazil_shp = [
+    country.geometry
+    for country in country_shp_reader.records()
+    if "Brazil" in country.attributes["NAME_LONG"]
+][0]
+
+cube = iris.load_cube(iris.sample_data_path("air_temp.pp"))
+brazil_cube = mask_cube_from_shapefile(cube, brazil_shp)
+
+qplt.pcolormesh(brazil_cube)
+plt.show()

docs/src/userguide/subsetting_a_cube.rst

@@ -5,7 +5,8 @@ Subsetting a Cube
 =================
 
 The :doc:`loading_iris_cubes` section of the user guide showed how to load data into multidimensional Iris cubes.
-However it is often necessary to reduce the dimensionality of a cube down to something more appropriate and/or manageable.
+However it is often necessary to reduce the dimensionality of a cube down to something more appropriate and/or manageable,
+or only examine and analyse a subset of data in a dimension.
 
 Iris provides several ways of reducing both the amount of data and/or the number of dimensions in your cube depending on the circumstance.
 In all cases **the subset of a valid cube is itself a valid cube**.
@@ -329,6 +330,36 @@ 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
+--------------
+
+.. _masking-from-shapefile:
+
+Masking from a shapefile
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+Often we want to perform so kind of analysis over a complex geographical feature - only over land points or sea points:
+or over a continent, a country, a river watershed or administrative region. These geographical features can often be described by shapefiles.
+Shapefiles are a file format first developed for GIS software in the 1990s, and now `Natural Earth`_ maintain a large freely usable database of shapefiles of many geographical and poltical divisions,
+accessible via cartopy. Users may also provide their own custom shapefiles.
+
+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.
+
+First, we load the correct shapefile from NaturalEarth via the `Cartopy`_ instructions. Here we get one for Brazil.
+The `.geometry` attribute of the records in the reader contain the shapely polygon we're interested in - once we have those we just need to provide them to
+the :class:`iris.util.mask_cube_from_shapefile` function. Once plotted, we can see that only our area of interest remains in the data.
+
+
+.. plot:: userguide/plotting_examples/masking_brazil_plot.py
+   :include-source:
+
+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 is masked.
+
+.. 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 from this function.
+
+
 Cube Iteration
 --------------
 It is not possible to directly iterate over an Iris cube. That is, you cannot use code such as
@@ -440,3 +471,7 @@ Similarly, Iris cubes have indexing capability::
 	# Get the second element of the first dimension and all of the second dimension
 	# in reverse, by steps of two.
 	print(cube[1, ::-2])
+
+
+.. _Cartopy: https://scitools.org.uk/cartopy/docs/latest/tutorials/using_the_shapereader.html#id1
+.. _Natural Earth: https://www.naturalearthdata.com/

docs/src/whatsnew/latest.rst

@@ -67,6 +67,12 @@ This document explains the changes made to Iris for this release
    NetCDF chunking with the use of the :data:`iris.fileformats.netcdf.loader.CHUNK_CONTROL`
    context manager. (:pull:`5588`)
 
+#. `@acchamber`_ and `@trexfeathers`_ (reviewer) added
+   :func:`iris.util.mask_cube_from_shapefile`. This builds on the original work
+   of `@ckmo`_, `@david-bentley`_, `@jmendesmetoffice`_, `@evyve`_ and
+   `@pelson`_ for the UK Met Office **ASCEND** library. See
+   :ref:`masking-from-shapefile` for documentation. (:pull:`5470`)
+
 
 🐛 Bugs Fixed
 =============
@@ -110,7 +116,7 @@ This document explains the changes made to Iris for this release
 #. `@bouweandela`_ changed :func:`iris.coords.Coord.cell` so it does not realize
    all coordinate data and only loads a single cell instead. (:pull:`5693`)
 
-#. `@rcomer`_ and `@trexfeathers`_ (reviewer) modified 
+#. `@rcomer`_ and `@trexfeathers`_ (reviewer) modified
    :func:`~iris.analysis.stats.pearsonr` so it preserves lazy data in all cases
    and also runs a little faster.  (:pull:`5638`)
 
@@ -242,6 +248,10 @@ This document explains the changes made to Iris for this release
 .. _@scottrobinson02: https://github.com/scottrobinson02
 .. _@acchamber: https://github.com/acchamber
 .. _@fazledyn-or: https://github.com/fazledyn-or
+.. _@ckmo: https://github.com/ckmo
+.. _@david-bentley: https://github.com/david-bentley
+.. _@jmendesmetoffice: https://github.com/jmendesmetoffice
+.. _@evyve: https://github.com/evyve
 
 
 .. comment

lib/iris/_shapefiles.py

@@ -0,0 +1,243 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+
+# Much of this code is originally based off the ASCEND library, developed in
+# the Met Office by Chris Kent, Emilie Vanvyve, David Bentley, Joana Mendes
+# many thanks to them. Converted to iris by Alex Chamberlain-Clay
+
+
+from itertools import product
+import warnings
+
+import numpy as np
+import shapely
+import shapely.errors
+import shapely.geometry as sgeom
+import shapely.ops
+
+from iris.exceptions import IrisDefaultingWarning, IrisUserWarning
+
+
+def create_shapefile_mask(
+    geometry,
+    cube,
+    minimum_weight=0.0,
+):
+    """Make a mask for a cube from a shape.
+
+    Get the mask of the intersection between the
+    given shapely geometry and cube with x/y DimCoords.
+    Can take a minimum weight and evaluate area overlaps instead
+
+    Parameters
+    ----------
+    geometry : :class:`shapely.Geometry`
+    cube : :class:`iris.cube.Cube`
+        A :class:`~iris.cube.Cube` which has 1d x and y coordinates
+    minimum_weight : float, default 0.0
+        A float between 0 and 1 determining what % of a cell
+        a shape must cover for the cell to remain unmasked.
+        eg: 0.1 means that at least 10% of the shape overlaps the cell
+        to be unmasked.
+        Requires geometry to be a Polygon or MultiPolygon
+        Defaults to 0.0 (eg only test intersection)
+
+    Returns
+    -------
+    :class:`np.array`
+        An array of the shape of the x & y coordinates of the cube, with points
+        to mask equal to True
+
+    """
+    from iris.cube import Cube, CubeList
+
+    try:
+        msg = "Geometry is not a valid Shapely object"
+        if not shapely.is_valid(geometry):
+            raise TypeError(msg)
+    except Exception:
+        raise TypeError(msg)
+    if not isinstance(cube, Cube):
+        if isinstance(cube, CubeList):
+            msg = "Received CubeList object rather than Cube - \
+            to mask a CubeList iterate over each Cube"
+            raise TypeError(msg)
+        else:
+            msg = "Received non-Cube object where a Cube is expected"
+            raise TypeError(msg)
+    if minimum_weight > 0.0 and isinstance(
+        geometry,
+        (
+            sgeom.Point,
+            sgeom.LineString,
+            sgeom.LinearRing,
+            sgeom.MultiPoint,
+            sgeom.MultiLineString,
+        ),
+    ):
+        minimum_weight = 0.0
+        warnings.warn(
+            """Shape is of invalid type for minimum weight masking,
+            must use a Polygon rather than Line shape.\n
+              Masking based off intersection instead. """,
+            category=IrisDefaultingWarning,
+        )
+
+    # prepare 2D cube
+    y_name, x_name = _cube_primary_xy_coord_names(cube)
+    cube_2d = cube.slices([y_name, x_name]).next()
+    for coord in cube_2d.dim_coords:
+        if not coord.has_bounds():
+            coord.guess_bounds()
+    trans_geo = _transform_coord_system(geometry, cube_2d)
+
+    y_coord, x_coord = [cube_2d.coord(n) for n in (y_name, x_name)]
+    x_bounds = _get_mod_rebased_coord_bounds(x_coord)
+    y_bounds = _get_mod_rebased_coord_bounds(y_coord)
+    # prepare array for dark
+    box_template = [
+        sgeom.box(x[0], y[0], x[1], y[1]) for x, y in product(x_bounds, y_bounds)
+    ]
+    # shapely can do lazy evaluation of intersections if it's given a list of grid box shapes
+    # delayed lets us do it in parallel
+    intersect_template = shapely.intersects(trans_geo, box_template)
+    # we want areas not under shapefile to be True (to mask)
+    intersect_template = np.invert(intersect_template)
+    # now calc area overlaps if doing weights and adjust mask
+    if minimum_weight > 0.0:
+        intersections = np.array(box_template)[~intersect_template]
+        intersect_template[~intersect_template] = [
+            trans_geo.intersection(box).area / box.area <= minimum_weight
+            for box in intersections
+        ]
+    mask_template = np.reshape(intersect_template, cube_2d.shape[::-1]).T
+    return mask_template
+
+
+def _transform_coord_system(geometry, cube, geometry_system=None):
+    """Project the shape onto another coordinate system.
+
+    Parameters
+    ----------
+    geometry: :class:`shapely.Geometry`
+    cube: :class:`iris.cube.Cube`
+        :class:`~iris.cube.Cube` with the coord_system to be projected to and
+        a x coordinate
+    geometry_system: :class:`iris.coord_systems`, optional
+        A :class:`~iris.coord_systems` object describing
+        the coord_system of the shapefile. Defaults to None,
+        which is treated as GeogCS
+
+    Returns
+    -------
+    :class:`shapely.Geometry`
+        A transformed copy of the provided :class:`shapely.Geometry`
+
+    """
+    y_name, x_name = _cube_primary_xy_coord_names(cube)
+    import iris.analysis.cartography
+
+    DEFAULT_CS = iris.coord_systems.GeogCS(
+        iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
+    )
+    target_system = cube.coord_system()
+    if not target_system:
+        warnings.warn(
+            "Cube has no coord_system; using default GeogCS lat/lon",
+            category=IrisDefaultingWarning,
+        )
+        target_system = DEFAULT_CS
+    if geometry_system is None:
+        geometry_system = DEFAULT_CS
+    target_proj = target_system.as_cartopy_projection()
+    source_proj = geometry_system.as_cartopy_projection()
+
+    trans_geometry = target_proj.project_geometry(geometry, source_proj)
+    # A GeogCS in iris can be either -180 to 180 or 0 to 360. If cube is 0-360, shift geom to match
+    if (
+        isinstance(target_system, iris.coord_systems.GeogCS)
+        and cube.coord(x_name).points[-1] > 180
+    ):
+        # chop geom at 0 degree line very finely then transform
+        prime_meridian_line = shapely.LineString([(0, 90), (0, -90)])
+        trans_geometry = trans_geometry.difference(prime_meridian_line.buffer(0.00001))
+        trans_geometry = shapely.transform(trans_geometry, _trans_func)
+
+    if (not isinstance(target_system, iris.coord_systems.GeogCS)) and cube.coord(
+        x_name
+    ).points[-1] > 180:
+        # this may lead to incorrect masking or not depending on projection type so warn user
+        warnings.warn(
+            """Cube has x-coordinates over 180E and a non-standard projection type.\n
+              This may lead to incorrect masking. \n
+             If the result is not as expected, you might want to transform the x coordinate points of your cube to -180-180 """,
+            category=IrisUserWarning,
+        )
+    return trans_geometry
+
+
+def _trans_func(geometry):
+    """Pocket function for transforming the x coord of a geometry from -180 to 180 to 0-360."""
+    for point in geometry:
+        if point[0] < 0:
+            point[0] = 360 - np.abs(point[0])
+    return geometry
+
+
+def _cube_primary_xy_coord_names(cube):
+    """Return the primary latitude and longitude coordinate names, or long names, from a cube.
+
+    Parameters
+    ----------
+    cube : :class:`iris.cube.Cube`
+
+    Returns
+    -------
+    tuple of str
+        The names of the primary latitude and longitude coordinates
+
+    """
+    latc = (
+        cube.coords(axis="y", dim_coords=True)[0]
+        if cube.coords(axis="y", dim_coords=True)
+        else -1
+    )
+    lonc = (
+        cube.coords(axis="x", dim_coords=True)[0]
+        if cube.coords(axis="x", dim_coords=True)
+        else -1
+    )
+
+    if -1 in (latc, lonc):
+        msg = "Error retrieving 1d xy coordinates in cube: {!r}"
+        raise ValueError(msg.format(cube))
+
+    latitude = latc.name()
+    longitude = lonc.name()
+    return latitude, longitude
+
+
+def _get_mod_rebased_coord_bounds(coord):
+    """Take in a coord and returns a array of the bounds of that coord rebased to the modulus.
+
+    Parameters
+    ----------
+    coord : :class:`iris.coords.Coord`
+        An Iris coordinate with a modulus
+
+    Returns
+    -------
+    :class:`np.array`
+        A 1d Numpy array of [start,end] pairs for bounds of the coord
+
+    """
+    modulus = coord.units.modulus
+    # Force realisation (rather than core_bounds) - more efficient for the
+    #  repeated indexing happening downstream.
+    result = np.array(coord.bounds)
+    if modulus:
+        result[result < 0.0] = (np.abs(result[result < 0.0]) % modulus) * -1
+        result[np.isclose(result, modulus, 1e-10)] = 0.0
+    return result