POC splitting input AOI into within/intersecting tiles

aoi.geojson

intersect.geojson

tile_geometry.py

import fiona
import mercantile
from shapely.geometry import shape 


def process_tile(tile_list, aoi):
    # main function to compare a list of tiles to an input geometry
    # will eventually return a list of tiles completely within the aoi (all zoom levels possible)
    # and tiles that intersect the aoi (must be z12 because that's as low as we go)
    # we could keep subdividing to higher zooms, but our data is z12 and that's all we care about

    within_list = []
    intersect_list = []

    for t in tile_list:

        # a tile either is completely within, completely outside, or intersects
        tile_geom = shape(mercantile.feature(t)['geometry'])

        # if it's within, great- our work is done
        if aoi.contains(tile_geom):
            within_list.append(t)

        elif tile_geom.intersects(aoi):

            # if it intersects and is < z12, subdivide it and start the 
            # process again
            if t.z < 12:

                # find the four children of this tile and check them for within/intersect/outside-ness
                tile_children = mercantile.children(t)
                new_within_list, new_intersect_list = process_tile(tile_children, aoi)

                # add the results to our initial lists
                within_list.extend(new_within_list)
                intersect_list.extend(new_intersect_list)

            # if it intersects and is at z12, add it to our intersect list
            else:
                intersect_list.append(t)

        # and if it's outside our geometry, drop it
        else:
            pass

    return within_list, intersect_list


def write_tiles_to_geojson(tile_list, output_file):

    schema={'geometry': 'Polygon', 'properties': {'title': 'str'}}

    with fiona.open(output_file, 'w', driver='GeoJSON', schema=schema) as outfile:
        for t in tile_list:
            feat = mercantile.feature(t)
            outfile.write(feat)


def main():

    # read in aoi 
    src = fiona.open('aoi.geojson')
    geom = shape(src[0]['geometry'])

    # use bounds to find the smallest tile that completely contains our input aoi
    # not useful for AOIs that cross lat or lon 0 (returns tile [0, 0, 0])
    # but helpful for many AOIs
    # https://github.com/mapbox/mercantile/blob/master/docs/cli.rst#bounding-tile
    bbox = src.bounds
    bounding_tile = mercantile.bounding_tile(*bbox)

    # build within and intersect lists
    within_list, intersect_list = process_tile([bounding_tile], geom)

    # write to geojson
    write_tiles_to_geojson(within_list, 'within.geojson')
    write_tiles_to_geojson(intersect_list, 'intersect.geojson')


if __name__ == '__main__':
    main()

within.geojson