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

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ł.

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)

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)

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

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))

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ą