Masking land/ocean with shapely

A common task when plotting data on a map (or reduce the number of points in KDTree search ;-) is to mask either the land or the oceanic part of the data. In this post I will show how to do mask land using a shapefile and shapely.

Let's start by loading some global data,

In [3]:
import warnings
import iris

url = "data/ERAI.EP.1979-2011.nc"

# Suppress iris warnings.
with warnings.catch_warnings():
warnings.simplefilter("ignore")
cube = iris.load_raw(url)[0]

print(cube)

Budget Evaporation minus Precipitation / (kg/m2/s) (time: 396; latitude: 256; longitude: 512)
Dimension coordinates:
time                                          x              -               -
latitude                                      -              x               -
longitude                                     -              -               x
Attributes:
Acknowledgment: -
comment: -
contact: -
conventions: -
creation_date: Thu Feb  2 17:09:42 2012
delta_time: monthly
history: by John Fasullo, Thu Feb  2 17:09:42 2012
institution: -
product: -
references: -



and loading the lon, lat values into a 2D grid.

In [4]:
import numpy as np
from oceans import wrap_lon180

lon = wrap_lon180(cube.coord(axis='X').points)
lat = cube.coord(axis='Y').points

x, y = np.meshgrid(lon, lat)


Now we need a shapefile to use as mask. I will use the a map of Brazil that I fused all the polygons, with cascaded_union, to get a simple shape.

In [5]:
from cartopy.io import shapereader
from shapely.ops import cascaded_union

shp = shapereader.Reader('./data/BRA/BRA_adm2')

geoms = shp.geometries()
polygon = cascaded_union(list(geoms))


The last piece we need is a function that checks if a points is inside the polygon. This is very slow! I know that there are some faster way to do this out there, like matplotlib's path, but we will try those in another post.

In [6]:
from shapely.geometry import Point

def inpolygon(polygon, xp, yp):
return np.array([Point(x, y).intersects(polygon) for x, y in zip(xp, yp)],
dtype=np.bool)

In [7]:
mask = inpolygon(polygon, x.ravel(), y.ravel())


The inpolygon() returns a mask of the same size as the points. We can use it to plot all but the area inside the polygon. For example:

In [8]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

m = mask.reshape(x.shape)

dx = dy = 10
mm = m[::dx, ::dy]
xx = x[::dx, ::dy]
yy = y[::dx, ::dy]

projection = ccrs.PlateCarree()
fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
ax.plot(xx[~mm], yy[~mm], 'k.', alpha=0.25)
c = ax.coastlines()


Nice! With global shapefile we could have masked all the land and plotted only the data over the ocean. Another way to do this is:

In [9]:
Points = [Point(xp, yp) for xp, yp in zip(x.ravel(), y.ravel())]

Points = cascaded_union(Points)

geom = Points.intersection(polygon)
geom

Out[9]:

This is much faster and we get a shapely object back! So there are a few more thing we can do with that object. One thing we cannot do though is to get the points without a looping over it.

In [10]:
geom.xy  # :-(

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-10-d1c0d04ceb59> in <module>()
----> 1 geom.xy

/home/filipe/.virtualenvs/blog/lib/python2.7/site-packages/shapely/geometry/base.py in xy(self)
310     def xy(self):
311         """Separate arrays of X and Y coordinate values"""
--> 312         raise NotImplementedError
313
314     # Python feature protocol

NotImplementedError: 

But this works...

In [11]:
np.array([(points.x, points.y) for points in geom.geoms])[:5]

Out[11]:
array([[-73.82723999,  -7.36840677],
[-73.12411499,  -8.77191257],
[-73.12411499,  -8.07015991],
[-73.12411499,  -7.36840677],
[-73.12411499,  -6.66665411]])


And finally our masked out figure!

In [12]:
import numpy.ma as ma

arr = cube[0, ...].data

fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
cs = ax.pcolormesh(lon, lat, ma.masked_array(arr, ~m))
ax.set_extent(geom.bounds)

In [13]:
HTML(html)

Out[13]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.