使用shapefile在Python中查找点的封闭多边形[重复]答案

【问题标题】：Find enclosing ploygon for point in Python, using shapefile [duplicate]使用shapefile在Python中查找点的封闭多边形[重复]
【发布时间】：2017-09-26 08:24:11
【问题描述】：

我有一个带有一系列代表县的多边形的 shapefile。 我想知道 Python 中的一个函数/模块，它可以在给定单个点的坐标的情况下推导出“边界”多边形。

换句话说，一个函数使用一个点的纬度/经度坐标和 shapefile，并返回这个 shapefile 中包含该点的多边形。见图表：

如您所见，该点位于“蓝色”多边形内，这是我对任何给定点所需要的。

我知道可能没有内置函数可以执行此操作，但任何有关如何执行此操作的建议都非常好，谢谢！

【问题讨论】：

我认为另一个问题虽然相似，但正在查看一个点是否在 shapefile（所有多边形）的范围内，而不是找出一个点所在的特定多边形/对象 - 我但是会查看这些答案，看看它们是否可以适应。

标签： python polygon point shapefile

【解决方案1】：

理查德对类似问题的回答是最好的。

在这里找到：https://stackoverflow.com/a/13433127/7019148

复制如下以供将来参考：

#!/usr/bin/python
import ogr
from IPython import embed
import sys

drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file
ds_in = drv.Open("MN.shp")    #Get the contents of the shape file
lyr_in = ds_in.GetLayer(0)    #Get the shape file's first layer

#Put the title of the field you are interested in here
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm")

#If the latitude/longitude we're going to use is not in the projection
#of the shapefile, then we will get erroneous results.
#The following assumes that the latitude longitude is in WGS84
#This is identified by the number "4326", as in "EPSG:4326"
#We will create a transformation between this and the shapefile's
#project, whatever it may be
geo_ref = lyr_in.GetSpatialRef()
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)
ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)

def check(lon, lat):
    #Transform incoming longitude/latitude to the shapefile's projection
    [lon,lat,z]=ctran.TransformPoint(lon,lat)

    #Create a point
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    #Set up a spatial filter such that the only features we see when we
    #loop through "lyr_in" are those which overlap the point defined above
    lyr_in.SetSpatialFilter(pt)

    #Loop through the overlapped features and display the field of interest
    for feat_in in lyr_in:
        print lon, lat, feat_in.GetFieldAsString(idx_reg)

#Take command-line input and do all this
check(float(sys.argv[1]),float(sys.argv[2]))
#check(-95,47)

【讨论】：