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,

```
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)
```

and loading the `lon`

, `lat`

values into a 2D grid.

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

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

```
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)
```

```
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:

```
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:

```
Points = [Point(xp, yp) for xp, yp in zip(x.ravel(), y.ravel())]
Points = cascaded_union(Points)
geom = Points.intersection(polygon)
geom
```

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.

```
geom.xy # :-(
```

But this works...

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

And finally our masked out figure!

```
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)
```

```
HTML(html)
```