Cut a multispectral TIFF into smaller parts

Question

I have a multispectral image with 6 bands that I want to cut into smaller parts. It seems like a trivial task but I can't make it work.

Approach 1: using rasterio.windows. I choose how many tiles I want, calculate how many pixels the individual tiles have and then read the corresponding part of the original image and write it to a new file:

with rasterio.open('path/raster.tif') as file:
    src = file.read()

bands, height, width = src.shape

grid_size = 5
grid_width = int(width/grid_size)
grid_height = int(height/grid_size)


with rasterio.open('path/raster.tif') as file:
    profile = file.profile
    profile.update(
    width=grid_width,
    height=grid_height)

    for i in range(grid_size):
        for j in range(grid_size):
            window = int(width*i/grid_size), int(height*j/grid_size), grid_width, grid_height
            
            src = file.read(window=Window(window[0], window[1], window[2], window[3]))
            
            with rasterio.open('data/cut_{}_{}.tif'.format(i,j), 'w', **profile) as sink:
                sink.write(src)

Problem: This correctly cuts the raster into pieces (in this code it is 5*5 = 25 tiles) but I don't know how to correctly set their coordinates (the "transform" parameter in profile). In my code, all the tiles are placed in the upper left corner of the original TIFF which makes sense because I use the original TIFF's profile for writing the tiles.

Approach 2: using gdal2tiles

import gdal2tiles

gdal2tiles.generate_tiles('data/clipped_data/kuopio.tif', 'data/')

Problem: I get an error: ERROR 6: PNG driver doesn't support 6 bands . So it seems this is made only for standard 3 band images. I haven't found an option to change that.

Question: How do I correctly cut a multispectral raster into a number of smaller rasters? I am using Python so this answer is not helpful for me.

Answer 1

You can set the output transform using the source dataset window_transform method:

etc...
window = int(width*i/grid_size), int(height*j/grid_size), grid_width, grid_height
transform = file.window_transform(window)
profile.update({'transform': transform})

Answer 2

You can use https://gdal.org/programs/gdal_translate.html. Perhaps it is easies to use option

-srcwin <xoff> <yoff> <xsize> <ysize>

Selects a subwindow from the source image for copying based on pixel/line location.

You can write the first 1000 by 1000 pixel part into GeoTIFF with a simple command

gdal_translate -of GTiff -srcwin 0 0 1000 1000 input.tif output_tile_1.tif

You can call gdal_translate executable with a shell script or from Python or you can use plain Python and utilise gdal.Translate. You can find a code example from https://github.com/OSGeo/gdal/blob/master/autotest/utilities/test_gdal_translate_lib.py#L209

Not so much code is needed

 ds = gdal.Open('../gcore/data/byte.tif')
 ds = gdal.Translate('tmp/test10.tif', ds, srcWin=[0, 0, 1, 1])