Convertir geoTiff proyectado a WGS84 con GDAL y Python

16

Disculpas si la siguiente pregunta es algo estúpida, pero solo soy MUY nuevo en todo este asunto SIG.

Estoy tratando de convertir algunas imágenes geoTiff proyectadas a WGS84 usando gdal en python. He encontrado una publicación 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 GeoTiff WGS84, ¿es esta la mejor manera de hacerlo? ¿Existe una función que haga tal tarea en 1 paso?

¡Gracias!

JT.WK
fuente

Respuestas:

21

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

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

Eso se puede invocar fácilmente a través de secuencias de comandos para trabajos por lotes. Esto deformará el ráster a una nueva cuadrícula, y hay otras opciones para controlar los detalles.

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

El sistema de coordenadas de 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 se puede establecer explícitamente como el destino.

Puede ser más fácil hacer el trabajo si no tiene que usar la API de GDAL a través de Python.

Aquí hay un tutorial para deformar con la API, dice que el soporte para Python está incompleto (no sé los detalles): http://www.gdal.org/warptut.html

mdsumner
fuente
1
Gracias por la gran respuesta mdsumner! Si bien se supone que la solución que busco está escrita en python, tal vez pueda saltearme este bucle en un script de python.
JT.WK
16

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

Entonces, si te gusta Python, como a mí, puedes ejecutar la herramienta de línea de comando 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 use herramientas de multiprocesamiento , use toda la potencia de su máquina para ejecutar grandes tareas.

Pablo
fuente
Gracias pablo ¡Las herramientas de multiprocesamiento son algo que definitivamente consideraré!
JT.WK
En gdal 2.0 hay soporte nativo para multiprocesamiento, solo agregue -wo NUM_THREADS = ALL_CPUS a su comando gdalwarp.
valveLondon
33
os.syses un módulo incorporado; quieres el os.systemcomando
Mike T
1
Mi respuesta está desactualizada, subprocess.call es un mejor enfoque.
Pablo