webentwicklung-frage-antwort-db.com.de

Prüfen Sie, ob ein Geopoint mit Längen- und Breitengrad in einem Shapefile enthalten ist

Wie kann ich überprüfen, ob sich ein Geopoint im Bereich eines bestimmten Shapefiles befindet? 

Ich habe es geschafft, ein Shapefile in Python zu laden, komme aber nicht weiter.

20
Gerald Bäck

Dies ist eine Anpassung der Antwort von Yosukesabai.

Ich wollte sicherstellen, dass der Punkt, den ich suchte, im selben Projektionssystem wie das Shapefile war, also habe ich dafür Code hinzugefügt.

Ich konnte nicht verstehen, warum er einen Test mit ply = feat_in.GetGeometryRef() durchführte (in meinem Test schien das alles genauso gut zu funktionieren), also entfernte ich das.

Ich habe auch den Kommentar verbessert, um besser zu erklären, was los ist (so wie ich es verstehe).

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

Diese Seite , diese Seite und diese Seite waren hilfreich bei der Projektionsprüfung. EPSG: 4326

16
Richard

Eine weitere Option ist die Verwendung von Shapely (einer Python-Bibliothek, die auf GEOS, der Engine für PostGIS, basiert) und Fiona (die hauptsächlich zum Lesen/Schreiben von Dateien dient):

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

Beachten Sie, dass Punkt-in-Polygon-Tests teuer sein können, wenn das Polygon groß/kompliziert ist (z. B. Shapefiles für einige Länder mit extrem unregelmäßigen Küstenlinien). In einigen Fällen kann es hilfreich sein, Bounding-Boxen zu verwenden, um Dinge schnell auszuschließen, bevor der intensivere Test durchgeführt wird:

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

if bounding_box.contains(point):
    ...

Denken Sie schließlich daran, dass das Laden und Parsen großer/unregelmäßiger Shapefiles einige Zeit dauert (leider sind diese Arten von Polygonen oft auch teuer im Speicher).

29
Clint Harris

Hier ist eine einfache Lösung basierend auf pyshp und shapely .

Nehmen wir an, dass Ihr Shapefile nur ein Polygon enthält (Sie können es jedoch problemlos für mehrere Polygone anpassen):

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

ich habe fast genau das gemacht, was du gestern mit gdal ogr mit python-bindung gemacht hast. Es sah so aus.

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

Eine Möglichkeit, dies zu tun, besteht darin, die ESRI-Shape-Datei mit der OGR-Bibliothek http://www.gdal.org/ogr zu lesen und dann die GEOS-Geometrie - Bibliothek zu verwenden http://trac.osgeo.org/geos/ um den Punkt-in-Polygon-Test durchzuführen. Dies erfordert einige C/C++ - Programmierung.

Es gibt auch eine Python-Schnittstelle zu GEOS unter http://sgillies.net/blog/14/python-geos-module/ (die ich noch nie verwendet habe). Vielleicht willst du das?

Eine andere Lösung ist die Verwendung der Bibliothek http://geotools.org/ . . Dies ist in Java.

Ich habe dazu auch meine eigene Java-Software (die Sie Von http://www.mapyrus.org plus jts.jar von http: //www.vividsolutions herunterladen können. de/products.asp ). Sie benötigen nur einen Textbefehl Datei inside.mapyrus, der Die folgenden Zeilen enthält, um zu prüfen, ob ein Punkt innerhalb des ersten Polygons Der ESRI-Shape-Datei liegt:

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

Und laufen mit:

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

Es wird eine 1 gedruckt, wenn der Punkt innen ist, andernfalls 0.

Wenn Sie diese Frage unter https://gis.stackexchange.com/ posten, erhalten Sie möglicherweise einige gute Antworten.

1
Simon C

Wenn Sie herausfinden möchten, welches Polygon (aus einem Shapefile voll davon) einen bestimmten Punkt enthält (und Sie auch eine Reihe von Punkten haben), ist der schnellste Weg die Verwendung von Postgis. Ich habe tatsächlich eine auf Fiona basierende Version implementiert, die Antworten hier verwendet, aber es war schmerzhaft langsam (ich verwendete zuerst Multiprocessing und überprüfte zuerst den Begrenzungsrahmen). 400 Minuten Verarbeitungszeit = 50.000 Punkte. Das Verwenden von Postgis dauerte weniger als 10 Sekunden. B-Baum-Indizes sind effizient!

shp2pgsql -s 4326 shapes.shp > shapes.sql

Dadurch wird eine SQL-Datei mit den Informationen aus den Shapefiles generiert, eine Datenbank mit Postgis-Unterstützung erstellt und diese SQL ausgeführt. Erstellen Sie einen Gist-Index für die Spalte "Geom". Dann, um den Namen des Polygons zu finden:

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