-
Notifications
You must be signed in to change notification settings - Fork 0
/
tilepyramidcut.py
193 lines (153 loc) · 6.2 KB
/
tilepyramidcut.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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
#!/usr/bin/env python3
"""
Prunes and copies portions of a tile pyramid such as one used in Google Earth or Mapbox.
The boundary polygon can be defined in DeepZoom (www.deepzoom.com) as one or more routes.
The routes will be converted to shapely Polygons and used as bondaries for the prune or
copy operation.
"""
__author__ = "Jay Borseth"
__version__ = "0.1.0"
__license__ = "MIT"
import argparse
import json
import os
import shutil
import numpy as np
from logzero import logger
from pyproj import Transformer
from shapely import from_geojson
from shapely.geometry import MultiPolygon, Polygon, shape
from tilematrix import Tile, TilePyramid
def dir_is_empty(path):
""" returns True if diretory exists and is empty, False if not empty, and None if doesn't exist """
if os.path.isdir(path):
with os.scandir(path) as it:
if any(it):
return False
return True
return None
def dir_delete_if_empty(args, tp: TilePyramid, tile: Tile, path: str):
""" Delete the directory if empty, and walk up the pyramid recursively deleting any empty ancestor directories.
return True if directory was deleted """
if os.path.isdir(path):
with os.scandir(path) as it:
if not any(it):
# dir is empty
os.rmdir(path)
if args.verbose >= 1:
logger.info(f"rmdir: {path}")
# recursively walk up the pyramid checking if parent is empty
parent = tile.get_parent()
if parent and parent.zoom >= args.zbegin:
parent_path = tuple_to_path(args.src, parent)
parent_dir = os.path.dirname(parent_path)
while dir_delete_if_empty(args, tp, parent, parent_dir):
return True
return False
def tuple_to_path (root, zyx):
""" Given a root directory and a zyx tuple, create the .png path """
z, y, x = zyx
full = os.path.join(root, f"Z{z}/{y}/{x}.png")
return full
def geojson_to_multipolygon (path):
""" Given a JSON file containing containing LatLng LINESTIRNGS OR POLYGONS,
as either a DeepZoom feature list or else a normal geojson feature list,
convert it to a shapely multipolygon in ESPG:3857 (UTM meters) units."""
with open(path, "r") as read_file:
geo = json.load(read_file)
if "featureCollection" in geo:
# geojson is embedded in a DeepZoom trip
featureCollection = geo["featureCollection"]
else:
# regular geoson
featureCollection = geo
features = featureCollection["features"]
t4326 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
allPolygons = []
for feature in features:
xx = list(map(lambda p:p[0], feature["geometry"]["coordinates"]))
yy = list(map(lambda p:p[1], feature["geometry"]["coordinates"]))
# following not needed: shapely autoconverts LINESTRING to POLYGON
# xx.append(xx[0])
# yy.append(yy[0])
mxx, myy = t4326.transform(xx, yy)
meters = [(i, j) for i, j in zip(mxx, myy)]
# print (meters)
polym = Polygon(meters)
allPolygons.append(polym)
multiPolygon = MultiPolygon(allPolygons)
return multiPolygon
def test (args):
""" LatLng to UTM and back """
tp = TilePyramid("mercator")
lon = 95.99
lat = 15.15
t4326 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
utm = t4326.transform(lon, lat)
print (utm)
t3857 = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)
lonlat = t3857.transform(utm[0], utm[1])
print (lonlat)
def op(args):
""" The point of it all """
tp = TilePyramid("mercator")
multiPolygon = geojson_to_multipolygon(args.polyfile)
count = 0
for zoom in range(args.zbegin, args.zend+1):
logger.info (f"Zoom: {zoom}")
tiles = tp.tiles_from_geom(multiPolygon, zoom)
for tile in tiles:
src_path = tuple_to_path(args.src, tile)
if os.path.isfile(src_path):
count = count + 1
if args.op == "prune":
if (args.verbose > 1):
print ("prune: ", src_path)
os.remove(src_path)
dir_path = os.path.dirname(src_path)
dir_delete_if_empty(args, tp, tile, dir_path)
elif args.op == "copy":
if (args.verbose > 1):
print ("copy: ", src_path)
dst_path = tuple_to_path(args.dst, tile)
dst_dir = os.path.dirname(dst_path)
os.makedirs(dst_dir, exist_ok = True)
shutil.copy2(src_path, dst_path)
else:
if (args.verbose > 1):
print ("missing: ", src_path)
logger.info(f"TOTAL {args.op}: {count}")
def main(args):
logger.info("tilepyramidcut")
if not os.path.isdir(args.src):
logger.error (f"Source directory (--src) {args.src} does not exist")
if not os.path.isfile(args.polyfile):
logger.error (f"Polygon file (--polyfile) {args.polyfile} does not exist")
# "prune" doesn't use a dst directory, but copy requires one
if args.op == "copy" and not os.path.isdir(args.dst):
logger.error (f"Destination directory (--dst) {args.src} does not exist")
logger.info(args)
op(args)
if __name__ == "__main__":
parser = argparse.ArgumentParser()
# Required positional argument
parser.add_argument("--src", help="source directory to process", dest="src", required=True)
parser.add_argument("--dst", action="store", dest="dst", help="destination directory for copy op")
parser.add_argument("--polyfile", action="store", dest="polyfile", required=True, type=str, help="json file containing geojson geometry of Polygon(s) to prune or copy")
parser.add_argument("--op", action="store", choices=["prune", "copy"], required=True, help="the operation to perform")
parser.add_argument("--zbegin", action="store", dest="zbegin", required=True, type=int, help= "start zoom")
parser.add_argument("--zend", action="store", dest="zend", required=True, type=int, help="end zoom")
# Optional verbosity counter (eg. -v, -vv, -vvv, etc.)
parser.add_argument(
"-v",
"--verbose",
action="count",
default=0,
help="Verbosity (-v, -vv, etc)")
# Specify output of "--version"
parser.add_argument(
"--version",
action="version",
version="%(prog)s (version {version})".format(version=__version__))
args = parser.parse_args()
main(args)