Преобразование спроецированного геоТиффа в WGS84 с помощью GDAL и Python

16

Извиняюсь, если следующий вопрос несколько глуп, но я только ОЧЕНЬ новичок во всей этой ГИС.

Я пытаюсь преобразовать некоторые проецируемые изображения GeoTiff в WGS84, используя gdal в python. Я нашел сообщение, в котором описывается процесс преобразования точек в проецируемых геотиффах с использованием чего-то похожего на следующее:

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) 

У меня вопрос: если я хочу преобразовать эти точки и создать новый файл WGS84 GeoTiff, это лучший способ сделать это? Существует ли функция, которая будет выполнять такие задачи, как задача за 1 шаг?

Благодарность!

JT.WK
источник

Ответы:

21

Один более простой способ - использовать инструменты командной строки GDAL:

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

Это можно легко вызвать с помощью сценариев для пакетных заданий. Это переместит растр к новой сетке, и есть другие опции для управления деталями.

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

Целевая (t_srs) система координат может быть PROJ.4, как показано, или фактическим файлом с WKT, среди других опций. Предполагается, что система координат источника (-s_srs) известна, но может быть задана явно как цель.

Это может быть легче выполнить, если вам не нужно использовать GDAL API через Python.

Здесь есть руководство по деформации с помощью API, в котором говорится, что поддержка Python неполна (я не знаю подробностей): http://www.gdal.org/warptut.html

mdsumner
источник
1
Спасибо за отличный ответ mdsumner! Хотя решение, которое я ищу, должно быть написано на python, возможно, мне удастся просто зациклить эту команду в скрипте python.
JT.WK
16

Как сказал mdsumner, гораздо проще использовать командную строку, чем привязки python, если вы не хотите выполнять очень сложные задачи.

Итак, если вам нравится python, как и я, вы можете запустить инструмент командной строки с помощью:

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

или перебрать список файлов:

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

И даже использовать многопроцессорные инструменты используют всю мощь вашей машины для выполнения больших задач.

Pablo
источник
Спасибо, Пабло. Мультипроцессорные инструменты - это то, что я обязательно изучу!
JT.WK
В gdal 2.0 есть встроенная поддержка многопроцессорности - просто добавьте -wo NUM_THREADS = ALL_CPUS к вашей команде gdalwarp.
valveLondon
3
os.sysэто встроенный модуль; Вы хотите os.systemкоманду
Майк Т
1
Мой ответ устарел, subprocess.call - лучший подход.
Пабло