Here is my question:
- the 2-d numpy array data represent some property of each grid space
- the shapefile as the administrative division of the study area(like a city).
For example:
http://i4.tietuku.com/84ea2afa5841517a.png
The whole area has 40x40 grids network, and I want to extract the data inside the purple area. In other words , I want to mask the data outside the administrative
boundary into np.nan.
My early attempt
I label the grid number and select the specific array data into np.nan.
http://i4.tietuku.com/523df4783bea00e2.png
value[0,:] = np.nan
value[1,:] = np.nan
.
.
.
.
Can Someone show me a easier method to achieve the target?
Add
Found an answer here which can plot the raster data into shapefile, but the data itself doesn't change.
Update -2016-01-16
I have already solved this problem inspired by some answers.
Someone which are interested on this target, check this two posts which I have asked:
1. Testing point with in/out of a vector shapefile
2. How to use set clipped path for Basemap polygon
The key step was to test the point within/out of the shapefile which I have already transform into shapely.polygon.
Step 1. Rasterize shapefile
Create a function that can determine whether a point at coordinates (x, y)
is or is not in the area. See here for more details on how to rasterize your shapefile into an array of the same dimensions as your target mask
def point_is_in_mask(mask, point):
# this is just pseudocode
return mask.contains(point)
Step 2. Create your mask
mask = np.zeros((height, width))
value = np.zeros((height, width))
for y in range(height):
for x in range(width):
if not point_is_in_mask(mask, (x, y)):
value[y][x] = np.nan
Best is to use matplotlib
:
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)
"""
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
alternatively, you may use shapely
contains method as suggested in the above answer. You may speed-up calculations by recursively sub-dividing the space, as indicated in this gist (but matplotlib solution was 1.5 times faster in my tests):
https://gist.github.com/perrette/a78f99b76aed54b6babf3597e0b331f8