Splitting a GeoTIFF raster file into tiles in meters

Question

I am new to GIS in general and I have started working with GDAL and Rasterio. What I am looking for is a method in GDAL or Rasterio that mimics the ArcGIS's SplitRaster function when the parameter tile_size is in meters (i.e. "unit" flag set to metres). I have seen multiple similar questions asked but I can only see solutions that would split/tile a large GeoTIFF by pixels.

In summary,

My input: A large GeoTIFF file (Planetscope)

Expected output: Bunch of GeoTIFF files (preserving geo-references) that are tiled of size 100m x 100m.

I have tried gdal_translate and techniques described in this question but I need tiling dimensions to be in meters not pixels.

Edit (Added full gdalinfo of the raster below):

Driver: GTiff/GeoTIFF
Files: Chicago.tif
Size is 10134, 9211
Coordinate System is:
PROJCS["NAD83 / UTM zone 16N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","26916"]]
Origin = (433874.999999512860086,4646969.999883539974689)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  433875.000, 4646970.000) ( 87d47'53.12"W, 41d58'19.02"N)
Lower Left  (  433875.000, 4619337.000) ( 87d47'42.00"W, 41d43'23.11"N)
Upper Right (  464277.000, 4646970.000) ( 87d25'52.21"W, 41d58'26.10"N)
Lower Right (  464277.000, 4619337.000) ( 87d25'46.20"W, 41d43'30.13"N)
Center      (  449076.000, 4633153.500) ( 87d36'48.39"W, 41d50'55.12"N)
Band 1 Block=10134x1 Type=Float32, ColorInterp=Gray

Hi, I edited the question to include the full gdalinfo of the raster. Thanks.. — CellDamag3, Jul 12 '20 at 22:03
How would you want to deal with the fact the 100m is not evenly divisible by your pixel size - closest (99m tile), closest that's at least 100m (102m tile), something else? — user2856, Jul 13 '20 at 02:47

score 2 · Accepted Answer · answered Jul 13 '20 at 05:49

Pretty much the same answer I gave in the linked question, the only difference is the width and height of the tiles may be specified in the units of the raster (metres in your case), then adjusted to convert them to the closest number of columns and rows by dividing by the desired height and width by pixel size.

import os.path
from itertools import product
import rasterio as rio
from rasterio import windows
infile = '/tmp/test.tif'
out_path = '/tmp/tiles'
output_filename = 'tile_{}-{}.tif'
def get_tiles(ds, width=256, height=256, map_units=False):
if map_units:
    # Get pixel size
    px, py = ds.transform.a, -ds.transform.e
    width, height = int(width / px + 0.5) , int(height / px + 0.5)

ncols, nrows = ds.meta['width'], ds.meta['height']

offsets = product(range(0, ncols, width), range(0, nrows, height))
big_window = windows.Window(col_off=0, row_off=0, width=ncols, 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(infile) as inds:
    tile_width, tile_height = 100, 100
meta = inds.meta.copy()

for window, transform in get_tiles(inds, tile_width, tile_height, map_units=True):

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

This is very useful. Really appreciate it. To answer your question, any rounding or clipping is fine for my use case in the event the corner pixels are not divisible by the pixel count. If you have any suggestions of throwing in some kind of checks for that (in the function you provided), that would be extremely valuable for me.
One last clarification, in def get_tiles(ds, width=256, height=256, map_units=False) you did not mean to hardcode width and height to 256 if I understood correctly. Again, a ton of gratitude. — CellDamag3, Jul 13 '20 at 21:26
Yes, the width=256, height=256 defaults are intentional (as that's the tile size I nearly always use). They're not hardcoded, they're defaults and can be overridden. So if I called get_tiles(ds) I would get 256x256 pixel tiles, but if I call get_tiles(ds, 100, 100, True) I would get 100x100 metre tiles (assuming the raster CRS uses metres). Re. your case of 3m pixel resolution and 100m tiles, you will get 99m tiles with this code as 100 is not evenly divisible by 3 and 99 is the closest evenly divisible number to 100. — user2856, Jul 13 '20 at 22:30
Excellent. You've been invaluable. Much appreciation. Thanks. — CellDamag3, Jul 13 '20 at 22:35

Splitting a GeoTIFF raster file into tiles in meters

1 Answers1