检查经纬度的地理点是否在shapefile中 [英] Check if a geopoint with latitude and longitude is within a shapefile

问题描述

如何检查地址是否位于给定shapefile的区域内？

我在python中加载了一个shapefile文件，但无法进一步处理。

解决方案

我想确保我搜索的点与shapefile在同一个投影系统中，所以我已经添加了代码。

我无法理解他为什么要对 ply = feat_in.GetGeometryRef（ ）（在我的测试中，如果没有它，工作似乎也能正常工作），所以我删除了它。

我也改进了评论以更好地解释发生了什么（据我了解）。

 ＃！/ usr / bin / python 
 import ogr 
从IPython导入嵌入
导入sys 
 
 drv = ogr.GetDriverByName（'ESRI Shapefile'）＃我们将加载一个形状文件
 ds_in = drv.Open （MN.shp）＃获取形状文件的内容
 lyr_in = ds_in.GetLayer（0）＃获取形状文件的第一层
 
＃将您感兴趣的字段标题放在这里
 idx_reg = lyr_in.GetLayerDefn（）。GetFieldIndex（P_Loc_Nm）
 
＃如果纬度/经度我们将要使用的不是在shapefile的投影
＃中，那么我们将得到错误的结果。 
＃以下假设纬度经度为WGS84 
＃这是由数字4326标识的，如EPSG：4326
＃我们将在此和shapefile的
 #project，无论它可能是
 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）：
＃将传入经度/纬度转换为shapefile的投影
 [lon， lat，z] = ctran.TransformPoint（lon，lat）
 
＃创建一个点
 pt = ogr.Geometry（ogr.wkbPoint）
 pt.SetPoint_2D（0，lon ，lat）
 
＃设置一个空间过滤器，使我们在
 #loop到lyr_in时看到的唯一特征是与
 lyr_in上方定义的点重叠的特征。 SetSpatialFilter（pt）
 
＃通过重叠的要素循环显示feat_in i中感兴趣的字段
 n lyr_in：
 print lon，lat，feat_in.GetFieldAsString（idx_reg）
 
＃进行命令行输入并执行所有操作
 check（float（sys.argv [1] ），float（sys.argv [2]））
 #check（-95,47）
  

此网站，本站，本网站对投影检查有帮助。 EPSG：4326

How can I check if a geopoint is within the area of a given shapefile?

I managed to load a shapefile in python, but can't get any further.

解决方案

This is an adaptation of yosukesabai's answer.

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.

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)

This site, this site, and this site were helpful regarding the projection check. EPSG:4326

这篇关于检查经纬度的地理点是否在shapefile中的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！