it-swarm-id.com

Periksa apakah geopoint dengan lintang dan bujur ada di dalam shapefile

Bagaimana saya bisa mengecek apakah geopoint berada di dalam area shapefile tertentu? 

Saya berhasil memuat shapefile dengan python, tetapi tidak bisa melanjutkan.

20
Gerald Bäck

Ini adalah adaptasi dari jawaban yosukesabai.

Saya ingin memastikan bahwa titik yang saya cari berada di sistem proyeksi yang sama dengan shapefile, jadi saya telah menambahkan kode untuk itu.

Saya tidak bisa mengerti mengapa dia melakukan pengujian isi pada ply = feat_in.GetGeometryRef() (dalam pengujian saya, hal-hal tampaknya bekerja dengan baik tanpa itu), jadi saya menghapusnya.

Saya juga meningkatkan komentar untuk lebih menjelaskan apa yang terjadi (seperti yang saya mengerti).

#!/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)

Situs ini , situs ini , dan situs ini sangat membantu mengenai pemeriksaan proyeksi. EPSG: 4326

16
Richard

Pilihan lain adalah dengan menggunakan Shapely (pustaka Python berdasarkan GEOS, mesin untuk PostGIS) dan Fiona (yang pada dasarnya untuk membaca/menulis file):

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

Perhatikan bahwa melakukan pengujian titik-dalam-poligon bisa mahal jika poligon itu besar/rumit (mis., Shapefile untuk beberapa negara dengan garis pantai yang sangat tidak teratur). Dalam beberapa kasus dapat membantu menggunakan kotak pembatas untuk dengan cepat mengesampingkan hal-hal sebelum melakukan tes yang lebih intensif:

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

if bounding_box.contains(point):
    ...

Terakhir, perlu diingat bahwa perlu waktu untuk memuat dan menguraikan shapefile besar/tidak beraturan (sayangnya, jenis-jenis poligon itu seringkali mahal untuk disimpan dalam memori juga).

29
Clint Harris

Berikut adalah solusi sederhana berdasarkan pyshp dan rupawan .

Mari kita asumsikan bahwa shapefile Anda hanya berisi satu poligon (tetapi Anda dapat dengan mudah beradaptasi dengan beberapa poligon):

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)
7
chilladx
2
ViennaMike

saya melakukan hampir persis apa yang Anda lakukan kemarin menggunakan gdal ogr dengan python binding Terlihat seperti ini.

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))
2
yosukesabai

Salah satu cara untuk melakukan ini adalah membaca file ESRI Shape menggunakan OGR Library http://www.gdal.org/ogr dan kemudian gunakan GEOS geometry Library http://trac.osgeo.org/geos/ untuk melakukan tes point-in-polygon. Ini membutuhkan beberapa pemrograman C/C++.

Ada juga antarmuka python ke GEOS di http://sgillies.net/blog/14/python-geos-module/ (yang belum pernah saya gunakan). Mungkin itu yang Anda inginkan?

Solusi lain adalah dengan menggunakan http://geotools.org/ library. Itu di Jawa.

Saya juga punya perangkat lunak Java sendiri untuk melakukan ini (yang dapat Anda unduh Dari http://www.mapyrus.org plus jts.jar dari http://www.vividsolutions.com/products .asp ). Anda hanya perlu perintah teks File inside.mapyrus yang berisi Baris berikut untuk memeriksa apakah suatu titik terletak di dalam poligon pertama Dalam file ESRI Shape:

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

Dan jalankan dengan:

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

Ini akan mencetak 1 jika intinya ada di dalam, 0 sebaliknya.

Anda mungkin juga mendapatkan beberapa jawaban yang bagus jika Anda memposting pertanyaan ini di https://gis.stackexchange.com/

1
Simon C

Jika Anda ingin mengetahui poligon mana (dari bentuk penuh mereka) yang mengandung titik tertentu (dan Anda juga memiliki banyak titik), cara tercepat menggunakan postgis. Saya benar-benar mengimplementasikan versi berbasis fiona, menggunakan jawaban di sini, tetapi sangat lambat (saya menggunakan multiprocessing dan memeriksa kotak pembatas terlebih dahulu). 400 menit pemrosesan = 50k poin. Menggunakan postgis, itu membutuhkan waktu kurang dari 10 detik. Indeks B tree efisien!

shp2pgsql -s 4326 shapes.shp > shapes.sql

Itu akan menghasilkan file sql dengan informasi dari shapefile, membuat database dengan dukungan postgis dan menjalankan sql itu. Buat indeks Gist pada kolom geom. Kemudian, untuk menemukan nama poligon:

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