2

I have been using this answer on GIS SE to split 3-band (RGB) imagery into 256x256 3-band image tiles:

import os, gdal

in_path = '/path/to/indata/'
input_filename = 'dtm_5.tif'

out_path = '/path/to/output_folder/'
output_filename = 'tile_'

tile_size_x = 256
tile_size_y = 256

ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize

for i in range(0, xsize, tile_size_x):
    for j in range(0, ysize, tile_size_y):
        com_string = "gdal_translate -of GTIFF -srcwin " + str(i)+ ", " + str(j) + ", " + str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(i) + "_" + str(j) + ".tif"
        os.system(com_string)

Is there a comparable Rasterio or numpy approach?

  • 1
    Just a suggestion but using gdal_translate with os.system will work much faster if you use subprocess.Popen(['gdal_translate','-of','GTIFF'.... you can generate a bunch of them and have them working simultaneously up to the bottleneck of HDD access. With trying to do this directly with GDAL you need to create a new dataset on each iteration, calculate the upper left and set the geotransform then read from source / write to target.. it's just easier to do it using gdal_translate. – Michael Stimson Jun 7 '18 at 5:04
2

Below is a simple example (rasterio 1.0.0 or later, won't work in 0.3.6). There might be better/simpler ways (and there is an easier way if your raster is internally tiled and the tile block sizes match your desired output tile size).

The rasterio docs have some examples of concurrent processing if you want to go down that road.

import os
from itertools import product
import rasterio as rio
from rasterio import windows

in_path = '/path/to/indata/'
input_filename = 'dtm_5.tif'

out_path = '/path/to/output_folder/'
output_filename = 'tile_{}-{}.tif'

def get_tiles(ds, width=256, height=256):
    nols, nrows = ds.meta['width'], ds.meta['height']
    offsets = product(range(0, nols, width), range(0, nrows, height))
    big_window = windows.Window(col_off=0, row_off=0, width=nols, height=nrows)
    for col_off, row_off in  offsets:
        window =windows.Window(col_off=col_off, row_off=row_off, width=width, height=height).intersection(big_window)
        transform = windows.transform(window, ds.transform)
        yield window, transform


with rio.open(os.path.join(in_path, input_filename)) as inds:
    tile_width, tile_height = 256, 256

    meta = inds.meta.copy()

    for window, transform in get_tiles(inds):
        print(window)
        meta['transform'] = transform
        meta['width'], meta['height'] = window.width, window.height
        outpath = os.path.join(out_path,output_filename.format(int(window.col_off), int(window.row_off)))
        with rio.open(outpath, 'w', **meta) as outds:
            outds.write(inds.read(window=window))
  • Dang. I get AttributeError: module 'rasterio.windows' has no attribute 'Window'. – François Leblanc Oct 30 '18 at 15:25
  • Turns out there was a recent API break in rasterio from version 0.36.0 to 1.0.0. As of writing, if you're using Anaconda make sure you install from conda-forge. – François Leblanc Oct 30 '18 at 15:40
  • 1
    @FrancoisLeblanc My answer specifically notes in the first sentence that it is for rasterio 1.0 or later not 0.3.6. – user2856 Oct 30 '18 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.