检查带纬度和经度的地理点是否在形状文件内

Question

我该如何检查一个地理点是否在某个给定的形状文件区域内呢？

我已经在Python中加载了一个形状文件，但不知道接下来该怎么做。

Answer 1

这里有一个简单的解决方案，使用了 pyshp 和 shapely 这两个库。

假设你的 shapefile 里只包含一个多边形（不过你可以很容易地调整代码来处理多个多边形）：

import shapefile
from shapely.geometry import shape, Point

# read your shapefile
r = shapefile.Reader("your_shapefile.shp")

# get the shapes
shapes = r.shapes()

# build a shapely polygon from your shape
polygon = shape(shapes[0])    

def check(lon, lat):
    # build a shapely point from your geopoint
    point = Point(lon, lat)

    # the contains function does exactly what you want
    return polygon.contains(point)

Answer 2

另一个选择是使用Shapely（这是一个基于GEOS的Python库，GEOS是PostGIS的引擎）和Fiona（主要用于读取和写入文件）：

import fiona
import shapely

with fiona.open("path/to/shapefile.shp") as fiona_collection:

    # In this case, we'll assume the shapefile only has one record/layer (e.g., the shapefile
    # is just for the borders of a single country, etc.).
    shapefile_record = fiona_collection.next()

    # Use Shapely to create the polygon
    shape = shapely.geometry.asShape( shapefile_record['geometry'] )

    point = shapely.geometry.Point(32.398516, -39.754028) # longitude, latitude

    # Alternative: if point.within(shape)
    if shape.contains(point):
        print "Found shape for point."

需要注意的是，如果多边形很大或者很复杂（比如某些国家的海岸线非常不规则的形状文件），进行点在多边形内的测试可能会很耗费资源。在某些情况下，使用边界框可以快速排除一些情况，然后再进行更深入的测试，这样会更有效率：

minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)

if bounding_box.contains(point):
    ...

最后，请记住，加载和解析大型或不规则的形状文件需要一些时间（不幸的是，这些类型的多边形在内存中占用的空间也通常很大）。

Answer 3

这是对yosukesabai回答的一个改编。

我想确保我查找的点和这个shapefile使用的是同一个投影系统，所以我加了一些代码来实现这一点。

我不太明白他为什么要在ply = feat_in.GetGeometryRef()上进行包含测试（在我的测试中，去掉这个测试似乎也能正常工作），所以我把它删掉了。

我还改进了注释，以便更好地解释发生了什么（根据我的理解）。

#!/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)

关于投影检查，这个网站、这个网站和这个网站都提供了很大的帮助。EPSG:4326

检查带纬度和经度的地理点是否在形状文件内

7 个回答

撰写回答