Memisahkan raster menjadi potongan yang lebih kecil menggunakan GDAL?

Saya memiliki raster (sebenarnya USGS DEM) dan saya perlu membaginya menjadi potongan-potongan kecil seperti gambar di bawah ini. Itu dicapai dalam ArcGIS 10.0 menggunakan alat Split Raster. Saya ingin metode FOSS untuk melakukan ini. Saya telah melihat GDAL, berpikir pasti itu akan melakukannya (entah bagaimana dengan gdal_translate), tetapi tidak dapat menemukan apa pun. Pada akhirnya, saya ingin dapat mengambil raster dan mengatakan seberapa besar (4KM dengan 4KM potongan) Saya ingin dibagi menjadi.

gdal_translate akan bekerja menggunakan opsi -srcwin atau -projwin.

-srcwin xoff yoff xsize ysize: Memilih subwindow dari gambar sumber untuk disalin berdasarkan lokasi piksel / garis.

-projwin ulx uly lrx lry: Memilih subwindow dari gambar sumber untuk disalin (seperti -srcwin) tetapi dengan sudut yang diberikan dalam koordinat georeferensi.

Anda harus menemukan lokasi piksel / garis atau koordinat sudut dan kemudian mengulangi nilainya dengan gdal_translate. Sesuatu seperti python cepat dan kotor di bawah ini akan berfungsi jika menggunakan nilai pixel dan -srcwin cocok untuk Anda, akan sedikit lebih banyak pekerjaan untuk diselesaikan dengan koordinat.

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)

Solusi saya, berdasarkan yang dari @wwnick membaca dimensi raster dari file itu sendiri, dan mencakup seluruh gambar dengan membuat ubin tepi lebih kecil jika diperlukan:

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)

Ada skrip python yang dibundel khusus untuk retiling raster, 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

misalnya:

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

Untuk @ Harun yang bertanya:

Saya berharap untuk menemukan versi gdalwarp dari jawaban @ wwnick yang menggunakan opsi -multi untuk operasi multicore dan multithreaded yang ditingkatkan

Penafian Sedikit

Ini menggunakan gdalwarp, meskipun saya tidak sepenuhnya yakin akan ada banyak peningkatan kinerja. Sejauh ini saya sudah terikat I / O - menjalankan skrip ini pada raster besar yang memotongnya menjadi banyak bagian yang lebih kecil tampaknya tidak intensif CPU, jadi saya menganggap hambatannya adalah menulis ke disk. Jika Anda berencana memproyeksikan ulang ubin secara bersamaan atau yang serupa, maka ini mungkin berubah. Ada tips menyetel di sini . Permainan singkat tidak menghasilkan perbaikan bagi saya, dan CPU sepertinya tidak pernah menjadi faktor pembatas.

Di samping penafian, inilah skrip yang akan digunakan gdalwarpuntuk membagi raster menjadi beberapa ubin yang lebih kecil. Mungkin ada beberapa kerugian karena pembagian lantai tetapi ini bisa diatasi dengan memilih jumlah ubin yang Anda inginkan. Ini akan menjadi n+1tempat nangka yang Anda bagi untuk mendapatkan tile_widthdan tile_heightvariabel.

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

Anda dapat menggunakan r.tile dari GRASS GIS. r.tile menghasilkan peta raster yang terpisah untuk setiap ubin dengan nama peta bernomor berdasarkan awalan yang ditentukan pengguna. Lebar ubin (kolom) dan tinggi ubin (baris) dapat ditentukan.

Menggunakan API Python sesi-rumput hanya beberapa baris kode Python yang diperlukan untuk memanggil fungsionalitas r.tile dari luar, yaitu untuk menulis skrip mandiri. Menggunakan r.external dan r.external.out juga tidak ada duplikasi data yang terjadi selama langkah pemrosesan GRASS GIS.

Kode palsu:

1. inisialisasi sesi rumput
2. tentukan format output dengan r.external.out
3. menghubungkan file input dengan r.external
4. jalankan r.tile yang menghasilkan petak dalam format yang ditentukan di atas
5. batalkan tautan r.external.out
6. tutup sesi rumput