Splitting raster into smaller chunks using GDAL?

My solution, based on the one from @wwnick reads the raster dimensions from the file itself, and covers the whole image by making the edge tiles smaller if needed:

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)

For @Aaron who asked:

I am hoping to find a gdalwarp version of @wwnick's answer that utilizes the -multi option for enhanced multicore and multithreaded operations

Slight Disclaimer

This uses gdalwarp, though I'm not totally convinced there will be much performance gain. So far I've been I/O bound - running this script on a large raster cutting it into many smaller parts doesn't seem CPU intensive, so I assume the bottleneck is writing to disk. If you're planning on simultaneously re-projecting the tiles or something similar, then this might change. There are tuning tips here. A brief play didn't yield any improvement for me, and CPU never seemed to be the limiting factor.

Disclaimer aside, here's a script that will use gdalwarp to split up a raster into several smaller tiles. There might be some loss due to the floor division but this can be taken care of by picking the number of tiles you want. It will be n+1 where n is the number you divide by to get the tile_width and tile_height variables.

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

You can use r.tile of GRASS GIS. r.tile generates a separate raster map for each tile with numbered map names based on the user defined prefix. Width of tiles (columns) and height of tiles (rows) can be defined.

Using the grass-session Python API only a few lines of Python code are needed to call the r.tile functionality from outside, i.e. to write a standalone script. Using r.external and r.external.out also no data duplication occurs during the GRASS GIS processing step.

Pseudo code:

1. initialize grass-session
2. define output format with r.external.out
3. link input file with r.external
4. run r.tile which generates the tiles in the format defined above
5. unlink r.external.out
6. close grass-session

There's a bundled python script specifically for retiling rasters, 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

e.g:

gdal_retile.py -ps 512 512 -targetDir C:exampledir some_dem.tif

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)