python - Testing point with in/out of a vector shapefile -

here question.

1. intro

• a shapefile in polygon type represent study area

http://i8.tietuku.com/08fdccbb7e11c0a9.png

• some point located in whole rectangle map

http://i8.tietuku.com/877f87022bf817b8.png

i want test whether each point located within/out polygon , further operation(for example, sum grid point amount within study area)

2. idea

i have 2 methods information on stack overflow.

2.1 idea

rasterize shapefile raster file , test.

i haven't done yet, have asked 1 question here , answer.

2.2 idea b

i have tried using poly.contain() test scatter point's location, result wasn't match reality.

3. code based on idea b:

for example:

• original data represent pt(a pandas dataframe) contain 1000 grids x,y.
• shapefile shown study area, want filter original data leaving point within area.

3.1 preparation

# map 4 boundaries xc1,xc2,yc1,yc2 = 113.49805889531724,115.5030664238035,37.39995194888143,38.789235929357105 # grid definition lon_grid  = np.linspace(x_map1,x_map2,38) lat_grid  = np.linspace(y_map1,y_map2,32) 

3.1 preparation

# generate (lon,lat)    xx = lon_grid[pt.x.iloc[:].as_matrix()] yy = lat_grid[pt.y.iloc[:].as_matrix()]  sh = (len(xx),2) data = np.zeros(len(xx)*2).reshape(*sh) in range(0,len(xx),1):     data[i] = np.array([xx[i],yy[i]])  # reading shapefile  map = basemap(llcrnrlon=x_map1,llcrnrlat=y_map1,urcrnrlon=x_map2,\               urcrnrlat=y_map2) map.readshapefile('/xx,'xx') 

3.2 test

patches=[] info, shape in zip(map.xxx_info, map.xxx):     x,y=zip(*shape)     patches.append(polygon(np.array(shape), true) ) poly in patches:      mask = np.array([poly.contains_point(xy) xy in data]) 

• then, have numpy array mask value 0,1 represent within/out.
• combine mask pt ==> pt = pt[[pt.mask == 1]], can filter points

but problem using poly,contains_point(xy), couldn't results match attempt.

an example idea 2

sum value 0,1:

unique, counts = np.unique(mask, return_counts=true)       print np.asarray((unique, counts)).t #result:   > [[0 7]     [1 3]] 

http://i4.tietuku.com/7d156db62c564a30.png

from unique value, there must 3 point within shapefile area, result shows 1 point besides.

another test 40 points

http://i4.tietuku.com/5fc12514265b5a50.png

4. problem

the result wrong, , haven't figured out.
think problem may happen 2 reasons:

• the polygon shapefile wrong(a simple polygon don't think problem remains here).
• using poly.contains_point(xy) incorrect.

add 2016-01-16

thanks answer, reason found out shapefile itself.
when change shapely.polygon, works well.

here code , result

c =    fiona.open("xxx.shp") pol = c.next() geom = shape(pol['geometry']) poly_data = pol["geometry"]["coordinates"][0] poly = polygon(poly_data) ax.add_patch(plt.polygon(poly_data))  xx = lon_grid[pt_select.x.iloc[:].as_matrix()] yy = lat_grid[pt_select.y.iloc[:].as_matrix()]  sh = (len(xx),2) points = np.zeros(len(xx)*2).reshape(*sh) in range(0,len(xx),1):     points[i] = np.array([xx[i],yy[i]]) mask = np.array([poly.contains(point(x, y)) x, y in points])  ax.plot(points[:, 0], points[:, 1], "rx") ax.plot(points[mask, 0], points[mask, 1], "ro")     

http://i4.tietuku.com/8d895efd3d9d29ff.png

you can use shapely:

import numpy np shapely.geometry import polygon, point  poly_data = [[0, 0], [0, 1], [1, 0], [0.2, 0.5]] poly = polygon(poly_data)  points = np.random.rand(100, 2)  mask = np.array([poly.contains(point(x, y)) x, y in points]) 

and here plot code:

import pylab pl

fig, ax = pl.subplots() ax.add_patch(pl.polygon(poly_data)) ax.plot(points[:, 0], points[:, 1], "rx") ax.plot(points[mask, 0], points[mask, 1], "ro") 

the output:

you can use multipoint speed calculation:

from shapely.geometry import polygon, multipoint  poly_data = [[0, 0], [0, 1], [1, 0], [0.2, 0.5]] poly = polygon(poly_data) points = np.random.rand(100, 2) inside_points = np.array(multipoint(points).intersection(poly)) 

you can use polygon.contains_point() in matplotlib:

poly = pl.polygon(poly_data) mask = [poly.contains_point(p) p in points]