如何检查具有纬度和经度的地理点是否在Shapefile中?

内容来源于 Stack Overflow,并遵循CC BY-SA 3.0许可协议进行翻译与使用

  • 回答 (2)
  • 关注 (0)
  • 查看 (110)

如何检查地理位置是否在给定形状文件的区域内?我设法在python中加载一个shapefile,但是不能再继续了。

提问于
用户回答回答于

我想确保我搜索的点与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)
用户回答回答于

另一种选择是使用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."

注意,如果多边形很大/很复杂,做点多边形测试可能会很昂贵(例如,对于一些海岸线极其不规则的国家来说,形状文件)。在某些情况下,在进行更密集的测试之前,使用边界框可以快速排除事情:

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

if bounding_box.contains(point):
    ...

最后,请记住,加载和解析大型/不规则形状文件需要一些时间(不幸的是,这些类型的多边形在内存中也常常很昂贵)。

扫码关注云+社区