Gdal: obcinanie rastra innym rastrem

14

Piszę proste narzędzie do przycinania partii wielopasmowych plików rastrowych geotiff do tego samego (mniejszego) obszaru. Korzystając z gdalwarp, mogę z łatwością przyciąć plik za pomocą pliku kształtowego przycinającego jeden wielokąt:

gdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif

Jednak rzeczywisty obszar, do którego chcę przyciąć, będzie zawsze początkowo definiowany przez inny plik rastrowy geotiff, a nie plik kształtu. Byłoby miło, gdybym mógł użyć zakresu tego rastra jako pliku przycinającego, ale nie jestem pewien, jak to zrobić. Nic dziwnego, że następujące nie działa (nie powoduje błędu, po prostu nic nie produkuje):

gdalwarp -cutline clipper.tif-crop_to_cutline input.tif output.tif

Moje pytanie brzmi: czy istnieje sposób na dostarczenie rastra gdalwarp -cutline? Alternatywnie, czy istnieje inna funkcja gdal, która może przyciąć raster przy użyciu innego rastra? Jeśli żadne z nich nie jest możliwe, czy istnieje bardzo prosty sposób na utworzenie pliku kształtu z pojedynczym wielokątem zdefiniowanym przez zasięg rastra?

Ten kod zostanie zawinięty w bardziej rozbudowany skrypt Pythona, więc mogę używać albo narzędzi gdal z wiersza poleceń, albo dowolnego powiązania Pythona dla gdal.

Na marginesie, wiem, że z łatwością mogę po prostu wyciąć plik kształtu, który obejmuje zasięg mojego rastra w QGIS. Mogę to zrobić, jeśli nie znajdę prostego rozwiązania, ale ostatecznie skończę przy użyciu tego narzędzia na dziesiątkach, jeśli nie setkach obszarów w ramach dużej zautomatyzowanej analizy, więc wolałbym nie mieć żmudnego krok ręczny, nawet jeśli jest to bardzo łatwe.

Joe
źródło

Odpowiedzi:

11

Nie wiem, czy można przyciąć raster innym rastrem, ale możesz użyć gdaltindex do zbudowania pliku kształtu z rozmiarem twojego rastra.

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

lejedi76
źródło
4
gdaltindex działa idealnie, aby utworzyć przycinający plik kształtu z mojego początkowego rastra. Aby rozwiązać problem, używam gdaltindex clipper.shp clipper.tif, a następniegdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif
Joe
Korzystałem z tego podejścia, ale okazało się, że czasami był to jeden piksel poza wersją przyciętą. Myślę, że łatwiej jest obliczyć swój cel, a następnie użyć odpowiedzi Xaviera poniżej, a następnie użyć gdalwarp i podać -te_srs do obsługi niedopasowanych CRS.
Jon
7

W przypadku nieregularnych wielokątów i zakładając, że plik rastrowy geotiff jest plikiem binarnym, możesz użyć GDAL_Calc :

GDAL_Calc.py -A Mask.tif -B CutBigImageToClip.tif --outfile=SmallerFile.tif --NoDataValue=0 --Calc="B*(A>0)" 

To zapytanie wypełni 0, gdzie Mask.tif <= 0 i BigImage, gdzie maska> 0. Aby to zrobić, oba rastry muszą mieć ten sam rozmiar komórki, wiersze i kolumny. Aby wyodrębnić te same zakresy, użyj GDAL_Translate z -projwin ulx uly lrx lryopcją (pole ma współrzędne rzutowane), ale upewnij się, że pole projwin nie wystaje poza krawędzie żadnego rastra.

GDAL_Translate -of GTIFF -projwin ulx uly lrx lry BigImageToClip.tif CutBigImageToClip.tif

Wartości zastępcze dla pola projwin pochodzącego z Maski.

Michael Stimson
źródło
1
+1 To jest przydatna informacja, ale myślę, że mogę rozwiązać problem w kilku krokach, korzystając z odpowiedzi @ lejedi.
Joe
4

Rozwiązanie bezpośrednio w Pythonie, bez tworzenia kształtu:

import gdal
from gdalconst import GA_ReadOnly

data = gdal.Open('img_mask.tif', GA_ReadOnly)
geoTransform = data.GetGeoTransform()
minx = geoTransform[0]
maxy = geoTransform[3]
maxx = minx + geoTransform[1] * data.RasterXSize
miny = maxy + geoTransform[5] * data.RasterYSize
call('gdal_translate -projwin ' + ' '.join([str(x) for x in [minx, maxy, maxx, miny]]) + ' -of GTiff img_orig.tif img_out.tif', shell=True)
XavierCLL
źródło
1
Uwaga: to rozwiązanie działa tylko wtedy, gdy są w tej samej SRS.
Skylion
@ Skylion Ale możesz łatwo to wyjaśnić, włączając opcję -te_srs, chociaż zamiast tego musisz także użyć gdalwarp z opcją -te.
Jon