Convertendo o geoTiff projetado para WGS84 com GDAL e Python

16

Desculpas se a pergunta a seguir é um pouco estúpida, mas eu sou muito MUITO novo nessa coisa toda de GIS.

Estou tentando converter algumas imagens geoTiff projetadas para WGS84 usando gdal em python. Eu encontrei um post que descreve o processo para transformar pontos nos GeoTiffs projetados usando algo semelhante ao seguinte:

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) 

Minha pergunta é: se eu quiser converter esses pontos e criar um novo arquivo WGS84 GeoTiff, essa é a melhor maneira de fazer isso? Existe uma função que funcione como tarefa em 1 etapa?

Obrigado!

JT.WK
fonte

Respostas:

21

Uma maneira mais simples seria usar as ferramentas de linha de comando GDAL:

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

Isso pode ser chamado com bastante facilidade por meio de script para tarefas em lote. Isso distorcerá a varredura para uma nova grade, e existem outras opções para controlar os detalhes.

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

O sistema de coordenadas de destino (t_srs) pode ser PROJ.4, como mostrado, ou um arquivo real com o WKT, entre outras opções. O sistema de coordenadas de origem (-s_srs) é considerado conhecido, mas pode ser definido explicitamente como o destino.

Pode ser mais fácil concluir o trabalho se você não precisar usar a API GDAL via Python.

Há um tutorial para distorcer a API aqui, que diz que o suporte ao Python está incompleto (não sei os detalhes): http://www.gdal.org/warptut.html

mdsumner
fonte
1
Obrigado pela ótima resposta mdsumner! Embora a solução que eu busco deva ser escrita em python, talvez eu consiga fazer um loop desse comando em um script python.
JT.WK
16

Como o mdsumner disse, é muito mais fácil usar a linha de comando do que as ligações python, a menos que você queira executar tarefas muito complexas.

Portanto, se você gosta de python, como eu, execute a ferramenta de linha de comando com:

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

ou itere através de uma lista de arquivos:

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

E mesmo as ferramentas de multiprocessamento usam todo o poder da sua máquina para executar grandes tarefas.

Pablo
fonte
Obrigado Pablo. As ferramentas de multiprocessamento são algo que eu definitivamente irei analisar!
JT.WK
No gdal 2.0, há suporte nativo para o multiprocessamento - basta anexar -nUM_THREADS = ALL_CPUS ao seu comando gdalwarp.
usar o seguinte comando
3
os.sysé um módulo embutido; você deseja que o os.systemcomando
Mike T
1
Minha resposta está desatualizada, subprocess.call é uma abordagem melhor.
Pablo