본문 바로가기
파이썬

파이썬 위도와 경도가있는 지리 점이 shapefile 내에 있는지 확인

by º기록 2020. 9. 27.
반응형

지오 포인트가 주어진 shapefile 영역 내에 있는지 어떻게 확인할 수 있습니까?

파이썬에서 shapefile을로드 할 수 있었지만 더 이상 얻을 수 없습니다.

 

해결 방법

 

이것은 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)


 

참조 페이지 https://stackoverflow.com/questions/7861196

 

 

반응형

댓글