Ich habe mehrere Codebeispiele mit Bibliotheken wie Shapefile, Fiona und Ogr ausprobiert, um zu überprüfen, ob ein Punkt (x, y) innerhalb der Grenzen eines mit ArcMap erstellten Multipolygons (und damit im Shapefile-Format) liegt. Keines der Beispiele funktioniert jedoch gut mit Multipolygonen, obwohl sie mit regulären Einzelpolygon-Shapefiles gut funktionieren. Einige Ausschnitte, die ich ausprobiert habe, sind unten aufgeführt:
# First example using shapefile and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile
polygon = shapefile.Reader('shapefile.shp')
polygon = polygon.shapes()
shpfilePoints = []
for shape in polygon:
shpfilePoints = shape.points
polygon = shpfilePoints
poly = Polygon(poly)
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
# Second example using ogr and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
from osgeo import ogr, gdal
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open("shapefile.shp", 0)
layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
feature = layer.GetFeature(index)
geometry = feature.GetGeometryRef()
polygon = Polygon(geometry)
print 'polygon points =', polygon # this prints 'multipoint' + all the points fine
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
Das erste Beispiel funktioniert gut mit jeweils einem einzelnen Polygon. Wenn ich jedoch einen Punkt innerhalb einer der Formen in meinem Multipolygon-Shapefile eingebe, wird "out" zurückgegeben, obwohl er in einen der vielen Teile fällt.
Für das zweite Beispiel erhalte ich den Fehler "Objekt vom Typ 'Geometrie' hat kein len ()", von dem ich annehme, dass das Geometriefeld nicht als normale, indizierte Liste / Array gelesen werden kann.
Ich habe außerdem versucht, den tatsächlichen Punkt im Polygoncode wie hier vorgeschlagen zu ersetzen , um sicherzustellen, dass nicht dieser Teil des Codes nicht funktioniert hat. Und während die Beispiele dieses Links mit einfachen Polygon-Shapefiles gut funktionieren, kann ich mein komplexes Multipolygon nicht richtig testen.
Ich kann mir also keine andere Möglichkeit vorstellen, um zu testen, ob ein Punkt über Python in ein Multipolygon-Shapefile fällt ... Vielleicht fehlen mir noch andere Bibliotheken?
quelle
polygon = Polygon(geometry)
eine Art Try-Schleife ersetzen, in die umgeschaltet wird,polygon = MultiPolygon(geometry)
wenn dieser Fehler auftritt.Antworten:
Shapefiles haben keinen Typ MultiPolygon (Typ = Polygon), aber sie unterstützen sie trotzdem (alle Ringe werden in einem Feature gespeichert = Liste der Polygone, siehe Konvertieren eines großen Multipolygons in Polygone )
Das Problem
Wenn ich ein MultiPolygon-Shapefile öffne, lautet die Geometrie "Polygon".
Lösung 1 mit Fiona
Fiona interpretiert die Funktion als MultiPolygon und Sie können die in Python für effizientere räumliche Verknüpfung vorgestellte Lösung ohne QGIS, ArcGIS, PostGIS usw. anwenden. (1)
Lösung 2 mit pyshp (Shapefile) und dem Protokoll geo_interface (GeoJSON like)
Dies ist eine Ergänzung zur Antwort von xulnik.
Lösung 3 mit ogr und dem geo_interface- Protokoll ( Python Geo_interface-Anwendungen )
Lösung 4 mit GeoPandas wie in Efficient Spatial Join in Python ohne QGIS, ArcGIS, PostGIS usw. (2)
Die Punkte 1,3,5,6 liegen innerhalb der Grenzen des MultiPolygons
quelle
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)
Lösung 2 auf? Ich kann keine Aufrufe einer shape () -Methode von erhaltenshapefile.py
. Ich habe es sogar versuchtshapefile.Shape()
; Es gibt eine Klasse dafür, aber es funktioniert nicht.within()
Methode?from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
)File "C:\WinPython\python-3.6.5.amd64\lib\site-packages\geopandas\tools\sjoin.py", line 43, in sjoin if left_df.crs != right_df.crs:
AttributeError: 'MultiPolygon' object has no attribute 'crs'
Das Problem in Ihrem ersten Beispiel liegt in dieser Schleife:
Es werden nur die letzten Feature-Punkte angehängt. Ich habe meinen Ansatz mit diesem Shapefile ausprobiert:
Ich habe Ihren Code geändert in:
Der obige Code wurde in der Python-Konsole von QGIS ausgeführt und das Ergebnis war:
Es funktioniert perfekt und jetzt können Sie überprüfen, ob ein Punkt (x, y) innerhalb der Grenzen jedes Features liegt.
quelle