检查具有纬度和经度的地理位置是否在shapefile中

Question

如何检查地理点是否在给定shapefile的区域内？

我设法在python中加载一个shapefile,但无法进一步.

Answer 1

另一个选择是使用Shapely(基于GEOS的Python库,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."

请注意,如果多边形较大/较复杂(例如,某些具有极不规则海岸线的国家/地区的shapefile),则执行多边形点测试可能会很昂贵.在某些情况下,在进行更密集的测试之前,它可以帮助使用边界框来快速排除问题:

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

if bounding_box.contains(point):
    ...

最后,请记住,加载和解析大型/不规则的shapefile需要一些时间(不幸的是,这些类型的多边形通常也很难保存在内存中).

Answer 2

这是对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

Answer 3

这是一个基于pyshp和匀称的简单解决方案.

假设你的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)