输入/输出矢量shapefile的测试点 [英] Testing point with in/out of a vector shapefile
问题描述
这是我的问题.
- 多边形类型的shapefile代表研究区域
http://i8.tietuku.com/08fdccbb7e11c0a9.png
- 位于整个矩形地图中的某个点
http://i8.tietuku.com/877f87022bf817b8.png
我想测试每个点是否位于多边形内/外,并做一些进一步的操作(例如,求和研究区域内的网格点数量)
I want to test whether each point was located within/out the polygon and do some further operation(for example, sum the grid point amount within the study area)
由于有堆栈溢出的信息,我有两种方法.
I have two methods all thanks to the information on stack overflow.
将shapefile栅格化为栅格文件,然后进行测试.
Rasterize the shapefile into raster file and then test.
I haven't done that yet, but I have asked one question here and get an answer.
我尝试使用poly.contain()
测试散点的位置,但结果与实际情况不符.
I have tried to using poly.contain()
to test the scatter point's location, but the result wasn't match with the reality.
例如:
- 原始数据由 pt (一个熊猫数据框)表示,其中包含1000个网格X,Y. 我已经显示的
- shapefile是研究区域,我想过滤原始数据,只保留该区域内的点.
- Original data represent by pt(a pandas Dataframe) which contain 1000 grids X,Y.
- shapefile I already shown was the study area, I want to filter the original data leaving only point within this area.
# 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])
- 然后,我有一个numpy数组蒙版,值0,1代表内部/外部.
- 将蒙版组合成 pt ==> pt = pt [[pt.mask == 1]],我可以过滤点
- Then, I have a numpy array mask with value 0,1 represent the within/out.
- Combine mask into pt ==> pt = pt[[pt.mask == 1]], I can filter the points
但是问题出在使用
poly,contains_point(xy)
,我无法获得与尝试匹配的结果.
But the problem is using
poly,contains_point(xy)
, I couldn't get the results match with my attempt.
我的想法2的示例
总和为0,1:
An example for my idea 2
sum the 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
从唯一值开始,shapefile区域内必须有3个点,但结果除此以外还显示1个点.
From unique value, there must by 3 point within the shapefile area, but the result shows one point besides.
另一项40分的测试
http://i4.tietuku.com/5fc12514265b5a50.png
结果是错误的,我还没有弄清楚.
但是我认为该问题可能是由两个原因引起的:
The result was wrong, and I haven't figured it out.
But I think the problem may happen by two reasons:
- 多边形shapefile错误(我认为问题不存在的简单多边形)
- 使用
poly.contains_point(xy)
不正确.
- the polygon shapefile was wrong(a simple polygon which I don't think the problem remains here).
- Using
poly.contains_point(xy)
incorrect.
感谢您的回答,我发现的原因是shapefile本身.
当我将其更改为shape.polygon时,它可以很好地工作.
Thanks for the answer, the reason I found out was the shapefile itself.
When I change it into shapely.polygon, it works well.
这是我的代码和结果
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
推荐答案
您可以匀称使用:
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
import pylab as 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加快计算速度:
You can also use MultiPoint to speed the 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))
您还可以在matplotlib中使用Polygon.contains_point()
:
you can also use Polygon.contains_point()
in matplotlib:
poly = pl.Polygon(poly_data)
mask = [poly.contains_point(p) for p in points]
这篇关于输入/输出矢量shapefile的测试点的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!