如何检查地理位置是否在给定shapefile的区域内?
我设法用python加载了一个shapefile,但是不能再继续了。
发布于 2012-11-18 01:49:28
这是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)
发布于 2013-09-12 03:11:50
另一种选择是使用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需要一些时间(不幸的是,这些类型的多边形在内存中保存的代价也很高)。
发布于 2017-03-10 02:08:55
让我们假设您的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)
https://stackoverflow.com/questions/7861196
复制相似问题