Kiểm tra xem một địa chỉ có vĩ độ và kinh độ nằm trong một shapefile

Question

Làm cách nào để kiểm tra xem một địa chỉ trong vùng của một shapefile đã cho? Tôi maneged để tải một shapefile trong python, nhưng không thể nhận được bất kỳ thêm.

Answer 1

Đây là sự thích nghi của câu trả lời của yosukesabai.

Tôi muốn đảm bảo rằng điểm tôi đang tìm kiếm nằm trong cùng hệ thống chiếu với shapefile, vì vậy tôi đã thêm mã cho điều đó.

Tôi không thể hiểu lý do tại sao anh ấy đang thực hiện thử nghiệm chứa trên ply = feat_in.GetGeometryRef() (trong thử nghiệm của tôi dường như không hoạt động tốt), vì vậy tôi đã xóa nó.

Tôi cũng đã cải thiện nhận xét để giải thích rõ hơn những gì đang diễn ra (như tôi đã hiểu).

#!/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, và this site là hữu ích liên quan đến việc kiểm tra chiếu. EPSG:4326

Answer 2

Một cách để làm điều này là để đọc các tập tin Shape ESRI sử dụng OGR thư viện http://www.gdal.org/ogr và sau đó sử dụng GEOS hình học thư viện http://trac.osgeo.org/geos/ để làm bài kiểm tra point-in-polygon. Điều này yêu cầu một số chương trình C/C++.

Ngoài ra còn có giao diện python cho GEOS tại http://sgillies.net/blog/14/python-geos-module/ (mà tôi chưa từng sử dụng). Có lẽ đó là những gì bạn muốn?

Một giải pháp khác là sử dụng thư viện http://geotools.org/. Đó là trong Java.

Tôi cũng có phần mềm Java của riêng mình để thực hiện việc này (bạn có thể tải xuống từ http://www.mapyrus.org cộng jts.jar từ http://www.vividsolutions.com/products.asp). Bạn chỉ cần một lệnh văn bản tập tin inside.mapyrus chứa những dòng sau để kiểm tra xem một điểm đặt bên trong đa giác đầu tiên trong ESRI Shape file:

dataset "shapefile", "us_states.shp" 
fetch 
print contains(GEOMETRY, -120, 46) 

Và chạy với:

java -cp mapyrus.jar:jts-1.8.jar org.mapyrus.Mapyrus inside.mapyrus 

Nó sẽ in 1 nếu điểm bên trong, 0 nếu không.

Bạn cũng có thể nhận được một số câu trả lời tốt nếu bạn đăng câu hỏi này trên https://gis.stackexchange.com/

Answer 3

Thanh toán http://geospatialpython.com/2011/01/point-in-polygon.html và http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

Answer 4

tôi đã gần như chính xác những gì bạn đang làm hôm qua sử dụng OGR GDAL với python ràng buộc. Nó trông như thế này

import ogr 

# load the shape file as a layer 
drv = ogr.GetDriverByName('ESRI Shapefile') 
ds_in = drv.Open("./shp_reg/satreg_etx12_wgs84.shp") 
lyr_in = ds_in.GetLayer(0) 

# field index for which i want the data extracted 
# ("satreg2" was what i was looking for) 
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("satreg2") 


def check(lon, lat): 
    # create point geometry 
    pt = ogr.Geometry(ogr.wkbPoint) 
    pt.SetPoint_2D(0, lon, lat) 
    lyr_in.SetSpatialFilter(pt) 

    # go over all the polygons in the layer see if one include the point 
    for feat_in in lyr_in: 
    # roughly subsets features, instead of go over everything 
    ply = feat_in.GetGeometryRef() 

    # test 
    if ply.Contains(pt): 
     # TODO do what you need to do here 
     print(lon, lat, feat_in.GetFieldAsString(idx_reg)) 

Answer 5

Một lựa chọn khác là sử dụng kiểu dáng cân đối (một thư viện Python dựa trên GEOS, động cơ cho PostGIS) và Fiona (đó là cơ bản để đọc/ghi tập tin):

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." 

Lưu ý rằng làm như point-in kiểm tra -polygon có thể tốn kém nếu đa giác là lớn/phức tạp (ví dụ, shapefiles cho một số quốc gia có đường bờ biển cực kỳ bất thường). Trong một số trường hợp nó có thể giúp để sử dụng hộp bounding để nhanh chóng loại trừ những điều trên trước khi thực hiện thử nghiệm chuyên sâu hơn:

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

if bounding_box.contains(point): 
    ... 

Cuối cùng, hãy nhớ rằng phải mất một thời gian để tải và phân tích cú pháp/shapefile bất thường lớn (không may, những loại đa giác này thường đắt tiền để giữ trong bộ nhớ, quá).

Answer 6

Nếu bạn muốn tìm hiểu đa giác nào (từ một hình dạng đầy đủ) chứa một điểm nhất định (và bạn cũng có một số điểm), cách nhanh nhất là sử dụng postgis. Tôi thực sự thực hiện một phiên bản dựa trên fiona, bằng cách sử dụng các câu trả lời ở đây, nhưng nó đã rất chậm (tôi đã sử dụng đa xử lý và kiểm tra hộp giới hạn đầu tiên). 400 phút xử lý = 50k điểm. Sử dụng postgis, mất ít hơn 10 giây. Chỉ số cây B hiệu quả!

shp2pgsql -s 4326 shapes.shp > shapes.sql 

Điều đó sẽ tạo ra một file sql với các thông tin từ các shapefile, tạo ra một cơ sở dữ liệu với sự hỗ trợ PostGIS và chạy sql đó. Tạo chỉ mục gist trên cột geom. Sau đó, để tìm tên của đa giác:

sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));" 
cur.execute(sql,(x,y)) 

Answer 7

Dưới đây là một giải pháp đơn giản dựa trên pyshp và shapely.

Giả sử rằng shapefile của bạn chỉ chứa một đa giác (nhưng bạn có thể dễ dàng thích ứng cho nhiều đa giác):

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)