矢量 shapefile 的 in/out 测试点

Question

这是我的问题。

1。简介

多边形类型的 shapefile 表示研究区域

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

某个点位于整个矩形地图中

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

我想测试每个点是否位于within/out多边形上，并做一些进一步的操作（例如，求和研究区域内的网格点数量）

2。我的想法

多亏了堆栈溢出的信息，我有两个方法。

2.1 想法A

将shapefile栅格化为栅格文件，然后进行测试。

我还没有这样做，但我问了一个问题并得到了答案。

2.2 想法B

我试过用poly.contain()测试散点的位置，结果与实际不符。

3。我的代码基于想法 B：

例如：

原始数据由pt（一个pandas Dataframe）表示，其中包含1000个网格X，Y。
shapefile我已经显示的是研究区域，我想过滤原始数据只留下这个区域内的点。

3.1 准备

# map four 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 准备

# 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)
for i in range(0,len(xx),1):
    data[i] = np.array([xx[i],yy[i]])

# reading the shapefile
              
map = Basemap(llcrnrlon=x_map1,llcrnrlat=y_map1,urcrnrlon=x_map2,\
              urcrnrlat=y_map2)
map.readshapefile('/xx,'xx')

3.2 测试

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

然后，我有一个数值为 0,1 的 numpy 数组掩码，代表 within/out.
将mask合并成pt ==> pt = pt[[pt.mask == 1]],我可以过滤点

But the problem is using poly,contains_point(xy), I couldn't get the results match with my attempt.

我的想法 2 的示例

求和值 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

根据唯一值，shapefile 区域内必须有 3 个点，但结果显示除此之外还有一个点。

再考40分

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

4。我的问题

结果不对，还没想好
但我认为这个问题可能有两个原因：

多边形 shapefile 错误（我认为问题不存在的简单多边形）。
使用 poly.contains_point(xy) 不正确。

添加 2016-01-16

感谢您的回答，我发现的原因是 shapefile 本身。
当我把它改成shapely.polygon时，效果很好。

这是我的代码和结果

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)
for i in range(0,len(xx),1):
    points[i] = np.array([xx[i],yy[i]])
mask = np.array([poly.contains(Point(x, y)) for 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

Answer 1

你可以使用 shapely:

import numpy as np
from 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)) for x, y in points])

这里是剧情代码：

将 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")

输出：

您还可以使用 MultiPoint 来加速计算：

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))

你也可以在matplotlib中使用Polygon.contains_point():

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

矢量 shapefile 的 in/out 测试点

Testing point with in/out of a vector shapefile

python

numpy

matplotlib

shapefile

matplotlib-basemap

1。简介

2。我的想法

3。我的代码基于想法 B：

4。我的问题

添加 2016-01-16