Created
June 3, 2018 23:35
-
-
Save mappingvermont/c9a0008ae91e19f2a5631014b8f52508 to your computer and use it in GitHub Desktop.
POC splitting input AOI into within/intersecting tiles
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"type": "FeatureCollection", | |
"features": [ | |
{ | |
"type": "Feature", | |
"properties": {}, | |
"geometry": { | |
"type": "Polygon", | |
"coordinates": [ | |
[ | |
[ | |
-76.2890625, | |
-2.460181181020993 | |
], | |
[ | |
-62.57812500000001, | |
-15.28418511407642 | |
], | |
[ | |
-62.57812500000001, | |
-23.563987128451217 | |
], | |
[ | |
-47.109375, | |
-20.3034175184893 | |
], | |
[ | |
-48.515625, | |
-6.315298538330033 | |
], | |
[ | |
-60.1171875, | |
5.266007882805498 | |
], | |
[ | |
-73.125, | |
2.8113711933311403 | |
], | |
[ | |
-76.2890625, | |
-2.460181181020993 | |
] | |
] | |
] | |
} | |
} | |
] | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment