split.py

import sys

if (len(sys.argv) != 3):
    print('Usage: '+sys.argv[0]+' geotiff_image.tif block_size')
    sys.exit()

filename = sys.argv[1]
block_size = int(sys.argv[2])

from osgeo import gdal

driver = gdal.GetDriverByName("GTiff")
def save(block, name):
    outData = driver.Create(name, block.shape[1], block.shape[0])
    outData.GetRasterBand(1).WriteArray(block)

data = gdal.Open(filename)
# print(data.RasterCount, data.RasterXSize, data.RasterYSize)
from math import ceil
n = int(ceil(data.RasterYSize / block_size))
m = int(ceil(data.RasterXSize / block_size))
for i in range(n):
    for j in range(m):
        x = i * block_size
        y = j * block_size
        w = min(block_size, data.RasterXSize - x)
        h = min(block_size, data.RasterYSize - y)
        block = data.ReadAsArray(x, y, w, h)
        save(block, 'out'+'_'+str(j)+'_'+str(i)+'.tif')