Python: 如何检查形状区域是否包含点？

Question

我有一个 shapefile 的新加坡，我是这样想象的：

import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
### Plot shapefile
sf = shapefile.Reader("myShape.shp")
recs    = sf.records()
shapes  = sf.shapes()
Nshp    = len(shapes)
cns     = []
for nshp in xrange(Nshp):
    cns.append(recs[nshp][1])
cns = array(cns)
cm    = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
#   -- plot --
fig,ax = plt.subplots(figsize=(20,10))

for nshp in xrange(Nshp):
    ptchs   = []
    pts     = array(shapes[nshp].points)
    prt     = shapes[nshp].parts
    par     = list(prt) + [pts.shape[0]]
    for pij in xrange(len(prt)):
        ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
        ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1))     
ax.set_xlim(103.6,104.1)
ax.set_ylim(1.15,1.48)

现在我有一个坐标在 dataframe 中的点列表，我想检查这些点在 shapefile 的哪个区域。

mypoints
    lon         lat
0   103.619740  1.280485
1   103.622632  1.268944
2   103.622632  1.274714
3   103.622632  1.277600
4   103.622632  1.280485

特别是从 shapefile 我可以提取特定区域的信息，例如 name。事实上，在 recs 我有这个信息。例如：

recs[0]:

[42,
 1,
 'STRAITS VIEW',
 'SVSZ01',
 'Y',
 'STRAITS VIEW',
 'SV',
 'CENTRAL REGION',
 'CR',
 '21451799CA1AB6EF',
 datetime.date(2016, 5, 11),
 30832.9017,
 28194.0843,
 5277.76082729,
 1127297.23737]

在这种情况下，STRAITS VIEW 是区域的名称。最后我想要一个 dataframe 像：

mypoints
    lon         lat       name
0   103.619740  1.280485  STRAITS VIEW
1   103.622632  1.268944  STRAITS VIEW
2   103.622632  1.274714  DOWNTOWN
3   103.622632  1.277600  BEDOK
4   103.622632  1.280485  CHANGI

Answer 1

要检查坐标是否位于某个多边形内，您可以使用 GDAL 中的 Intersects 方法。此方法比较两个 OGRGeometry 对象和 return 指示它们是否拦截的布尔值。

有了这个，您可以遍历您的点，并针对每个点测试它是否 Intersects() 与您的所有多边形区域。它应该看起来像：

for point_geom in your_points:
    #iterate over your area polygons and check 
    for area_geom in your_areas:
        if area_geom.Intersects(point_geom):
            #then you have a match, save or do as you want

注意： 如果您的观点是 而不是 OGRGeometry（您的多边形很可能与您正在阅读的一样来自 Shapefile），您可以按照建议创建此类几何体 here:

from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(1198054.34, 648493.09)

对于多边形：

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
ring.AddPoint(1161053.0218226474, 667456.2684348812)
#etc...add all your points, in a loop would be better
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

Python: 如何检查形状区域是否包含点？

Python: how to check if shape areas contain points?

python

shapefile