Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
vector outline to raster mask
import numpy as np
def outline_to_mask(line, x, y):
"""Create mask from outline contour
Parameters
----------
line: array-like (N, 2)
x, y: 1-D grid coordinates (input for meshgrid)
Returns
-------
mask : 2-D boolean array (True inside)
Examples
--------
>>> from shapely.geometry import Point
>>> poly = Point(0,0).buffer(1)
>>> x = np.linspace(-5,5,100)
>>> y = np.linspace(-5,5,100)
>>> mask = outline_to_mask(poly.boundary, x, y)
"""
import matplotlib.path as mplp
mpath = mplp.Path(line)
X, Y = np.meshgrid(x, y)
points = np.array((X.flatten(), Y.flatten())).T
mask = mpath.contains_points(points).reshape(X.shape)
return mask
def _grid_bbox(x, y):
dx = dy = 0
return x[0]-dx/2, x[-1]+dx/2, y[0]-dy/2, y[-1]+dy/2
def _bbox_to_rect(bbox):
l, r, b, t = bbox
return Polygon([(l, b), (r, b), (r, t), (l, t)])
def shp_mask(shp, x, y, m=None):
"""Use recursive sub-division of space and shapely contains method to create a raster mask on a regular grid.
Parameters
----------
shp : shapely's Polygon (or whatever with a "contains" method and intersects method)
x, y : 1-D numpy arrays defining a regular grid
m : mask to fill, optional (will be created otherwise)
Returns
-------
m : boolean 2-D array, True inside shape.
Examples
--------
>>> from shapely.geometry import Point
>>> poly = Point(0,0).buffer(1)
>>> x = np.linspace(-5,5,100)
>>> y = np.linspace(-5,5,100)
>>> mask = shp_mask(poly, x, y)
"""
rect = _bbox_to_rect(_grid_bbox(x, y))
if m is None:
m = np.zeros((y.size, x.size), dtype=bool)
if not shp.intersects(rect):
m[:] = False
elif shp.contains(rect):
m[:] = True
else:
k, l = m.shape
if k == 1 and l == 1:
m[:] = shp.contains(Point(x[0], y[0]))
elif k == 1:
m[:, :l//2] = shp_mask(shp, x[:l//2], y, m[:, :l//2])
m[:, l//2:] = shp_mask(shp, x[l//2:], y, m[:, l//2:])
elif l == 1:
m[:k//2] = shp_mask(shp, x, y[:k//2], m[:k//2])
m[k//2:] = shp_mask(shp, x, y[k//2:], m[k//2:])
else:
m[:k//2, :l//2] = shp_mask(shp, x[:l//2], y[:k//2], m[:k//2, :l//2])
m[:k//2, l//2:] = shp_mask(shp, x[l//2:], y[:k//2], m[:k//2, l//2:])
m[k//2:, :l//2] = shp_mask(shp, x[:l//2], y[k//2:], m[k//2:, :l//2])
m[k//2:, l//2:] = shp_mask(shp, x[l//2:], y[k//2:], m[k//2:, l//2:])
return m
@arunganesan

This comment has been minimized.

Copy link

@arunganesan arunganesan commented May 13, 2019

👍 👍 👍 Super fast code. Thank you

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.