Разбить растр на более мелкие куски, используя GDAL?

У меня есть растр (на самом деле USGS DEM), и мне нужно разделить его на более мелкие фрагменты, как показано на рисунке ниже. Это было достигнуто в ArcGIS 10.0 с использованием инструмента Split Raster. Я хотел бы метод FOSS, чтобы сделать это. Я посмотрел на GDAL, думая, что он обязательно это сделает (каким-то образом с gdal_translate), но ничего не могу найти. В конечном счете, я хотел бы иметь возможность взять растр и сказать, на какой большой размер (4KM на 4KM) я бы хотел, чтобы он разделился.

gdal_translate будет работать с использованием параметров -srcwin или -projwin.

-srcwin xoff yoff xsize ysize: выбирает подокно из исходного изображения для копирования на основе местоположения пикселей / линий.

-projwin ulx uly lrx lry: выбирает подокно из исходного изображения для копирования (например, -srcwin), но с углами, заданными в координатах с географической привязкой.

Вам нужно будет найти расположение пикселей или линий или угловые координаты, а затем перебрать значения с помощью gdal_translate. Что-то вроде быстрого и грязного питона, приведенного ниже, будет работать, если использовать значения пикселей и -srcwin подходит для вас, будет немного больше работы, чтобы разобраться с координатами.

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, считывает растровые измерения из самого файла и покрывает все изображение, уменьшая при необходимости краевые плитки:

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)

Для связывания растров в комплекте есть скрипт на python : 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

например:

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

Для @ Аарона, который спросил:

Я надеюсь найти gdalwarp версию ответа @ wwnick, которая использует опцию -multi для улучшенных многоядерных и многопоточных операций.

Небольшой отказ от ответственности

Это использует gdalwarp, хотя я не совсем уверен, что будет большой прирост производительности. До сих пор я был связан с вводом / выводом - выполнение этого сценария на большом растре, разрезании его на множество более мелких частей, не требует большой загрузки процессора, поэтому я предполагаю, что узким местом является запись на диск. Если вы планируете одновременно перепроектировать плитки или что-то подобное, это может измениться. Есть настройки советы здесь . Короткая игра не принесла мне никакого улучшения, и процессор никогда не казался ограничивающим фактором.

Отказ от ответственности, вот скрипт, который будет использоваться gdalwarpдля разделения растра на несколько меньших тайлов. Может быть некоторая потеря из-за разделения пола, но об этом можно позаботиться, выбрав желаемое количество плиток. Это будет n+1где nэто число разделить на , чтобы получить tile_widthи 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))

Вы можете использовать рутил GRASS GIS. r.tile генерирует отдельную растровую карту для каждой плитки с пронумерованными именами карт на основе определенного пользователем префикса. Ширина плиток (столбцов) и высота плиток (рядов) могут быть определены.

При использовании Python API для сессии с использованием травы требуется всего несколько строк кода Python для вызова функции r.tile извне, то есть для написания автономного скрипта. Использование r.external и r.external.out также не приводит к дублированию данных на этапе обработки GRASS GIS.

Псевдокод:

1. инициализировать сессию травы
2. определить выходной формат с помощью r.external.out
3. связать входной файл с r.external
4. запустить r.tile, который генерирует плитки в формате, определенном выше
5. отсоединить r.external.out
6. закрывать траву