So fügen Sie in GDAL Raster unterschiedlicher Größe hinzu, sodass das Ergebnis nur im Schnittbereich angezeigt wird

11

Ich schreibe eine Python-Methode, die zwei Raster hinzufügt und eine einzelne Rasterausgabe generiert. Aus Gründen, die außerhalb meiner Kontrolle liegen, sind die Ausmaße der Eingabe-Raster unterschiedlich, sie überschneiden sich jedoch.

Ist es möglich, ausschließlich die Eingabe-Rasterpixel zu bearbeiten, die in den 2 überlappenden Bereichen überlappt sind, um meine Ausgabe so zu generieren, dass die Ausgabe-Raster-Ausdehnung nur der Schnittbereich der beiden Eingaben ist?

python raster gdal Reich
quelle

Antworten:

12

Als erstes müssen Sie das überlappende Rechteck in räumlichen Koordinaten bestimmen. Dazu erhalten Sie die Geotransform für jedes Quellbild:

gt1 = ds1.GetGeoTransform()
# r1 has left, top, right, bottom of dataset's bounds in geospatial coordinates.
r1 = [gt1[0], gt1[3], gt1[0] + (gt1[1] * ds1.RasterXSize), gt1[3] + (gt1[5] * ds1.RasterYSize)]

# Do the same for dataset 2 ...

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Konvertieren Sie dann dieses Rechteck für jedes Bild in Pixel, indem Sie die oberen und linken Koordinaten subtrahieren und durch die Pixelgröße dividieren und aufrunden.

Von hier aus können Sie ReadRaster()jedes Bild aufrufen und ihm die soeben berechneten Pixelbereiche geben:

band.ReadRaster(px1[0], px1[1], px1[2] - px1[0], px1[3] - px1[1], px1[2] - px1[0], px1[3] - px1[1],
                # <band's datatype here>
)

Ich bin ein bisschen müde. Wenn das nicht viel Sinn macht, lass es mich wissen!

MerseyViking
quelle
Funktioniert dies auch für Raster mit unterschiedlichen Projektionen (auch Koordinatenreferenzsysteme / Raumbezugssysteme genannt)? Und selbst wenn die Projektionen gleich sind: Funktioniert das auch, wenn gt1[1]und gt2[1](oder gt1[5]und gt2[5]) entgegengesetzte Vorzeichen haben? (Welches würde eines der Raster vertikal oder horizontal drehen, denke ich.) Oder ob abs(gt1[2])und abs(gt1[4])größer als abs(gt1[1])und abs(gt1[5])aber abs(gt2[2])und abs(gt2[4])kleiner als abs(gt2[1])und sind abs(gt2[5])(was wahrscheinlich eines der Raster diagonal drehen würde)?
das-g
6

Das dritte Schnittelement sollte min sein (r1 [2], r2 [2]):

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Außerdem würde ich eine Logik empfehlen, um zu überprüfen, ob sich die Datensätze tatsächlich überschneiden.

Dies ist ein Ansatz:

if (intersection[2] < intersection[0]) or (intersection[1] < intersection[3]):
    intersection = None
David Shean
quelle