Konvertieren von projiziertem geoTiff nach WGS84 mit GDAL und Python

16

Entschuldigung, wenn die folgende Frage etwas dumm ist, aber ich bin nur SEHR neu in dieser ganzen GIS-Sache.

Ich versuche, einige projizierte GeoTiff-Bilder mit gdal in Python nach WGS84 zu konvertieren. Ich habe einen Beitrag gefunden, in dem der Prozess zum Transformieren von Punkten in projizierten GeoTiffs mit etwas ähnlichem wie dem Folgenden beschrieben wird:

from osgeo import osr, gdal

# get the existing coordinate system
ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs) 

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 

#get the coordinates in lat long
latlong = transform.TransformPoint(x,y) 

Meine Frage ist, ob ich diese Punkte konvertieren und eine neue WGS84-GeoTiff-Datei erstellen möchte. Ist dies der beste Weg, dies zu tun? Gibt es eine Funktion, die solche Aufgaben in einem Schritt erledigt?

Vielen Dank!

JT.WK
quelle

Antworten:

21

Eine einfachere Möglichkeit wäre die Verwendung der GDAL-Befehlszeilentools:

gdalwarp infile.tif outfile.tif -t_srs "+proj=longlat +ellps=WGS84"

Das kann einfach genug per Scripting für Batch-Jobs aufgerufen werden. Hierdurch wird das Raster in ein neues Raster verschoben, und es gibt andere Optionen zum Steuern der Details.

http://www.gdal.org/gdalwarp.html

Das Zielkoordinatensystem (t_srs) kann unter anderem PROJ.4 wie gezeigt oder eine tatsächliche Datei mit WKT sein. Das Quellkoordinatensystem (-s_srs) wird als bekannt vorausgesetzt, kann jedoch wie das Ziel explizit festgelegt werden.

Dies ist möglicherweise einfacher, wenn Sie die GDAL-API nicht über Python verwenden müssen.

Es gibt hier ein Tutorial zum Warping mit der API. Es besagt, dass die Unterstützung für Python unvollständig ist (ich kenne die Details nicht): http://www.gdal.org/warptut.html

mdsumner
quelle
1
Danke für die tolle Antwort mdsumner! Während die Lösung, nach der ich strebe, in Python geschrieben werden soll, kann ich vielleicht durch einfaches Schleifen dieses Befehls in einem Python-Skript davonkommen.
JT.WK
16

Wie mdsumner sagte, ist es viel einfacher, die Befehlszeile als die Python-Bindungen zu verwenden, es sei denn, Sie möchten sehr komplexe Aufgaben ausführen .

Wenn Sie Python mögen, wie ich, können Sie das Kommandozeilen-Tool folgendermaßen ausführen:

import os  
os.sys('gdalwarp infile.tif outfile.tif -t_srs "+proj=longlat +ellps=WGS84"')

oder durchlaufen Sie eine Liste von Dateien:

listoffiles=['file1, file2, file3, filen']  
for file in listoffiles:
    os.system('gdalwarp %s %s -t_srs "+proj=longlat +ellps=WGS84"' %(file,'new_'+file))  

Und selbst Multiprocessing- Tools nutzen die volle Leistung Ihrer Maschine, um große Aufgaben auszuführen.

Pablo
quelle
Vielen Dank, Pablo. Das Multiprocessing-Tool werde ich mir auf jeden Fall ansehen!
JT.WK
In gdal 2.0 gibt es native Unterstützung für Multiprocessing - fügen Sie einfach -wo NUM_THREADS = ALL_CPUS an Ihren gdalwarp-Befehl an.
valveLondon
3
os.sysist ein eingebautes Modul; Sie wollen den os.systemBefehl
Mike T
1
Meine Antwort ist veraltet, subprocess.call ist ein besserer Ansatz.
Pablo