Fishnet split of geometry via shapely

split_geom.py

import os, sys
import argparse
from shapely.geometry import shape, mapping, box
import json

def split_geom(geometry, output_dir, base_name, cols, rows):
    bounds  = geometry.bounds
    xmin = bounds[0]
    xmax =  bounds[2]
    ymin = bounds[1]
    ymax = bounds[3]
    dx = (xmax - xmin) / cols
    dy = (ymax - ymin) / rows
    cut_geoms = {}
    for row in range(0, rows):
        for col in range(0, cols):
            b = box(xmin + (col * dx),
                    ymax - ((row + 1) * dy),
                    xmin + ((col + 1) * dx),
                    ymax - (row * dy))
            g = geometry.intersection(b)
            if not g.is_empty:
                cut_geoms[(col, row)] = g
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
    for (col, row), geom in cut_geoms.items():
        name = "{}_{}_{}".format(base_name, col, row)
        p = os.path.join(output_dir, "{}.geojson".format(name))
        fc = {}
        fc['type'] = 'FeatureCollection'
        fc['features'] = [{ 'type': 'Feature',
                            'properties': { 'name': name },
                            'geometry': mapping(geom) }]
        with open(p, 'w') as f:
            f.write(json.dumps(fc, indent=4))


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Split up an AOI.')
    parser.add_argument('cols', metavar='COLS', type=int,
                        help='Number of cols in the split grid.')
    parser.add_argument('rows', metavar='ROWS', type=int,
                        help='Number of rows in the split grid.')
    parser.add_argument('input_path', metavar='INPUT_PATH')
    parser.add_argument('output_dir', metavar='OUTPUT_DIR')

    args = parser.parse_args()

    base_name = os.path.splitext(os.path.basename(args.input_path))[0]

    gj = json.loads(open(args.input_path).read())

    features = gj['features']
    if not len(features) == 1:
        print('Feature collection must only contain one feature')
        sys.exit(1)

    geom = shape(features[0]['geometry'])

    split_geom(geom, args.output_dir, base_name, args.cols, args.rows)

This is awesome, was looking for GeoPandas replacement of Arcpy's Fishnet tool and this nails it.