使用GDAL将栅格分割成较小的块？

我有一个栅格（实际上是USGS DEM），我需要将其分成较小的块，如下图所示。这是在ArcGIS 10.0中使用“分割栅格”工具完成的。我想要一个FOSS方法来做到这一点。我看过GDAL，可以肯定会这样做（用gdal_translate以某种方式），但找不到任何东西。最终，我希望能够拍摄栅格并说出它要分成多大（4KM x 4KM块）。

我有一个使用subprocess.Popen的实用程序，可同时运行多个gdal转换，用于使用渔网将大栅格提取到图块，特别是在输入和/或输出高度压缩（例如LZW或放气GeoTiff）时特别有用），如果两者都不被高度压缩，则该过程将在HDD访问上达到高峰，并且不会比一次运行一个过程快得多。不幸的是，由于严格的命名约定，它的通用性不足以共享，但无论如何还是值得深思的。 GDALWarp的-multi选项通常会引起麻烦，并且仅使用2个线程（一次读取，一次写入），而并非全部可用。

#1 楼

gdal_translate将使用-srcwin或-projwin选项工作。

-srcwin xoff yoff xsize ysize：
从源图像中选择一个子窗口以基于像素/行的位置进行复制。

-projwin ulx uly lrx lry：从源图像中选择一个子窗口
进行复制（如-srcwin），但其转角以地理参考
坐标给出。 >

您需要提供像素/线位置或角坐标，然后使用gdal_translate遍历这些值。如果使用像素值并且-srcwin适合您，则类似于下面的快速且肮脏的python会起作用，使用坐标进行排序会更加麻烦。

 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)

嗨，当我尝试对带有Geotiff图像的-projwin选项进行警告时，会收到警告：“警告：计算出的-srcwin -3005000 1879300 50 650完全落在栅格范围之外。但是，继续进行下去”我不确定我在哪里做错了，似乎并没有使用其地理参考坐标。

–ncelik
16-10-26在9:06

@ncelik这可能是因为您在projwin中使用单元格坐标，而应该使用srcwin。如果您遇到困难，请发布一个包含所有相关信息的新问题，以便我们为您的具体问题提供建议。

–迈克尔·斯廷森（Michael Stimson）
18-10-18在4:52

#2 楼

我的解决方案基于@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)

我认为应该是sys.argv [1]，上面写着sys.argv [2]，对吗？

– oskarlin
17年7月30日17:00

我相信sys.argv [2]用作输出文件的前缀。超级有帮助-感谢@Ries！

– Charlie Lefrak
17年8月8日在14:40

#3 楼

有一个专门用于重新整理栅格的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

#4 楼

对于@Aaron，他问：

我希望找到@wwnick答案的gdalwarp版本，该版本利用-multi选项来增强多核和多线程操作

轻微的免责声明

它使用gdalwarp，尽管我不完全相信会有很大的性能提升。到目前为止，我一直受I / O限制-在较大的栅格上运行此脚本并将其切成许多较小的部分，似乎并不占用大量CPU，因此我认为瓶颈正在写入磁盘。如果您打算同时重新投影图块或类似图像，则可能会发生变化。这里有调优技巧。简短的播放对我没有任何改善，CPU似乎从来都不是限制因素。

除免责声明外，这里的脚本将使用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))

#5 楼

您可以使用GRASS GIS的r.tile。 r.tile根据用户定义的前缀为每个带有编号地图名称的图块生成单独的栅格地图。可以定义图块的宽度（列）和图块的高度（行）。

使用草期Python API，只需几行Python代码即可从外部调用r.tile功能，即编写一个独立的脚本。
使用r.external和r.external.out，在GRASS GIS处理步骤中也不会发生数据重复。

伪代码：

使用r.external.out初始化草稿会话
定义输出格式
带有r.external的链接输入文件
运行r.tile，以这种格式生成图块上面定义的
取消链接r.external.out
关闭grass-session

#6 楼

import gdal
import os

def create_tiles(tile_size, input_filename, in_path='/media/Data/avinash/input/',
           
out_path='/home/nesac/PycharmProjects/GIS_Data_Automation/data/tiles/'):

ds = gdal.Open(in_path + input_filename)

band = ds.GetRasterBand(1)

output_filename = 'tile_'

xsize, ysize = (band.XSize, band.YSize)
tile_size_x, tile_size_y = tile_size
tile_list = {}

complete_x = xsize // tile_size_x
complete_y = ysize // tile_size_y

residue_x = xsize % tile_size_x
residue_y = ysize % tile_size_y

# for part A
for i in range(complete_x):
    for j in range(complete_y):
        Xmin = i * tile_size_x
        Xmax = i * tile_size_x + tile_size_x - 1
        Ymin = j * tile_size_y
        Ymax = j * tile_size_y + tile_size_y - 1
        # do patch creation here
        com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
            tile_size_x) + ", " + str(tile_size_y) + " " + \
                     str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
            Xmin) + "_" + str(Ymin) + ".tif"
        os.system(com_string)

x_residue_count = 1
y_residue_count = 1

# for part B
for j in range(complete_y):
    Xmin = tile_size_x * complete_x
    Xmax = tile_size_x * complete_x + residue_x - 1
    Ymin = j * tile_size_y
    Ymax = j * tile_size_y + tile_size_y - 1
    # do patch creation here
    com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
        residue_x) + ", " + str(tile_size_y) + " " + \
                 str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
        Xmin) + "_" + str(Ymin) + ".tif"
    os.system(com_string)

# for part C
for i in range(complete_x):
    Xmin = i * tile_size_x
    Xmax = i * tile_size_x + tile_size_x - 1
    Ymin = tile_size_y * complete_y
    Ymax = tile_size_y * complete_y + residue_y - 1
    com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
        tile_size_x) + ", " + str(residue_y) + " " + \
                 str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
        Xmin) + "_" + str(Ymin) + ".tif"
    os.system(com_string)

# for part D
Xmin = complete_x * tile_size_x
Ymin = complete_y * tile_size_y
com_string = "gdal_translate -of GTIFF -srcwin " + str(Xmin) + ", " + str(Ymin) + ", " + str(
    residue_x) + ", " + str(residue_y) + " " + \
             str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(
    Xmin) + "_" + str(Ymin) + ".tif"
os.system(com_string)

编程黑洞网