Mask cube with features

Question

I want to plot data from a global cube, but only for a list of countries. So I select a subcube according to the countries' "bounding box".

So far so good. What I'm looking for is an easy way to mask out all points of a cube which do not fall in any of my countries (which are represented as features), so that only those points of the cube which lie within any of my features are plotted.

Any idea is greatly appreciated =)

Answer 1

You can achieve this directly at the plotting stage rather than masking the cube within iris. I've approached this by setting the clip path of the artist returned by pcolor. The method is to create a list of geometries from features (in this case countries from Natural Earth, they could be from a shapefile) then transform these geometries into a matplotlib path which the image can be clipped to. I'll detail this method, and hopefully this will be enough to get you started:

I first defined a function to retrieve the Shapely geometries corresponding to given country names, the geometries come from the Natural Earth 110m administrative boundaries shapefile, access through the cartopy interface.

I then defined a second function which is a wrapper around the iris.plot.pcolor function which makes the plot and clips it to the given geometries.

Now all I need to do is set up the plot as normal, but use the plotting wrapper instead of directly calling the iris.plot.pcolor function.

Here is a complete example:

import cartopy.crs as ccrs
from cartopy.io.shapereader import natural_earth, Reader
from cartopy.mpl.patch import geos_to_path
import iris
import iris.plot as iplt
import matplotlib.pyplot as plt
from matplotlib.path import Path


def get_geometries(country_names):
    """
    Get an iterable of Shapely geometries corrresponding to given countries.

    """
    # Using the Natural Earth feature interface provided by cartopy.
    # You could use a different source, all you need is the geometries.
    shape_records = Reader(natural_earth(resolution='110m',
                                         category='cultural',
                                         name='admin_0_countries')).records()
    geoms = []
    for country in shape_records:
        if country.attributes['name_long'] in country_names:
            try:
                geoms += country.geometry
            except TypeError:
                geoms.append(country.geometry)
    return geoms, ccrs.PlateCarree()._as_mpl_transform


def pcolor_mask_geoms(cube, geoms, transform):
    path = Path.make_compound_path(*geos_to_path(geoms))
    im = iplt.pcolor(cube)
    im.set_clip_path(path, transform=transform)


# First plot the full map:
cube = iris.load_cube(iris.sample_data_path('air_temp.pp'))
plt.figure(figsize=(12, 6))
ax1 = plt.axes(projection=ccrs.PlateCarree())
ax1.coastlines()
iplt.pcolor(cube)

# Now plot just the required countries:
plt.figure(figsize=(12, 6))
ax2 = plt.axes(projection=ccrs.PlateCarree())
ax2.coastlines()
countries = [
    'United States',
    'United Kingdom',
    'Saudi Arabia',
    'South Africa',
    'Nigeria']
geoms, transform = get_geometries(countries)
pcolor_mask_geoms(cube, geoms, transform(ax2))

plt.show()

The results of which look like this:

If you want to use iris.plot.pcolormesh instead you will need to modify the plotting function a little bit. This is dues to a workaround for a matplotlib issue that is currently included in cartopy. The modified version would look like this:

def pcolor_mask_geoms(cube, geoms, transform):
    path = Path.make_compound_path(*geos_to_path(geoms))
    im = iplt.pcolormesh(cube)
    im.set_clip_path(path, transform=transform)
    try:
        im._wrapped_collection_fix.set_clip_path(path, transform)
    except AttributeError:
        pass