chunk_contours.py

"""Splits rasters into chunks and creates contours"""

import os
import argparse
import uuid
import re
from osgeo import gdal

def chunk():
    """Splits rasters into 5000 x 5000 chunks."""
    parser = argparse.ArgumentParser()
    parser.add_argument("input")
    parser.add_argument("-s", type=int, help="tile size")
    parser.add_argument("-x", type=int, help="starting width")
    parser.add_argument("-y", type=int, help="starting height")
    parser.add_argument("-c", type=str, help="contour distance")
    args = parser.parse_args()

    dset = gdal.Open(args.input)

    width = dset.RasterXSize
    height = dset.RasterYSize

    print width, 'x', height
    uid = uuid.uuid4().hex
    
    if args.s:
        tilesize = args.s
    else:
        tilesize = 5000

    if args.x:
        _w1 = int(round(args.x / tilesize) * tilesize)
    else:
        _w1 = 0

    if args.y:
        _h1 = int(round(args.y / tilesize) * tilesize)
    else:
        _h1 = 0
        
    if args.c:
        contour = args.c
    else:
        contour = '10'

    for i in range(_w1, width, tilesize):
        for j in range(_h1, height, tilesize):
            _w = min(i + tilesize, width) - i
            _h = min(j + tilesize, height) - j
            grid = uid + "_" + str(i) + "_" + str(j) + ".tif"
            shp = re.sub(r'tif$', 'shp', grid, 0, re.MULTILINE)
            gdaltran_string = "gdal_translate -of GTIFF -srcwin " + str(i) + ", " + str(j)\
              + ", " + str(_w) + ", " + str(_h) + " " + args.input + " " + grid
            os.system(gdaltran_string)

            contour_string = 'gdal_contour -a elev ' + grid + ' ' + shp + ' -i ' + contour
            os.system(contour_string)
            os.remove(grid)

chunk()