GeoJSON WebMercator Tile Covering

gistfile1.txt

import json
import tiletanic
from shapely import geometry, wkb, ops
from osgeo import ogr, osr
import pyproj
from functools import partial

# Cover a geometry with WebMercator tiles at zoom level 9.

# Should make a pull request to add something like this to the Tiletanic CLI

# Read in GeoJSON and get web mercator covering:
g = geometry.shape(json.load(open('aoi.geojson', 'r'))['features'][0]['geometry'])
og = ogr.Geometry(wkb=wkb.dumps(g))
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4326)
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(3857)
xform = osr.CoordinateTransformation(src_srs, dst_srs)
og.Transform(xform)
wmg = wkb.loads(og.ExportToWkb())
tiles = list(tiletanic.tilecover.cover_geometry(tiletanic.tileschemes.WebMercator(), wmg, 9))

# Create GeoJSON FeatureCollection of the covering:
projection = partial(pyproj.transform, pyproj.Proj(init='epsg:3857'), pyproj.Proj(init='epsg:4326'))
features = []
for t in tiles:
    b=scheme.bbox(t)
    # Dont forget to reproject, as GeoJSON should be in wgs84
    g = geometry.mapping(ops.transform(projection, geometry.box(b[0],b[1],b[2],b[3], False)))
    features.append( {"type": "Feature", "properties": {"quadkey": scheme.quadkey(t)}, "geometry": g})
json.dump({"type": "FeatureCollection", "features": features}, open("DYNAMIC_NorthAmerica_37_mercator_covering.geojson", "w"))