4

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.

Community
  • 1
  • 1
Han Zhengzu
  • 3,694
  • 7
  • 44
  • 94

2 Answers2

3

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
mirosval
  • 6,671
  • 3
  • 32
  • 46
3

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

Mahé
  • 445
  • 4
  • 9