Conversion de geoTiff projeté en WGS84 avec GDAL et Python

16

Toutes mes excuses si la question suivante est quelque peu stupide, mais je ne suis que TRÈS novice dans ce domaine SIG.

J'essaie de convertir des images geoTiff projetées en WGS84 en utilisant gdal en python. J'ai trouvé un article qui décrit le processus de transformation des points dans les GeoTiffs projetés en utilisant quelque chose de similaire à ce qui suit:

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) 

Ma question est, si je veux convertir ces points et créer un nouveau fichier GeoTiff WGS84, est-ce la meilleure façon de procéder? Existe-t-il une fonction qui fera comme une tâche en 1 étape?

Merci!

JT.WK
la source

Réponses:

21

Une façon plus simple serait d'utiliser les outils de ligne de commande GDAL:

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

Cela peut être invoqué assez facilement via un script pour les travaux par lots. Cela déformera le raster vers une nouvelle grille, et il existe d'autres options pour contrôler les détails.

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

Le système de coordonnées cible (t_srs) peut être PROJ.4 comme indiqué, ou un fichier réel avec WKT, entre autres options. Le système de coordonnées source (-s_srs) est supposé connu, mais peut être défini explicitement comme la cible.

Cela pourrait être plus facile de faire le travail si vous n'avez pas à utiliser l'API GDAL via Python.

Il y a un tutoriel pour la déformation avec l'API ici, il dit que le support pour Python est incomplet (je ne connais pas les détails): http://www.gdal.org/warptut.html

mdsumner
la source
1
Merci pour la grande réponse mdsumner! Alors que la solution que je recherche est censée être écrite en python, je peux peut-être m'en tirer en bouclant simplement cette commande dans un script python.
JT.WK
16

Comme l'a dit mdsumner, il est beaucoup plus facile d'utiliser la ligne de commande que les liaisons python, sauf si vous souhaitez exécuter des tâches très complexes.

Donc, si vous aimez python, comme moi, vous pouvez exécuter l'outil de ligne de commande avec:

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

ou parcourez une liste de fichiers:

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

Et même utiliser des outils de multitraitement utiliser toute la puissance de votre machine pour exécuter de grandes tâches.

Pablo
la source
Merci Pablo. Les outils de multitraitement sont quelque chose que je vais certainement étudier!
JT.WK
Dans gdal 2.0, il existe une prise en charge native du multitraitement - il suffit d'ajouter -wo NUM_THREADS = ALL_CPUS à votre commande gdalwarp.
valveLondon
3
os.sysest un module intégré; vous voulez la os.systemcommande
Mike T
1
Ma réponse est obsolète, subprocess.call est une meilleure approche.
Pablo