-
Notifications
You must be signed in to change notification settings - Fork 47
/
encoder.py
440 lines (348 loc) · 14.5 KB
/
encoder.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
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
from math import fabs
from past.builtins import long, unicode
from numbers import Number
from shapely.geometry.base import BaseGeometry
from shapely.geometry.multipolygon import MultiPolygon
from shapely.geometry.polygon import orient, Polygon
from shapely.wkb import loads as load_wkb
from shapely.wkt import loads as load_wkt
import sys
PY3 = sys.version_info[0] == 3
if PY3:
from .Mapbox import vector_tile_pb2_p3 as vector_tile
xrange = range
def apply_map(fn, x):
return list(map(fn, x))
else:
from .Mapbox import vector_tile_pb2 as vector_tile
def apply_map(fn, x):
return map(fn, x)
# tiles are padded by this number of pixels for the current zoom level
padding = 0
cmd_bits = 3
tolerance = 0
CMD_MOVE_TO = 1
CMD_LINE_TO = 2
CMD_SEG_END = 7
def transform(shape, func):
''' Ported from TileStache'''
construct = shape.__class__
if shape.type.startswith('Multi'):
parts = [transform(geom, func) for geom in shape.geoms]
return construct(parts)
if shape.type in ('Point', 'LineString'):
return construct(apply_map(func, shape.coords))
if shape.type == 'Polygon':
exterior = apply_map(func, shape.exterior.coords)
rings = [apply_map(func, ring.coords) for ring in shape.interiors]
return construct(exterior, rings)
if shape.type == 'GeometryCollection':
return construct()
raise ValueError('Unknown geometry type, "%s"' % shape.type)
class VectorTile:
"""
"""
def __init__(self, extents):
self.tile = vector_tile.tile()
self.extents = extents
def addFeatures(self, features, layer_name='',
quantize_bounds=None, y_coord_down=False):
self.layer = self.tile.layers.add()
self.layer.name = layer_name
self.layer.version = 2
self.layer.extent = self.extents
self.key_idx = 0
self.val_idx = 0
self.seen_keys_idx = {}
self.seen_values_idx = {}
for feature in features:
# skip missing or empty geometries
geometry_spec = feature.get('geometry')
if geometry_spec is None:
continue
shape = self._load_geometry(geometry_spec)
if shape is None:
raise NotImplementedError(
'Can\'t do geometries that are not wkt, wkb, or shapely '
'geometries')
if shape.is_empty:
continue
if quantize_bounds:
shape = self.quantize(shape, quantize_bounds)
if shape.type == 'MultiPolygon':
# If we are a multipolygon, we need to ensure that the
# winding orders of the consituent polygons are
# correct. In particular, the winding order of the
# interior rings need to be the opposite of the
# exterior ones, and all interior rings need to follow
# the exterior one. This is how the end of one polygon
# and the beginning of another are signaled.
shape = self.enforce_multipolygon_winding_order(shape)
elif shape.type == 'Polygon':
# Ensure that polygons are also oriented with the
# appropriate winding order. Their exterior rings must
# have a clockwise order, which is translated into a
# clockwise order in MVT's tile-local coordinates with
# the Y axis in "screen" (i.e: +ve down) configuration.
# Note that while the Y axis flips, we also invert the
# Y coordinate to get the tile-local value, which means
# the clockwise orientation is unchanged.
shape = self.enforce_polygon_winding_order(shape)
self.addFeature(feature, shape, y_coord_down)
def quantize(self, shape, bounds):
minx, miny, maxx, maxy = bounds
def fn(point):
x, y = point
xfac = self.extents / (maxx - minx)
yfac = self.extents / (maxy - miny)
x = xfac * (x - minx)
y = yfac * (y - miny)
return round(x), round(y)
return transform(shape, fn)
def enforce_multipolygon_winding_order(self, shape):
assert shape.type == 'MultiPolygon'
parts = []
for part in shape.geoms:
# see comment in shape.type == 'Polygon' above about why
# the sign here has to be -1.
part = self.enforce_polygon_winding_order(part)
parts.append(part)
oriented_shape = MultiPolygon(parts)
return oriented_shape
def enforce_polygon_winding_order(self, shape):
assert shape.type == 'Polygon'
def fn(point):
x, y = point
return round(x), round(y)
exterior = apply_map(fn, shape.exterior.coords)
rings = None
if len(shape.interiors) > 0:
rings = [apply_map(fn, ring.coords) for ring in shape.interiors]
return orient(Polygon(exterior, rings), sign=-1.0)
def _load_geometry(self, geometry_spec):
if isinstance(geometry_spec, BaseGeometry):
return geometry_spec
try:
return load_wkb(geometry_spec)
except:
try:
return load_wkt(geometry_spec)
except:
return None
def addFeature(self, feature, shape, y_coord_down):
f = self.layer.features.add()
fid = feature.get('id')
if fid is not None:
if isinstance(fid, Number) and fid >= 0:
f.id = fid
# properties
properties = feature.get('properties')
if properties is not None:
self._handle_attr(self.layer, f, properties)
f.type = self._get_feature_type(shape)
self._geo_encode(f, shape, y_coord_down)
def _get_feature_type(self, shape):
if shape.type == 'Point' or shape.type == 'MultiPoint':
return self.tile.Point
elif shape.type == 'LineString' or shape.type == 'MultiLineString':
return self.tile.LineString
elif shape.type == 'Polygon' or shape.type == 'MultiPolygon':
return self.tile.Polygon
def _encode_cmd_length(self, cmd, length):
return (length << cmd_bits) | (cmd & ((1 << cmd_bits) - 1))
def _chunker(self, seq, size):
return [seq[pos:pos + size] for pos in xrange(0, len(seq), size)]
def _can_handle_key(self, k):
return isinstance(k, str) or \
isinstance(k, str if PY3 else unicode)
def _can_handle_val(self, v):
if isinstance(v, str) or \
isinstance(v, str if PY3 else unicode):
return True
elif isinstance(v, bool):
return True
elif (isinstance(v, int) or
isinstance(v, int if PY3 else long)):
return True
elif isinstance(v, float):
return True
return False
def _can_handle_attr(self, k, v):
return self._can_handle_key(k) and \
self._can_handle_val(v)
def _handle_attr(self, layer, feature, props):
for k, v in props.items():
if self._can_handle_attr(k, v):
if not PY3 and isinstance(k, str):
k = k.decode('utf-8')
if k not in self.seen_keys_idx:
layer.keys.append(k)
self.seen_keys_idx[k] = self.key_idx
self.key_idx += 1
feature.tags.append(self.seen_keys_idx[k])
if v not in self.seen_values_idx:
self.seen_values_idx[v] = self.val_idx
self.val_idx += 1
val = layer.values.add()
if isinstance(v, bool):
val.bool_value = v
elif (isinstance(v, str)):
if PY3:
val.string_value = str(v)
else:
val.string_value = unicode(v, 'utf8')
elif (isinstance(v, str if PY3 else unicode)):
val.string_value = v
elif (isinstance(v, int)) or (
isinstance(v, int if PY3 else long)):
val.int_value = v
elif (isinstance(v, float)):
val.double_value = v
feature.tags.append(self.seen_values_idx[v])
def _handle_skipped_last(self, f, skipped_index, cur_x, cur_y, x_, y_):
last_x = f.geometry[skipped_index - 2]
last_y = f.geometry[skipped_index - 1]
last_dx = ((last_x >> 1) ^ (-(last_x & 1)))
last_dy = ((last_y >> 1) ^ (-(last_y & 1)))
dx = cur_x - x_ + last_dx
dy = cur_y - y_ + last_dy
x_ = cur_x
y_ = cur_y
f.geometry.__setitem__(skipped_index - 2, ((dx << 1) ^ (dx >> 31)))
f.geometry.__setitem__(skipped_index - 1, ((dy << 1) ^ (dy >> 31)))
def _parseGeometry(self, shape):
coordinates = []
lineType = "line"
polygonType = "polygon"
def _get_point_obj(x, y, cmd=CMD_MOVE_TO):
coordinate = {
'x': x,
'y': y,
'cmd': cmd
}
coordinates.append(coordinate)
def _get_arc_obj(arc, type):
length = len(arc.coords)
iterator = 0
cmd = CMD_MOVE_TO
while (iterator < length):
x = arc.coords[iterator][0]
y = arc.coords[iterator][1]
if iterator == 0:
cmd = CMD_MOVE_TO
elif iterator == length - 1 and type == polygonType:
cmd = CMD_SEG_END
else:
cmd = CMD_LINE_TO
_get_point_obj(x, y, cmd)
iterator = iterator + 1
if shape.type == 'GeometryCollection':
# do nothing
coordinates = []
elif shape.type == 'Point':
_get_point_obj(shape.x, shape.y)
elif shape.type == 'LineString':
_get_arc_obj(shape, lineType)
elif shape.type == 'Polygon':
rings = [shape.exterior] + list(shape.interiors)
for ring in rings:
_get_arc_obj(ring, polygonType)
elif shape.type == 'MultiPoint':
for point in shape.geoms:
_get_point_obj(point.x, point.y)
elif shape.type == 'MultiLineString':
for arc in shape.geoms:
_get_arc_obj(arc, lineType)
elif shape.type == 'MultiPolygon':
for polygon in shape.geoms:
rings = [polygon.exterior] + list(polygon.interiors)
for ring in rings:
_get_arc_obj(ring, polygonType)
else:
raise NotImplementedError("Can't do %s geometries" % shape.type)
return coordinates
def _geo_encode(self, f, shape, y_coord_down):
x_, y_ = 0, 0
cmd = -1
cmd_idx = -1
vtx_cmd = -1
prev_cmd = -1
skipped_index = -1
skipped_last = False
cur_x = 0
cur_y = 0
it = 0
length = 0
coordinates = self._parseGeometry(shape)
while (True):
if it >= len(coordinates):
break
c_it = coordinates[it]
x, y, vtx_cmd = c_it['x'], c_it['y'], c_it['cmd']
if vtx_cmd != cmd:
if cmd_idx >= 0:
f.geometry[cmd_idx] = self._encode_cmd_length(cmd, length)
cmd = vtx_cmd
length = 0
cmd_idx = len(f.geometry)
f.geometry.append(0) # placeholder added in first pass
if (vtx_cmd == CMD_MOVE_TO or vtx_cmd == CMD_LINE_TO):
if cmd == CMD_MOVE_TO and skipped_last and skipped_index > 1:
self._handle_skipped_last(
f, skipped_index, cur_x, cur_y, x_, y_)
# ensure that floating point values don't get truncated
if isinstance(x, float):
x = round(x)
if isinstance(y, float):
y = round(y)
x = int(x)
y = int(y)
if not y_coord_down:
y = self.extents - y
# Compute delta to the previous coordinate.
cur_x = int(x)
cur_y = int(y)
dx = cur_x - x_
dy = cur_y - y_
sharp_turn_ahead = False
if (it + 2 <= len(coordinates)):
next_coord = coordinates[it + 1]
if next_coord['cmd'] == CMD_LINE_TO:
next_x, next_y = next_coord['x'], next_coord['y']
next_dx = fabs(cur_x - int(next_x))
next_dy = fabs(cur_y - int(next_y))
if ((next_dx == 0 and next_dy >= tolerance) or
(next_dy == 0 and next_dx >= tolerance)):
sharp_turn_ahead = True
# Keep all move_to commands, but omit other movements
# that are not >= the tolerance threshold and should
# be considered no-ops.
# NOTE: length == 0 indicates the command has changed and will
# preserve any non duplicate move_to or line_to
if (length == 0 or sharp_turn_ahead or
fabs(dx) >= tolerance or fabs(dy) >= tolerance):
# Manual zigzag encoding.
f.geometry.append((dx << 1) ^ (dx >> 31))
f.geometry.append((dy << 1) ^ (dy >> 31))
x_ = cur_x
y_ = cur_y
skipped_last = False
length = length + 1
else:
skipped_last = True
skipped_index = len(f.geometry)
elif vtx_cmd == CMD_SEG_END:
if prev_cmd != CMD_SEG_END:
length = length + 1
else:
raise Exception("Unknown command type: '%s'" % vtx_cmd)
it = it + 1
prev_cmd = cmd
# at least one vertex + cmd/length
if skipped_last and skipped_index > 1:
# if we skipped previous vertex we just update it to the
# last one here.
self._handle_skipped_last(f, skipped_index, cur_x, cur_y, x_, y_)
# Update the last length/command value.
if cmd_idx >= 0:
f.geometry[cmd_idx] = self._encode_cmd_length(cmd, length)