Konwersja rzutowanego geoTiff na WGS84 za pomocą GDAL i Pythona

16

Przepraszam, jeśli poniższe pytanie jest nieco głupie, ale jestem BARDZO nowy w tej całej sprawie GIS.

Próbuję przekonwertować niektóre rzutowane obrazy geoTiff na WGS84 przy użyciu gdal w python. Znalazłem post opisujący proces przekształcania punktów w rzutowanych GeoTiff przy użyciu czegoś podobnego do następującego:

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) 

Moje pytanie brzmi: czy chcę przekonwertować te punkty i utworzyć nowy plik GeoTiff WGS84? Czy istnieje funkcja, która wykona zadanie takie jak 1 krok?

Dzięki!

JT.WK
źródło

Odpowiedzi:

21

Jednym prostszym sposobem byłoby użycie narzędzi wiersza poleceń GDAL:

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

Można to łatwo wywołać za pomocą skryptów dla zadań wsadowych. Spowoduje to wypaczenie rastra do nowej siatki i istnieją inne opcje kontrolowania szczegółów.

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

Docelowym układem współrzędnych (t_srs) może być PROJ.4, jak pokazano, lub rzeczywisty plik z WKT, między innymi opcjami. Zakłada się, że układ współrzędnych źródłowy (-s_srs) jest znany, ale może być ustawiony jawnie jak cel.

Wykonanie zadania może być łatwiejsze, jeśli nie musisz korzystać z interfejsu GDAL API przez Python.

Tutaj znajduje się samouczek dotyczący wypaczania z interfejsem API, który mówi, że obsługa Pythona jest niepełna (nie znam szczegółów): http://www.gdal.org/warptut.html

mdsumner
źródło
1
Dzięki za świetną odpowiedź mdsumner! Podczas gdy rozwiązanie, którego szukam, ma być napisane w pythonie, być może uda mi się zapętlić to polecenie w skrypcie Pythona.
JT.WK
16

Jak powiedział mdsumner, o wiele łatwiej jest korzystać z wiersza poleceń niż powiązania Pythona, chyba że chcesz wykonywać bardzo złożone zadania.

Jeśli więc lubisz Pythona, tak jak ja, możesz uruchomić narzędzie wiersza poleceń za pomocą:

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

lub iteruj po liście plików:

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

Nawet narzędzia wieloprocesorowe wykorzystują pełną moc maszyny do wykonywania dużych zadań.

Pablo
źródło
Dzięki Pablo. Narzędzia do wieloprzetwarzania to coś, na co na pewno przyjrzę się!
JT.WK
W gdal 2.0 dostępna jest natywna obsługa wielu procesów - wystarczy dołączyć -wo NUM_THREADS = ALL_CPUS do polecenia gdalwarp.
valveLondon
3
os.sysjest wbudowanym modułem; chcesz os.systempolecenia
Mike T.
1
Moja odpowiedź jest nieaktualna, subprocess.call to lepsze podejście.
Pablo