-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain.py
95 lines (77 loc) · 3.27 KB
/
main.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
import os
import json
import argparse
import geopandas as gpd
import supermercado as sm
import mercantile
import smart_open
parser = argparse.ArgumentParser(description="Generate tiles covering roads and save as GeoJSON.")
parser.add_argument("shapefile_path", type=str, help="Path to the shapefile containing road data")
parser.add_argument("zoom_level", type=int, help="Zoom level for tile generation")
args = parser.parse_args()
shapefile_path = args.shapefile_path
zoom_level = args.zoom_level
# Load shapefile
gdf = gpd.read_file(shapefile_path)
gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]
# Convert to WGS84 (EPSG:4326) if needed
if gdf.crs is not None and gdf.crs.to_epsg() != 4326:
gdf = gdf.to_crs(epsg=4326)
# Generate tile coverage
def generate_tiles(geometry, zoom):
"""Generate tile indices using mercantile instead of burntiles."""
tiles = set()
for geom in geometry:
if geom.is_empty:
continue
bbox = geom.bounds # Get bounding box of geometry
tile_bounds = mercantile.tiles(bbox[0], bbox[1], bbox[2], bbox[3], zoom)
tiles.update(tile_bounds)
return tiles
# Generate tiles
tiles = generate_tiles(gdf.geometry, zoom_level)
output_dir = f"tiles_z{zoom_level}/"
os.makedirs(output_dir, exist_ok=True)
geojson_tiles = {"type": "FeatureCollection", "features": []}
tile_list_path = os.path.join(output_dir, f"tile_list_z{zoom_level}.txt")
geojson_path = os.path.join(output_dir, f"tiles_z{zoom_level}.geojson")
with open(tile_list_path, "w") as f:
for tile in tiles:
tile_bounds = mercantile.bounds(*tile)
feature = {
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[tile_bounds.west, tile_bounds.south],
[tile_bounds.east, tile_bounds.south],
[tile_bounds.east, tile_bounds.north],
[tile_bounds.west, tile_bounds.north],
[tile_bounds.west, tile_bounds.south],
]
],
},
"properties": {"z": int(tile[2]), "x": int(tile[0]), "y": int(tile[1])},
}
geojson_tiles["features"].append(feature)
f.write(f"{tile[2]}/{tile[0]}/{tile[1]}\n")
with open(geojson_path, "w") as f:
json.dump(geojson_tiles, f, indent=2)
# Upload to AWS S3 using smart_open
s3_tile_list = f"s3://planet.openhistoricalmap.org/tile_coverage/tiles_boundary_{zoom_level}.list"
s3_geojson = f"s3://planet.openhistoricalmap.org/tile_coverage/tiles_boundary_{zoom_level}.geojson"
print(f"Uploading {tile_list_path} to {s3_tile_list}")
with smart_open.open(s3_tile_list, "w") as s3_file:
with open(tile_list_path, "r") as local_file:
s3_file.write(local_file.read())
print(f"Uploading {geojson_path} to {s3_geojson}")
with smart_open.open(s3_geojson, "w") as s3_file:
with open(geojson_path, "r") as local_file:
s3_file.write(local_file.read())
print(
f"Tile list available at: https://s3.amazonaws.com/planet.openhistoricalmap.org/tile_coverage/tiles_boundary_{zoom_level}.list"
)
print(
f"GeoJSON file available at: https://s3.amazonaws.com/planet.openhistoricalmap.org/tile_coverage/tiles_boundary_{zoom_level}.geojson"
)