Georeferenzierungsraster mit GDAL und Python?

10

Ich möchte ein Raster mit pythonund georeferenzieren GDAL. Mein aktueller Ansatz ist es, eine hässliche Liste von Bodenkontrollpunkten aufzurufen gdal_translateund zu gdalwarpverwenden os.system. Ich würde wirklich gerne einen Weg finden, dies von Haus aus zu tun python.

Dies ist der aktuelle Prozess, den ich verwende:

import os

os.system('gdal_translate -of GTiff -gcp 1251.92 414.538 -7.9164e+06 5.21094e+06 -gcp 865.827 107.699 -7.91651e+06 5.21104e+06  "inraster.tif" "outraster1.tif"')
os.system('gdalwarp -r bilinear -tps -co COMPRESS=NONE "outraster2.tif" "outraster3.tif"')

Es gibt eine frühere Frage und Antwort aus dem Jahr 2012, auf die gdal_translatenach dem Import zugegriffen werden kann gdal. Ich bin nicht sicher, ob es veraltet ist oder ob es falsch ist, aber wenn ich es ausführe, from osgeo import gdalsehe ich es nicht gdal.gdal_translateals Option.

Ich weiß nicht, ob es existiert, aber ich würde es lieben, wenn ich Raster auf pythonische Weise übersetzen und neu projizieren könnte. Zum Beispiel:

# translate
gcp_points = [(1251.92, 414.538), (-7.9164e+06, 5.21094e+06)]
gdal.gdal_translate(in_raster, gcp_points, out_raster1)

# warp
gdal.gdalwarp(out_raster1, out_raster2, 'bilinear', args*)

Ist ein solcher Ansatz möglich?

djq
quelle
1
gdal_translate und gdalwarp sowie andere gdal-Dienstprogramme sind keine Python-Pakete / -Module, sondern eigenständige ausführbare Dateien. Sie können os.system oder (vorzugsweise) subprocess.Popen verwenden, um sie aufzurufen.
user2856
@ Luke, gibt es also keine Möglichkeit, mit diesen GDAL-Dienstprogrammen zu interagieren? Ich würde gerne untersuchen, wie man einen Python-Wrapper für diese schreibt, wenn es keinen gibt.
DJQ
Sie können die gdal-Python-API verwenden, um fast alles zu tun, was die Dienstprogramme gdal_translate und gdalwarp können. Es ist nur etwas komplizierter. Schauen Sie sich gdal.AutoCreateWarpedVRT und gdal Driver.CreateCopy an.
user2856

Antworten:

17

Hier ist ein Beispiel, das ungefähr das tut, wonach Sie fragen. Die Hauptparameter sind das geotransformArray, mit dem gdal eine Rasterposition (Position, Pixelskala und Versatz) und den epsgCode der Projektion beschreibt. Damit sollte der folgende Code das Raster ordnungsgemäß georeferenzieren und seine Projektion angeben.

Ich habe nicht so viel getestet, aber es schien für mich zu funktionieren. Sie müssten sicherstellen, dass die von mir eingegebenen Koordinaten korrekt sind, und wahrscheinlich die Projektion und die Pixelskala / -größe ändern. Ich hoffe es hilft.

from osgeo import gdal, osr

src_filename ='/path/to/source.tif'
dst_filename = '/path/to/destination.tif'

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Specify raster location through geotransform array
# (uperleftx, scalex, skewx, uperlefty, skewy, scaley)
# Scale = size of one pixel in units of raster projection
# this example below assumes 100x100
gt = [-7916400, 100, 0, 5210940, 0, -100]

# Set location
dst_ds.SetGeoTransform(gt)

# Get raster projection
epsg = 3857
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()

# Set projection
dst_ds.SetProjection(dest_wkt)

# Close files
dst_ds = None
src_ds = None
gelbe Kappe
quelle
2
Es ist richtig, aber wenn Sie im Falle einer Rotation genauer sein möchten, müssen Sie die Affine Transformation verwenden und die richtigen Parameter mit dem Affine-Paket berechnen (hier herunterladen ). Verwendung: 'Affine.scale (PARAMETER) * Affine.rotation (ANGLE)'. PARAMETER sind aus dem Pixelverhältnis der beiden Bilder, Sie können 'gdalinfo' oder PIL verwenden, um die Pixelgröße zu kennen, ANGLE ist der Winkel.
CaMa