2

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
2
  • 1
    Hi, I edited the question to include the full gdalinfo of the raster. Thanks..
    – CellDamag3
    Jul 12, 2020 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, 2020 at 2:47

1 Answer 1

2

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))
3
  • 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, 2020 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, 2020 at 22:30
  • Excellent. You've been invaluable. Much appreciation. Thanks.
    – CellDamag3
    Jul 13, 2020 at 22: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.