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
Copy link

arunganesan commented May 13, 2019

👍 👍 👍 Super fast code. Thank you

@conorssmith
Copy link

conorssmith commented Aug 19, 2021

Thank you for this!

@perrette
Copy link
Author

perrette commented Aug 19, 2021

FYI I dont use this code anymore. Rasterio offers faster code to do the same thing, even though it's a little tricky to use. If you're interested I can search for relevant bits of code.

@CHIRANJIT007-HUE
Copy link

CHIRANJIT007-HUE commented Apr 22, 2022

@perrette would be very great if you could share any rasterio sample code which could use any shape file and can crop the data using that shape file only.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment