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/.

python4oceanographers

Turning ripples into waves

Masking land/ocean with shapely

Comments