问题描述
如何检查地理点是否在给定 shapefile 的区域内?
How can I check if a geopoint is within the area of a given shapefile?
我设法在 python 中加载了一个 shapefile,但无法继续.
I managed to load a shapefile in python, but can't get any further.
推荐答案
这是对yosukesabai的回答的改编.
This is an adaptation of yosukesabai's answer.
我想确保我正在搜索的点与 shapefile 在同一个投影系统中,因此我为此添加了代码.
I wanted to ensure that the point I was searching for was in the same projection system as the shapefile, so I've added code for that.
我不明白他为什么要对 ply = feat_in.GetGeometryRef() 进行包含测试(在我的测试中,没有它似乎也能正常工作),所以我删除了它.
I couldn't understand why he was doing a contains test on ply = feat_in.GetGeometryRef() (in my testing things seemed to work just as well without it), so I removed that.
我还改进了评论以更好地解释正在发生的事情(据我了解).
I've also improved the commenting to better explain what's going on (as I understand it).
#!/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
This site, this site, and this site were helpful regarding the projection check. EPSG:4326