Mask to coast

%matplotlib notebook
from cem_gridtools.cem_gridtools import Coast,CGridRect
import matplotlib.pyplot as plt
from shapely.geometry import MultiPolygon
import numpy as np

x = np.array([147.2599, 147.2938, 147.2990, 147.2651])
y = np.array([-43.0557, -43.0657, -43.0480, -43.0380])

CG = CGridRect(x,y,True)
CG.gen_gridrect(200)
CG.plot()

# Import and plot coastline on the same plot as above
fname = 'setas.cst'
C = Coast(fname)
ax = plt.gca()
C.plot(ax)

Polygon using indices (2279:2283) is invalid:Too few points in geometry component[145.299916 -42.180851]
Polygon using indices (16483:16487) is invalid:Too few points in geometry component[146.348145 -43.58688]
Polygon using indices (18964:18986) is invalid:Self-intersection[147.976090376623 -42.8566280649351]
Polygon using indices (18987:18991) is invalid:Self-intersection[148.045237 -42.653316]
Polygon using indices (76783:77384) is invalid:Self-intersection[146.347855 -43.58678]
Polygon using indices (77385:222908) is invalid:Self-intersection[145.543558333333 -42.3834755555556]
Polygon using indices (222909:223059) is invalid:Self-intersection[145.984111481604 -43.01193]

# Zoom into grid
bd=CG.get_grid_as_lines().buffer(0.01).bounds
ax.set_xlim([bd[0], bd[2]])
ax.set_ylim([bd[1], bd[3]])
ax.grid(True)

# We can now create a MultiPolygon geometry to help with the masking
polyFunc = C.get_polyons_as_generator()
p = next(polyFunc)
polys = []
while p is not None:
    polys.append(p)
    p = next(polyFunc)
# Looks like we still need a bit of a clean-up
MP = MultiPolygon(polys).buffer(0)

Polygon using indices (2279:2283) is invalid:Too few points in geometry component[145.299916 -42.180851]
Polygon using indices (16483:16487) is invalid:Too few points in geometry component[146.348145 -43.58688]
Polygon using indices (18964:18986) is invalid:Self-intersection[147.976090376623 -42.8566280649351]
Polygon using indices (18987:18991) is invalid:Self-intersection[148.045237 -42.653316]
Polygon using indices (76783:77384) is invalid:Self-intersection[146.347855 -43.58678]
Polygon using indices (77385:222908) is invalid:Self-intersection[145.543558333333 -42.3834755555556]
Polygon using indices (222909:223059) is invalid:Self-intersection[145.984111481604 -43.01193]

Ways to mask LAND cells

The easiest would be something like LandCells = GridMultiPolygon.intersection(CoastMultiPolygon). However, MultiPolygon is not valid for grids due to overlaps at edges so they have to be represented as a list of Polygons. In this case we have to loop over each grid cell and ask check whether it intersects with the coastline as follows:

# Check intersection and plot an X in the above plot
for grid_poly in CG.get_grid_as_polygons():
    if MP.intersects(grid_poly):
        c = grid_poly.centroid.coords[:][0]
        ax.plot(c[0],c[1],'kx')

Another option is to use the MultiPoint set of the grid cell centres. This has the performance and concise syntax advantage but will miss any intersections of coastlines that do not encompass the cell centres

land_points = CG.get_cell_centres_as_points().intersection(MP)

land_points

Other options could be, percentage overlap, crossings of the grid lines …