Conversión de geoTiff proyectado a WGS84 con GDAL y Python

Question

Pido disculpas si la siguiente pregunta es algo estúpida, pero soy MUY nuevo en todo esto del SIG.

Estoy intentando convertir algunas imágenes geoTiff proyectadas a WGS84 utilizando gdal en python. He encontrado un post que describe el proceso para transformar puntos dentro de GeoTiffs proyectados usando algo similar a lo siguiente:

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) 

Mi pregunta es, si quiero convertir estos puntos y crear un nuevo archivo WGS84 GeoTiff, ¿es esta la mejor manera de hacerlo? ¿Existe alguna función que haga esta tarea en un solo paso?

Gracias.

Answer 1

Una forma más sencilla sería utilizar las herramientas de línea de comandos de GDAL:

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

Esto se puede invocar fácilmente mediante scripts para trabajos por lotes. Esto deformará la trama en una nueva cuadrícula, y hay otras opciones para controlar los detalles.

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

El sistema de coordenadas destino (t_srs) puede ser PROJ.4 como se muestra, o un archivo real con WKT, entre otras opciones. El sistema de coordenadas de origen (-s_srs) se supone conocido, pero puede establecerse explícitamente como el de destino.

Eso podría ser más fácil de hacer el trabajo si no tienes que usar la API de GDAL a través de Python.

Hay un tutorial para hacer warping con la API aquí, dice que el soporte para Python es incompleto (no conozco los detalles): http://www.gdal.org/warptut.html

Answer 2

Como dijo mdsumner, es mucho más fácil usar la línea de comandos que los bindings de python, a menos que quieras ejecutar muy tareas complejas.

Así que, si te gusta python, como a mí, puedes ejecutar la herramienta de línea de comandos con:

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

o iterar a través de una lista de archivos:

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

E incluso utilizar multiprocesamiento utilizan toda la potencia de su máquina para ejecutar grandes tareas.