Raster mit GDAL in kleinere Chunks aufteilen?

18

Ich habe ein Raster (USGS DEM tatsächlich) und ich muss es in kleinere Stücke aufteilen, wie das Bild unten zeigt. Dies wurde in ArcGIS 10.0 mit dem Werkzeug "Raster teilen" durchgeführt. Ich möchte eine FOSS-Methode, um dies zu tun. Ich habe mir GDAL angesehen und gedacht, es würde es schaffen (irgendwie mit gdal_translate), aber ich kann nichts finden. Letztendlich möchte ich in der Lage sein, das Raster zu nehmen und zu sagen, in wie große (4KM mal 4KM große Stücke) ich es aufteilen möchte.

Bildbeschreibung hier eingeben

Chad Cooper
quelle
Ich habe ein Hilfsprogramm, das subprocess.Popen verwendet, um mehrere gdal-Übersetzungen gleichzeitig auszuführen, die ich zum Extrahieren eines großen Rasters in Kacheln mithilfe eines Fischnetzes verwende. Dies ist besonders nützlich, wenn die Eingabe und / oder Ausgabe stark komprimiert ist (z. B. LZW oder GeoTiff entleeren) ), wenn keiner der beiden stark komprimiert ist, ist der Prozess auf dem Festplattenzugriff am größten und nicht viel schneller als einer nach dem anderen. Leider ist es aufgrund strenger Namenskonventionen nicht generisch genug, um geteilt zu werden, aber es gibt trotzdem Anlass zum Nachdenken. Die Option -multi für GDALWarp verursacht häufig Probleme und verwendet nur 2 Threads (ein Lese-, ein Schreib-Thread), die nicht alle verfügbar sind.
Michael Stimson

Antworten:

18

gdal_translate funktioniert mit den Optionen -srcwin oder -projwin.

-srcwin xoff yoff xsize ysize: Wählt basierend auf der Pixel- / Linienposition ein Unterfenster aus dem Quellbild zum Kopieren aus.

-projwin ulx uly lrx lry: Wählt ein Unterfenster aus dem Quellbild zum Kopieren aus (wie -srcwin), wobei die Ecken in georeferenzierten Koordinaten angegeben sind.

Sie müssten sich die Pixel- / Linienpositionen oder Eckkoordinaten einfallen lassen und dann die Werte mit gdal_translate durchlaufen. Etwas wie das schnelle und schmutzige Python unten wird funktionieren, wenn die Verwendung von Pixelwerten und -srcwin für Sie geeignet ist, wird ein bisschen mehr Arbeit sein, um mit Koordinaten zu sortieren.

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
quelle
1
Hallo, wenn ich die Option -projwin mit einem Geotiff-Image ausprobiere, erhalte ich die Warnmeldung "Warnung: Berechnet -srcwin -3005000 1879300 50 650 liegt vollständig außerhalb des Rasterbereichs. Wenn ich jedoch weitermache", bin ich mir nicht sicher, wo ich falsch mache mit seinen georeferenzierten Koordinaten.
Ncelik
@ncelik Das liegt wahrscheinlich daran, dass Sie in Ihrem projwin Zellkoordinaten verwenden und stattdessen srcwin verwenden sollten. Wenn Sie auf Schwierigkeiten stoßen, stellen Sie bitte eine neue Frage mit allen relevanten Informationen, damit wir Vorschläge zu Ihrem spezifischen Problem machen können.
Michael Stimson
15

Meine Lösung, die auf der von @wwnick basiert, liest die Rastermaße aus der Datei selbst und deckt das gesamte Bild ab, indem die Randkacheln bei Bedarf verkleinert werden:

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
quelle
Ich denke, es sollte sys.argv [1] sein, wo sys.argv [2] steht, oder?
oskarlin
3
Ich glaube, sys.argv [2] wird als Präfix für die Ausgabedateien verwendet. Super hilfsbereit - danke @Ries!
Charlie Hofmann
4

Es gibt ein gebündeltes Python-Skript speziell für das Retiling von Rastern, 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

z.B:

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

Mikewatt
quelle
4

Für @Aaron, der gefragt hat:

Ich hoffe, eine gdalwarp-Version von @ wwnicks Antwort zu finden, die die -multi-Option für erweiterte Multicore- und Multithread-Operationen verwendet

Leichter Haftungsausschluss

Dies nutzt gdalwarp, obwohl ich nicht ganz davon überzeugt bin, dass es viel Leistungsgewinn geben wird. Bisher war ich an E / A gebunden - das Ausführen dieses Skripts auf einem großen Raster, das es in viele kleinere Teile zerlegt, scheint nicht CPU-intensiv zu sein, daher gehe ich davon aus, dass der Engpass beim Schreiben auf die Festplatte besteht. Wenn Sie vorhaben, die Kacheln oder ähnliches gleichzeitig erneut zu projizieren, kann sich dies ändern. Es gibt Tuning - Tipps hier . Ein kurzes Spiel brachte für mich keine Besserung und die CPU schien nie der limitierende Faktor zu sein.

Abgesehen vom Haftungsausschluss ist hier ein Skript, mit gdalwarpdem Sie ein Raster in mehrere kleinere Kacheln aufteilen können. Aufgrund der Bodenteilung kann es zu Verlusten kommen. Dies kann jedoch durch Auswahl der gewünschten Anzahl von Fliesen behoben werden. Es wird sein, n+1wo ndie Zahl ist, durch die Sie dividieren, um die Variablen tile_widthund zu erhalten tile_height.

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
quelle
3

Sie können r.tile verwenden von GRASS GIS verwenden. r.tile generiert eine separate Rasterkarte für jede Kachel mit nummerierten Kartennamen basierend auf dem benutzerdefinierten Präfix. Die Breite der Kacheln (Spalten) und die Höhe der Kacheln (Zeilen) können definiert werden.

Mit Hilfe des Gras-Sitzung Python API nur ein paar Zeilen Python - Code benötigt , um die r.tile Funktionalität von außen zu nennen, dh ein eigenständiges Skript zu schreiben. Bei Verwendung von r.external und r.external.out tritt während des GRASS GIS-Verarbeitungsschritts auch keine Datenverdoppelung auf.

Pseudocode:

  1. Grassitzung initialisieren
  2. Ausgabeformat mit r.external.out definieren
  3. Verknüpfung der Eingabedatei mit r.external
  4. Führen Sie r.tile aus, um die Kacheln im oben definierten Format zu generieren
  5. Verknüpfung aufheben r.external.out
  6. Grassession schließen
markusN
quelle