Feststellen, ob sich Shapefile und Raster in Python mithilfe von OGR / GDAL überlappen? [geschlossen]

9

Ich erstelle ein Skript in Python mit OGR / GDAL.

Ich habe eine Reihe von Shapefiles und eine Reihe von GeoTiff-Rasterdateien.

Ich möchte, dass mein Skript Shapefiles ignoriert, wenn sie sich nicht mit dem Rasterbereich überschneiden.

Das Shapefile ist kein Rechteck, daher kann ich die von layer.GetExtent () zurückgegebenen xmin / xmax-, ymin / ymax-Werte nicht einfach vergleichen. Ich brauche das eigentliche Polygon, das seine Gesamtform darstellt, und dann eine Möglichkeit, um festzustellen, ob sich dieses Polygon mit dem Rasterquadrat schneidet.

Ich dachte, ich könnte irgendwie alle Polygone im Shapefile zu einem Feature zusammenführen, dann die Geometrie dieses Features lesen und diese Informationen dann mit dem Rasterumfang vergleichen. Ich bin mir jedoch nicht sicher, wie ich dies ausführen soll.

  1. Wie extrahiere ich Randpolygoninformationen aus dem Shapefile?
  2. Wie kann festgestellt werden, ob dieses Polygon eine bestimmte quadratische Fläche schneidet?
python raster shapefile gdal ogr JFerg
quelle
Ich bin mit osgeo nicht vertraut, aber das bogenförmige Äquivalent würde (könnte) beinhalten: Rasterausdehnung lesen, Polygon-Ausdehnungsausdehnung im Speicher erstellen, Shapefiles durchlaufen, jedes auf Ausdehnungsrechteck ausschneiden, testen, ob sich etwas ergibt.
Phloem

Antworten:

17

Das folgende Skript bestimmt den Begrenzungsrahmen eines Rasters und erstellt basierend auf dem Begrenzungsrahmen eine Geometrie. 

import ogr, gdal

raster = gdal.Open('sample.tif')
vector = ogr.Open('sample.shp')

# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

Als nächstes wird die Geometrie des Vektorpolygons bestimmt. Dies beantwortet Ihre erste Frage.

# Get vector geometry
layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

Zuletzt wird die Geometrie des Vektors und des Rasters auf Schnittmenge (Rückgabe Trueoder False) getestet . Dies beantwortet Ihre zweite Frage. 

print rasterGeometry.Intersect(vectorGeometry)
ustroetz
quelle
2
Danke, genau das habe ich gesucht. Diese Antwort zeigt deutlich, wie Funktionen zwischen Geometrieobjekten erstellt, extrahiert und ausgeführt werden. Genau das habe ich gesucht.
JFerg
Diese Lösung testet, ob sich das Polygon FID = 0 mit dem Raster schneidet. Wie erhält man die Geometrie der Shapefile-Summe, wenn kein Polygon dies darstellt?
JFerg
1
BEARBEITEN: Die Verlängerung der Rechenzeit spielt keine Rolle, daher überprüfe ich, ob sich jedes Polygon im Shapefile jetzt schneidet.
JFerg
4

Ich finde die @ustroetz-Lösung sehr hilfreich, aber sie musste an zwei Stellen korrigiert werden. Erstens ist pixelHeight = transform [5] bereits ein negativer Wert, daher sollte die Gleichung lauten

yBottom = yTop+rows*pixelHeight

Zweitens muss die Reihenfolge der Punkte im Ring gegen den Uhrzeigersinn sein. Ich hatte Probleme damit. Richtige Reihenfolge ist:

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xLeft, yTop)
Pandza
quelle