Search code examples
pythonnumpymatplotlibmaskcartopy

How to add an outline to a mask with no internal lines in matplotlib


For a given 2D mask, I would like to draw its outline with no internal lines between adjacent grid cells. Highlighting one cell is straightforward: How to mark cells in matplotlib.pyplot.imshow (drawing cell borders) and highlighting many cells is also straightforward: Python matplotlib - add borders to grid plot based on value. However, I do not want internal boundaries within the mask to be separated by a line, so the second link above is not suitable.

Here is example code to generate a mask:

import matplotlib.pyplot as plt
import numpy as np

N = 50
xbound = np.linspace(0, 10, N + 1)
ybound = np.linspace(0, 10, N + 1)
x = (xbound[:-1] + xbound[1:]) / 2
y = (ybound[:-1] + ybound[1:]) / 2

X, Y = np.meshgrid(x, y)
Z = np.exp(-((X - 2.5)**2 + (Y - 2.5)**2)) + np.exp(-2 * ((X - 7.5)**2 + (Y - 6.5)**2))
mask = Z > 0.2

plt.imshow(mask, origin='lower', extent=(0, 10, 0, 10))

Which generates the image:

Simple 2D mask

I want the boundary to be around each of the pixelated yellow circles, and to follow the same edges as show above (i.e. be parallel to the x/y-axis) - I want to emphasize that the underlying data is gridded. Using:

plt.contour(x, y, mask, levels=[0.5])

comes close, but the contour is at 45° along the staircase edges.

Bonus points if the outline can be filled or shown using cartopy.


Solution

  • This can be done using shapely shapes and the union operation. The idea is to create a square (or rectangle etc.) for each cell where the mask is true, then combine contiguous squares together. These geometries can then be plotted or used to create matplotlib polygons. The following shows the outline in red (with the staircase edge as required):

    import shapely.geometry
    import shapely.ops
    
    geoms = []
    for yidx, xidx in zip(*np.where(mask)):
        geoms.append(shapely.geometry.box(xbound[xidx], ybound[yidx], xbound[xidx+1], ybound[yidx+1]))
    full_geom = shapely.ops.unary_union(geoms)
    
    for geom in full_geom.geoms:
        plt.plot(*geom.exterior.xy, linewidth=4, color='r')
    plt.imshow(mask, origin='lower', extent=(0, 10, 0, 10))
    

    mask with plotted outline

    It would be hard to fill each shape with a colour/pattern, so instead matplotlib polygons can be used:

    from matplotlib.patches import Polygon
    
    polygons = []
    for geom in full_geom.geoms:
        p = Polygon(list(zip(*geom.exterior.xy)), closed=True, facecolor='none', edgecolor='r', linestyle='-', hatch='//', linewidth=4)
        polygons.append(p)
    
    plt.figure()
    axes = plt.gca()
    for p in polygons:
        axes.add_patch(p)
    plt.imshow(mask, origin='lower', extent=(0, 10, 0, 10))
    

    mask with outline defined by (fillable) polygons

    Finally, the raw shapely polygons can be added directly to a cartopy map:

    import cartopy.crs as ccrs
    
    fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
    ax.add_geometries(
        full_geom.geoms,
        crs=ccrs.PlateCarree(),
        facecolor='white',
        edgecolor='red',
        linewidth=4,
        hatch='//',
    )
    ax.coastlines()
    ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
    ax.set_xlim((0, 10))
    ax.set_ylim((0, 10))
    

    mask with outline on cartopy map