Czy dzielisz raster na mniejsze części za pomocą GDAL?

18

Mam raster (właściwie USGS DEM) i muszę go podzielić na mniejsze części, jak pokazano na poniższym obrazku. Dokonano tego w ArcGIS 10.0 za pomocą narzędzia Split Raster. Chciałbym, aby to zrobić metoda FOSS. Patrzyłem na GDAL, myśląc, że na pewno by to zrobił (jakoś z gdal_translate), ale nic nie mogę znaleźć. Ostatecznie chciałbym móc wziąć raster i powiedzieć, jak duży (fragmenty 4KM po 4KM) chciałbym, aby się podzielił.

wprowadź opis zdjęcia tutaj

Chad Cooper
źródło
Mam narzędzie, które korzysta z podprocesu. Otwórz, aby uruchomić wiele tłumaczeń gdal w tym samym czasie, którego używam do wyodrębnienia dużego rastra na kafelki za pomocą kabaretki, szczególnie przydatne, jeśli dane wejściowe i / lub wyjściowe są silnie skompresowane (np. LZW lub deflate GeoTiff ), jeśli żadna z nich nie jest wysoce skompresowana, proces osiąga szczyt w dostępie do dysku twardego i nie jest dużo szybszy niż uruchamianie pojedynczo. Niestety nie jest wystarczająco ogólny, aby się nim dzielić z powodu sztywnych konwencji nazewnictwa, ale i tak jedzenie do przemyślenia. Opcja -multi dla GDALWarp często powoduje problemy i używa tylko 2 wątków (jeden odczyt, jeden zapis), nie wszystkie są dostępne.
Michael Stimson,

Odpowiedzi:

18

gdal_translate będzie działać przy użyciu opcji -srcwin lub -projwin.

-srcwin xoff yoff xsize ysize: Wybiera okno podrzędne z obrazu źródłowego do kopiowania na podstawie położenia piksela / linii.

-projwin ulx uly lrx lry: Wybiera okno podrzędne z obrazu źródłowego do kopiowania (jak -srcwin), ale z narożnikami podanymi we współrzędnych georeferencyjnych.

Musisz wymyślić lokalizacje pikseli / linii lub współrzędne narożne, a następnie zapętlić wartości za pomocą gdal_translate. Coś takiego jak szybki i brudny python poniżej zadziała, jeśli użycie wartości pikseli i -srcwin jest dla ciebie odpowiedni, będzie trochę więcej pracy, aby uporządkować współrzędne.

import os, gdal
from gdalconst import *

width = 512
height = 512
tilesize = 64

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(tilesize)+", " \
            +str(tilesize)+" utm.tif utm_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)
wwnick
źródło
1
Cześć, kiedy wypróbuję opcję -projwin z obrazem geotiff, pojawia się ostrzeżenie: „Ostrzeżenie: Computing -srcwin -3005000 1879300 50 650 jest całkowicie poza zasięgiem rastrowym. Kontynuowanie jednak” Nie jestem pewien, gdzie robię źle, wygląda na to, że nie używając współrzędnych georeferencyjnych.
ncelik
@ncelik prawdopodobnie dlatego, że używasz współrzędnych komórki w swoim projwinie i zamiast tego powinieneś użyć srcwin. Jeśli masz trudności, napisz nowe pytanie ze wszystkimi istotnymi informacjami, abyśmy mogli zaproponować konkretny problem.
Michael Stimson,
15

Moje rozwiązanie, oparte na tym z @wwnick, odczytuje wymiary rastrowe z samego pliku i obejmuje cały obraz, zmniejszając w razie potrzeby kafelki krawędzi:

import os, sys
from osgeo import gdal

dset = gdal.Open(sys.argv[1])

width = dset.RasterXSize
height = dset.RasterYSize

print width, 'x', height

tilesize = 5000

for i in range(0, width, tilesize):
    for j in range(0, height, tilesize):
        w = min(i+tilesize, width) - i
        h = min(j+tilesize, height) - j
        gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(w)+", " \
            +str(h)+" " + sys.argv[1] + " " + sys.argv[2] + "_"+str(i)+"_"+str(j)+".tif"
        os.system(gdaltranString)
Ries
źródło
Myślę, że powinien to być sys.argv [1], gdzie jest napisane sys.argv [2], prawda?
oskarlin
3
Uważam, że sys.argv [2] jest używany jako przedrostek dla plików wyjściowych. Super pomocny - dzięki @Ries!
Charlie Hofmann
4

W pakiecie znajduje się skrypt Pythona specjalnie do przeskalowywania rastrów, gdal_retile :

gdal_retile.py [-v] [-co NAME=VALUE]* [-of out_format] [-ps pixelWidth pixelHeight]
               [-overlap val_in_pixel]
               [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
                      CInt16/CInt32/CFloat32/CFloat64}]'
               [ -tileIndex tileIndexName [-tileIndexField tileIndexFieldName]]
               [ -csv fileName [-csvDelim delimiter]]
               [-s_srs srs_def]  [-pyramidOnly]
               [-r {near/bilinear/cubic/cubicspline/lanczos}]
               -levels numberoflevels
               [-useDirForEachRow]
               -targetDir TileDirectory input_files

na przykład:

gdal_retile.py -ps 512 512 -targetDir C:\example\dir some_dem.tif

mikewatt
źródło
4

Dla @Aaron, który zapytał:

Mam nadzieję znaleźć wersję gdalwarp odpowiedzi @ wwnick, która wykorzystuje opcję -multi do ulepszonych operacji wielordzeniowych i wielowątkowych

Nieznaczne zrzeczenie się odpowiedzialności

Wykorzystuje to gdalwarp, chociaż nie jestem do końca przekonany, że będzie znaczny wzrost wydajności. Do tej pory byłem związany z operacjami wejścia / wyjścia - uruchomienie tego skryptu na dużym rastrze pocięcie go na wiele mniejszych części nie wydaje się bardzo obciążające procesor, więc zakładam, że wąskim gardłem jest zapisywanie na dysk. Jeśli planujesz jednoczesne ponowne wyświetlanie kafelków lub czegoś podobnego, to może się to zmienić. Istnieją wskazówki strojenia tutaj . Krótka gra nie przyniosła mi żadnej poprawy, a procesor nigdy nie wydawał się czynnikiem ograniczającym.

Poza zrzeczeniem się odpowiedzialności, oto skrypt, który posłuży gdalwarpdo podzielenia rastra na kilka mniejszych płytek. Może wystąpić pewna strata z powodu podziału podłogi, ale można temu zaradzić, wybierając odpowiednią liczbę płytek. Będzie to n+1gdzie njest numer podzielić przez uzyskać tile_widthi tile_heightzmienne.

import subprocess
import gdal
import sys


def gdalwarp(*args):
    return subprocess.check_call(['gdalwarp'] + list(args))


src_path = sys.argv[1]
ds = gdal.Open(src_path)

try:
    out_base = sys.argv[2]
except IndexError:
    out_base = '/tmp/test_'

gt = ds.GetGeoTransform()

width_px = ds.RasterXSize
height_px = ds.RasterYSize

# Get coords for lower left corner
xmin = int(gt[0])
xmax = int(gt[0] + (gt[1] * width_px))

# get coords for upper right corner
if gt[5] > 0:
    ymin = int(gt[3] - (gt[5] * height_px))
else:
    ymin = int(gt[3] + (gt[5] * height_px))

ymax = int(gt[3])

# split height and width into four - i.e. this will produce 25 tiles
tile_width = (xmax - xmin) // 4
tile_height = (ymax - ymin) // 4

for x in range(xmin, xmax, tile_width):
    for y in range(ymin, ymax, tile_height):
        gdalwarp('-te', str(x), str(y), str(x + tile_width),
                 str(y + tile_height), '-multi', '-wo', 'NUM_THREADS=ALL_CPUS',
                 '-wm', '500', src_path, out_base + '{}_{}.tif'.format(x, y))
RoperMaps
źródło
3

Możesz użyć r. Płytki GRASS GIS. r.tile generuje osobną mapę rastrową dla każdej płytki z numerowanymi nazwami map na podstawie prefiksu zdefiniowanego przez użytkownika. Można zdefiniować szerokość płytek (kolumny) i wysokość płytek (wiersze).

Używając trawiastego interfejsu API Python wystarczy kilka linii kodu Python, aby wywołać funkcjonalność r.tile z zewnątrz, tj. Napisać samodzielny skrypt. Użycie r.external i r.external.out również nie powoduje duplikacji danych podczas etapu przetwarzania GRASS GIS.

Pseudo kod:

  1. zainicjuj sesję trawiastą
  2. zdefiniuj format wyjściowy za pomocą r.external.out
  3. połącz plik wejściowy z r.external
  4. uruchom r.tile, który generuje kafelki w zdefiniowanym powyżej formacie
  5. rozłącz r.external.out
  6. zamknij sesję trawiastą
markusN
źródło