Mengubah proyeksi geoTiff ke WGS84 dengan GDAL dan Python

16

Maaf jika pertanyaan berikut ini agak bodoh, tapi saya SANGAT baru dalam hal GIS ini.

Saya mencoba untuk mengubah beberapa gambar geoTiff yang diproyeksikan ke WGS84 menggunakan gdal di python. Saya telah menemukan pos yang menguraikan proses untuk mengubah poin dalam GeoTiff yang diproyeksikan menggunakan sesuatu yang mirip dengan yang berikut:

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) 

Pertanyaan saya adalah, jika saya ingin mengonversi titik-titik ini dan membuat file GeoTiff WGS84 baru, apakah ini cara terbaik untuk melakukannya? Apakah ada fungsi yang ada yang akan melakukan tugas seperti dalam 1 langkah?

Terima kasih!

JT.WK
sumber

Jawaban:

21

Satu cara yang lebih sederhana adalah dengan menggunakan alat baris perintah GDAL:

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

Itu bisa dipanggil cukup mudah melalui scripting untuk pekerjaan batch. Ini akan membelokkan raster ke kisi baru, dan ada opsi lain untuk mengontrol detailnya.

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

Sistem koordinat target (t_srs) dapat berupa PROJ.4 seperti yang ditunjukkan, atau file aktual dengan WKT, di antara opsi lain. Sistem koordinat sumber (-s_srs) diasumsikan diketahui, tetapi dapat diatur secara eksplisit seperti target.

Itu mungkin lebih mudah untuk menyelesaikan pekerjaan jika Anda tidak harus menggunakan API GDAL melalui Python.

Ada tutorial untuk warping dengan API di sini, dikatakan dukungan untuk Python tidak lengkap (saya tidak tahu detailnya): http://www.gdal.org/warptut.html

mdsumner
sumber
1
Terima kasih atas jawaban hebat mdsumner! Sementara solusi yang saya cari seharusnya ditulis dalam python, mungkin saya bisa lolos dengan hanya mengulang perintah ini dalam skrip python.
JT.WK
16

Seperti kata mdsumner, jauh lebih mudah untuk menggunakan baris perintah daripada binding python, kecuali jika Anda ingin mengeksekusi sangat tugas yang kompleks.

Jadi, jika Anda suka python, seperti saya, Anda dapat menjalankan alat baris perintah dengan:

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

atau beralih melalui daftar file:

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

Dan bahkan menggunakan alat multiprosesing menggunakan kekuatan penuh mesin Anda untuk menjalankan tugas besar.

Pablo
sumber
Terima kasih Pablo. Alat multi-pemrosesan adalah sesuatu yang saya pasti akan melihat ke dalamnya!
JT.WK
Di gdal 2.0 ada dukungan asli untuk multi-pemrosesan - tambahkan saja -DUA NUM_THREADS = ALL_CPUS ke perintah gdalwarp Anda.
valveLondon
3
os.sysadalah modul bawaan; Anda menginginkan os.systemperintah
Mike T
1
Jawaban saya sudah ketinggalan zaman, subprocess.call adalah pendekatan yang lebih baik.
Pablo