1 #pragma once
2
3 #include <mapbox/geojsonvt/convert.hpp>
4 #include <mapbox/geojsonvt/tile.hpp>
5 #include <mapbox/geojsonvt/types.hpp>
6 #include <mapbox/geojsonvt/wrap.hpp>
7
8 #include <chrono>
9 #include <cmath>
10 #include <map>
11 #include <unordered_map>
12
13 namespace mapbox {
14 namespace geojsonvt {
15
16 using geometry = mapbox::geometry::geometry<double>;
17 using feature = mapbox::geometry::feature<double>;
18 using feature_collection = mapbox::geometry::feature_collection<double>;
19 using geometry_collection = mapbox::geometry::geometry_collection<double>;
20 using geojson = mapbox::util::variant<geometry, feature, feature_collection>;
21
22 struct ToFeatureCollection {
operator ()mapbox::geojsonvt::ToFeatureCollection23 feature_collection operator()(const feature_collection& value) const {
24 return value;
25 }
operator ()mapbox::geojsonvt::ToFeatureCollection26 feature_collection operator()(const feature& value) const {
27 return { value };
28 }
operator ()mapbox::geojsonvt::ToFeatureCollection29 feature_collection operator()(const geometry& value) const {
30 return { { value } };
31 }
32 };
33
34 struct TileOptions {
35 // simplification tolerance (higher means simpler)
36 double tolerance = 3;
37
38 // tile extent
39 uint16_t extent = 4096;
40
41 // tile buffer on each side
42 uint16_t buffer = 64;
43 };
44
45 struct Options : TileOptions {
46 // max zoom to preserve detail on
47 uint8_t maxZoom = 18;
48
49 // max zoom in the tile index
50 uint8_t indexMaxZoom = 5;
51
52 // max number of points per tile in the tile index
53 uint32_t indexMaxPoints = 100000;
54 };
55
56 const Tile empty_tile{};
57
toID(uint8_t z,uint32_t x,uint32_t y)58 inline uint64_t toID(uint8_t z, uint32_t x, uint32_t y) {
59 return (((1ull << z) * y + x) * 32) + z;
60 }
61
geoJSONToTile(const geojson & geojson_,uint8_t z,uint32_t x,uint32_t y,const TileOptions & options=TileOptions (),bool wrap=false,bool clip=false)62 inline const Tile geoJSONToTile(const geojson& geojson_,
63 uint8_t z,
64 uint32_t x,
65 uint32_t y,
66 const TileOptions& options = TileOptions(),
67 bool wrap = false,
68 bool clip = false) {
69
70 const auto features_ = geojson::visit(geojson_, ToFeatureCollection{});
71 auto z2 = 1u << z;
72 auto tolerance = (options.tolerance / options.extent) / z2;
73 auto features = detail::convert(features_, tolerance);
74 if (wrap) {
75 features = detail::wrap(features, double(options.buffer) / options.extent);
76 }
77 if (clip) {
78 const double p = options.buffer / options.extent;
79
80 const auto left = detail::clip<0>(features, (x - p) / z2, (x + 1 + p) / z2, -1, 2);
81 features = detail::clip<1>(left, (y - p) / z2, (y + 1 + p) / z2, -1, 2);
82 }
83 return detail::InternalTile({ features, z, x, y, options.extent, tolerance }).tile;
84 }
85
86 class GeoJSONVT {
87 public:
88 const Options options;
89
GeoJSONVT(const mapbox::geometry::feature_collection<double> & features_,const Options & options_=Options ())90 GeoJSONVT(const mapbox::geometry::feature_collection<double>& features_,
91 const Options& options_ = Options())
92 : options(options_) {
93
94 const uint32_t z2 = 1u << options.maxZoom;
95
96 auto converted = detail::convert(features_, (options.tolerance / options.extent) / z2);
97 auto features = detail::wrap(converted, double(options.buffer) / options.extent);
98
99 splitTile(features, 0, 0, 0);
100 }
101
GeoJSONVT(const geojson & geojson_,const Options & options_=Options ())102 GeoJSONVT(const geojson& geojson_, const Options& options_ = Options())
103 : GeoJSONVT(geojson::visit(geojson_, ToFeatureCollection{}), options_) {
104 }
105
106 std::map<uint8_t, uint32_t> stats;
107 uint32_t total = 0;
108
getTile(const uint8_t z,const uint32_t x_,const uint32_t y)109 const Tile& getTile(const uint8_t z, const uint32_t x_, const uint32_t y) {
110
111 if (z > options.maxZoom)
112 throw std::runtime_error("Requested zoom higher than maxZoom: " + std::to_string(z));
113
114 const uint32_t z2 = 1u << z;
115 const uint32_t x = ((x_ % z2) + z2) % z2; // wrap tile x coordinate
116 const uint64_t id = toID(z, x, y);
117
118 auto it = tiles.find(id);
119 if (it != tiles.end())
120 return it->second.tile;
121
122 it = findParent(z, x, y);
123
124 if (it == tiles.end())
125 throw std::runtime_error("Parent tile not found");
126
127 // if we found a parent tile containing the original geometry, we can drill down from it
128 const auto& parent = it->second;
129
130 // drill down parent tile up to the requested one
131 splitTile(parent.source_features, parent.z, parent.x, parent.y, z, x, y);
132
133 it = tiles.find(id);
134 if (it != tiles.end())
135 return it->second.tile;
136
137 it = findParent(z, x, y);
138 if (it == tiles.end())
139 throw std::runtime_error("Parent tile not found");
140
141 return empty_tile;
142 }
143
getInternalTiles() const144 const std::unordered_map<uint64_t, detail::InternalTile>& getInternalTiles() const {
145 return tiles;
146 }
147
148 private:
149 std::unordered_map<uint64_t, detail::InternalTile> tiles;
150
151 std::unordered_map<uint64_t, detail::InternalTile>::iterator
findParent(const uint8_t z,const uint32_t x,const uint32_t y)152 findParent(const uint8_t z, const uint32_t x, const uint32_t y) {
153 uint8_t z0 = z;
154 uint32_t x0 = x;
155 uint32_t y0 = y;
156
157 const auto end = tiles.end();
158 auto parent = end;
159
160 while ((parent == end) && (z0 != 0)) {
161 z0--;
162 x0 = x0 / 2;
163 y0 = y0 / 2;
164 parent = tiles.find(toID(z0, x0, y0));
165 }
166
167 return parent;
168 }
169
splitTile(const detail::vt_features & features,const uint8_t z,const uint32_t x,const uint32_t y,const uint8_t cz=0,const uint32_t cx=0,const uint32_t cy=0)170 void splitTile(const detail::vt_features& features,
171 const uint8_t z,
172 const uint32_t x,
173 const uint32_t y,
174 const uint8_t cz = 0,
175 const uint32_t cx = 0,
176 const uint32_t cy = 0) {
177
178 const double z2 = 1u << z;
179 const uint64_t id = toID(z, x, y);
180
181 auto it = tiles.find(id);
182
183 if (it == tiles.end()) {
184 const double tolerance =
185 (z == options.maxZoom ? 0 : options.tolerance / (z2 * options.extent));
186
187 it = tiles
188 .emplace(id,
189 detail::InternalTile{ features, z, x, y, options.extent, tolerance })
190 .first;
191 stats[z] = (stats.count(z) ? stats[z] + 1 : 1);
192 total++;
193 // printf("tile z%i-%i-%i\n", z, x, y);
194 }
195
196 auto& tile = it->second;
197
198 if (features.empty())
199 return;
200
201 // if it's the first-pass tiling
202 if (cz == 0u) {
203 // stop tiling if we reached max zoom, or if the tile is too simple
204 if (z == options.indexMaxZoom || tile.tile.num_points <= options.indexMaxPoints) {
205 tile.source_features = features;
206 return;
207 }
208
209 } else { // drilldown to a specific tile;
210 // stop tiling if we reached base zoom
211 if (z == options.maxZoom)
212 return;
213
214 // stop tiling if it's our target tile zoom
215 if (z == cz) {
216 tile.source_features = features;
217 return;
218 }
219
220 // stop tiling if it's not an ancestor of the target tile
221 const double m = 1u << (cz - z);
222 if (x != static_cast<uint32_t>(std::floor(cx / m)) ||
223 y != static_cast<uint32_t>(std::floor(cy / m))) {
224 tile.source_features = features;
225 return;
226 }
227 }
228
229 const double p = 0.5 * options.buffer / options.extent;
230 const auto& min = tile.bbox.min;
231 const auto& max = tile.bbox.max;
232
233 const auto left = detail::clip<0>(features, (x - p) / z2, (x + 0.5 + p) / z2, min.x, max.x);
234
235 splitTile(detail::clip<1>(left, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
236 x * 2, y * 2, cz, cx, cy);
237 splitTile(detail::clip<1>(left, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
238 x * 2, y * 2 + 1, cz, cx, cy);
239
240 const auto right =
241 detail::clip<0>(features, (x + 0.5 - p) / z2, (x + 1 + p) / z2, min.x, max.x);
242
243 splitTile(detail::clip<1>(right, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
244 x * 2 + 1, y * 2, cz, cx, cy);
245 splitTile(detail::clip<1>(right, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
246 x * 2 + 1, y * 2 + 1, cz, cx, cy);
247
248 // if we sliced further down, no need to keep source geometry
249 tile.source_features = {};
250 }
251 };
252
253 } // namespace geojsonvt
254 } // namespace mapbox
255