diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ff7064c3..6e226650a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +# 2.29.0 + +* Add tippecanoe-overzoom tool + # 2.28.1 * Allow --set-attribute to override an existing attribute value diff --git a/Makefile b/Makefile index 9791de1d8..e98524bd8 100644 --- a/Makefile +++ b/Makefile @@ -29,17 +29,18 @@ else FINAL_FLAGS := -g $(WARNING_FLAGS) $(DEBUG_FLAGS) endif -all: tippecanoe tippecanoe-enumerate tippecanoe-decode tile-join unit tippecanoe-json-tool +all: tippecanoe tippecanoe-enumerate tippecanoe-decode tile-join unit tippecanoe-json-tool tippecanoe-overzoom docs: man/tippecanoe.1 -install: tippecanoe tippecanoe-enumerate tippecanoe-decode tile-join tippecanoe-json-tool +install: tippecanoe tippecanoe-enumerate tippecanoe-decode tile-join tippecanoe-json-tool tippecanoe-overzoom mkdir -p $(PREFIX)/bin mkdir -p $(MANDIR) cp tippecanoe $(PREFIX)/bin/tippecanoe cp tippecanoe-enumerate $(PREFIX)/bin/tippecanoe-enumerate cp tippecanoe-decode $(PREFIX)/bin/tippecanoe-decode cp tippecanoe-json-tool $(PREFIX)/bin/tippecanoe-json-tool + cp tippecanoe-overzoom $(PREFIX)/bin/tippecanoe-overzoom cp tile-join $(PREFIX)/bin/tile-join cp man/tippecanoe.1 $(MANDIR)/tippecanoe.1 @@ -57,7 +58,7 @@ C = $(wildcard *.c) $(wildcard *.cpp) INCLUDES = -I/usr/local/include -I. LIBS = -L/usr/local/lib -tippecanoe: geojson.o jsonpull/jsonpull.o tile.o pool.o mbtiles.o geometry.o projection.o memfile.o mvt.o serial.o main.o text.o dirtiles.o pmtiles_file.o plugin.o read_json.o write_json.o geobuf.o flatgeobuf.o evaluator.o geocsv.o csv.o geojson-loop.o json_logger.o visvalingam.o compression.o +tippecanoe: geojson.o jsonpull/jsonpull.o tile.o pool.o mbtiles.o geometry.o projection.o memfile.o mvt.o serial.o main.o text.o dirtiles.o pmtiles_file.o plugin.o read_json.o write_json.o geobuf.o flatgeobuf.o evaluator.o geocsv.o csv.o geojson-loop.o json_logger.o visvalingam.o compression.o clip.o $(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread tippecanoe-enumerate: enumerate.o @@ -75,6 +76,9 @@ tippecanoe-json-tool: jsontool.o jsonpull/jsonpull.o csv.o text.o geojson-loop.o unit: unit.o text.o $(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread +tippecanoe-overzoom: overzoom.o mvt.o clip.o + $(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread + -include $(wildcard *.d) %.o: %.c @@ -87,12 +91,12 @@ clean: rm -f ./tippecanoe ./tippecanoe-* ./tile-join ./unit *.o *.d */*.o */*.d tests/**/*.mbtiles tests/**/*.check indent: - clang-format -i -style="{BasedOnStyle: Google, IndentWidth: 8, UseTab: Always, AllowShortIfStatementsOnASingleLine: false, ColumnLimit: 0, ContinuationIndentWidth: 8, SpaceAfterCStyleCast: true, IndentCaseLabels: false, AllowShortBlocksOnASingleLine: false, AllowShortFunctionsOnASingleLine: false, SortIncludes: false}" $(C) $(H) + clang-format -i -style="{BasedOnStyle: Google, IndentWidth: 8, UseTab: Always, AllowShortIfStatementsOnASingleLine: false, ColumnLimit: 0, ContinuationIndentWidth: 8, SpaceAfterCStyleCast: true, IndentCaseLabels: false, AllowShortBlocksOnASingleLine: false, AllowShortFunctionsOnASingleLine: false, SortIncludes: false}" $(filter-out flatgeobuf.cpp,$(C)) $(H) TESTS = $(wildcard tests/*/out/*.json) SPACE = $(NULL) $(NULL) -test: tippecanoe tippecanoe-decode $(addsuffix .check,$(TESTS)) raw-tiles-test parallel-test pbf-test join-test enumerate-test decode-test join-filter-test unit json-tool-test allow-existing-test csv-test layer-json-test pmtiles-test decode-pmtiles-test +test: tippecanoe tippecanoe-decode $(addsuffix .check,$(TESTS)) raw-tiles-test parallel-test pbf-test join-test enumerate-test decode-test join-filter-test unit json-tool-test allow-existing-test csv-test layer-json-test pmtiles-test decode-pmtiles-test overzoom-test ./unit suffixes = json json.gz @@ -266,6 +270,22 @@ enumerate-test: cmp tests/ne_110m_admin_0_countries/out/enum.check tests/ne_110m_admin_0_countries/out/enum rm tests/ne_110m_admin_0_countries/out/enum.mbtiles tests/ne_110m_admin_0_countries/out/enum.check +overzoom-test: tippecanoe-overzoom + # Basic operation + ./tippecanoe-overzoom -o tests/pbf/13-1310-3166.pbf tests/pbf/11-327-791.pbf 11/327/791 13/1310/3166 + ./tippecanoe-decode tests/pbf/13-1310-3166.pbf 13 1310 3166 > tests/pbf/13-1310-3166.pbf.json.check + cmp tests/pbf/13-1310-3166.pbf.json.check tests/pbf/13-1310-3166.pbf.json + rm tests/pbf/13-1310-3166.pbf tests/pbf/13-1310-3166.pbf.json.check + # Different detail and buffer, and attribute stripping + ./tippecanoe-overzoom -d8 -b30 -y NAME -y name -y scalerank -o tests/pbf/13-1310-3166-8-30.pbf tests/pbf/11-327-791.pbf 11/327/791 13/1310/3166 + ./tippecanoe-decode tests/pbf/13-1310-3166-8-30.pbf 13 1310 3166 > tests/pbf/13-1310-3166-8-30.pbf.json.check + cmp tests/pbf/13-1310-3166-8-30.pbf.json.check tests/pbf/13-1310-3166-8-30.pbf.json + rm tests/pbf/13-1310-3166-8-30.pbf tests/pbf/13-1310-3166-8-30.pbf.json.check + # No features in child tile + ./tippecanoe-overzoom -o tests/pbf/14-2616-6331.pbf tests/pbf/11-327-791.pbf 11/327/791 14/2616/6331 + cmp tests/pbf/14-2616-6331.pbf /dev/null + rm tests/pbf/14-2616-6331.pbf + join-test: tile-join ./tippecanoe -q -f -z12 -o tests/join-population/tabblock_06001420.mbtiles -YALAND10:'Land area' -L'{"file": "tests/join-population/tabblock_06001420.json", "description": "population"}' ./tippecanoe -q -f -Z5 -z10 -o tests/join-population/macarthur.mbtiles -l macarthur tests/join-population/macarthur.json diff --git a/README.md b/README.md index 45d27315c..e7a71b151 100644 --- a/README.md +++ b/README.md @@ -960,3 +960,25 @@ Join block geometries to employment attributes: ``` $ tippecanoe-json-tool -c in_wac_S000_JT00_2015.csv tl_2010_18157_tabblock10.sort.json > blocks-wac.json ``` + +tippecanoe-overzoom +=================== + +The `tippecanoe-overzoom` utility creates a vector tile from one of its parent tiles, +clipping and scaling the geometry from the parent tile and excluding features that +are clipped away. The idea is that if you create very high resolution tiles +(using `--extra-detail`) at a moderate zoom level, you can use `tippecanoe-overzoom` +to turn those into moderate detail tiles at high zoom levels, for the benefit of +renderers that cannot internally overzoom high-resolution tiles without losing +some of the precision. Running: + + tippecanoe-overzoom -o out.mvt.gz inz/inx/iny outz/outx/outy in.mvt.gz + +reads tile `inz/inx/iny` of `in.mvt.gz` and produces tile `outz/outx/outy` of `out.mvt.gz`. + +### Options + + * `-b` *buffer*: Set the tile buffer in the output tile (default 5) + * `-d` *detail*: Set the detail of the output tile (default 12) + * `-y` *attribute*: Retain the specified *attribute* in the output features. All attributes that are not named in a `-y` option will be removed. + diff --git a/clip.cpp b/clip.cpp new file mode 100644 index 000000000..678f6808c --- /dev/null +++ b/clip.cpp @@ -0,0 +1,567 @@ +#include +#include +#include +#include +#include +#include "geometry.hpp" +#include "errors.hpp" + +drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + mapbox::geometry::point min(minx, miny); + mapbox::geometry::point max(maxx, maxy); + mapbox::geometry::box bbox(min, max); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + mapbox::geometry::linear_ring ring; + for (size_t k = i; k < j; k++) { + ring.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); + } + + mapbox::geometry::linear_ring lr = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox); + + if (lr.size() > 0) { + for (size_t k = 0; k < lr.size(); k++) { + if (k == 0) { + out.push_back(draw(VT_MOVETO, lr[k].x, lr[k].y)); + } else { + out.push_back(draw(VT_LINETO, lr[k].x, lr[k].y)); + } + } + + if (lr.size() > 0 && lr[0] != lr[lr.size() - 1]) { + out.push_back(draw(VT_LINETO, lr[0].x, lr[0].y)); + } + } + + i = j - 1; + } else { + fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); + exit(EXIT_IMPOSSIBLE); + } + } + + return out; +} + +drawvec simple_clip_poly(drawvec &geom, int z, int buffer) { + long long area = 1LL << (32 - z); + long long clip_buffer = buffer * area / 256; + + return simple_clip_poly(geom, -clip_buffer, -clip_buffer, area + clip_buffer, area + clip_buffer); +} + +drawvec clip_point(drawvec &geom, int z, long long buffer) { + long long min = 0; + long long area = 1LL << (32 - z); + + min -= buffer * area / 256; + area += buffer * area / 256; + + return clip_point(geom, min, min, area, area); +} + +drawvec clip_point(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].x >= minx && geom[i].y >= miny && geom[i].x <= maxx && geom[i].y <= maxy) { + out.push_back(geom[i]); + } + } + + return out; +} + +drawvec clip_lines(drawvec &geom, int z, long long buffer) { + long long min = 0; + long long area = 1LL << (32 - z); + min -= buffer * area / 256; + area += buffer * area / 256; + + return clip_lines(geom, min, min, area, area); +} + +drawvec clip_lines(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (i > 0 && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO) && geom[i].op == VT_LINETO) { + double x1 = geom[i - 1].x; + double y1 = geom[i - 1].y; + + double x2 = geom[i - 0].x; + double y2 = geom[i - 0].y; + + int c = clip(&x1, &y1, &x2, &y2, minx, miny, maxx, maxy); + + if (c > 1) { // clipped + out.push_back(draw(VT_MOVETO, x1, y1)); + out.push_back(draw(VT_LINETO, x2, y2)); + out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); + } else if (c == 1) { // unchanged + out.push_back(geom[i]); + } else { // clipped away entirely + out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); + } + } else { + out.push_back(geom[i]); + } + } + + return out; +} + +#define INSIDE 0 +#define LEFT 1 +#define RIGHT 2 +#define BOTTOM 4 +#define TOP 8 + +static int computeOutCode(double x, double y, double xmin, double ymin, double xmax, double ymax) { + int code = INSIDE; + + if (x < xmin) { + code |= LEFT; + } else if (x > xmax) { + code |= RIGHT; + } + + if (y < ymin) { + code |= BOTTOM; + } else if (y > ymax) { + code |= TOP; + } + + return code; +} + +int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax) { + int outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); + int outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); + int accept = 0; + int changed = 0; + + while (1) { + if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop + accept = 1; + break; + } else if (outcode0 & outcode1) { // Bitwise AND is not 0. Trivially reject and get out of loop + break; + } else { + // failed both tests, so calculate the line segment to clip + // from an outside point to an intersection with clip edge + double x = *x0, y = *y0; + + // At least one endpoint is outside the clip rectangle; pick it. + int outcodeOut = outcode0 ? outcode0 : outcode1; + + // Now find the intersection point; + // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0) + if (outcodeOut & TOP) { // point is above the clip rectangle + x = *x0 + (*x1 - *x0) * (ymax - *y0) / (*y1 - *y0); + y = ymax; + } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle + x = *x0 + (*x1 - *x0) * (ymin - *y0) / (*y1 - *y0); + y = ymin; + } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle + y = *y0 + (*y1 - *y0) * (xmax - *x0) / (*x1 - *x0); + x = xmax; + } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle + y = *y0 + (*y1 - *y0) * (xmin - *x0) / (*x1 - *x0); + x = xmin; + } + + // Now we move outside point to intersection point to clip + // and get ready for next pass. + if (outcodeOut == outcode0) { + *x0 = x; + *y0 = y; + outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); + changed = 1; + } else { + *x1 = x; + *y1 = y; + outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); + changed = 1; + } + } + } + + if (accept == 0) { + return 0; + } else { + return changed + 1; + } +} + +static void decode_clipped(mapbox::geometry::multi_polygon &t, drawvec &out) { + out.clear(); + + for (size_t i = 0; i < t.size(); i++) { + for (size_t j = 0; j < t[i].size(); j++) { + drawvec ring; + + for (size_t k = 0; k < t[i][j].size(); k++) { + ring.push_back(draw((k == 0) ? VT_MOVETO : VT_LINETO, t[i][j][k].x, t[i][j][k].y)); + } + + if (ring.size() > 0 && ring[ring.size() - 1] != ring[0]) { + fprintf(stderr, "Had to close ring\n"); + ring.push_back(draw(VT_LINETO, ring[0].x, ring[0].y)); + } + + double area = get_area(ring, 0, ring.size()); + + if ((j == 0 && area < 0) || (j != 0 && area > 0)) { + fprintf(stderr, "Ring area has wrong sign: %f for %zu\n", area, j); + exit(EXIT_IMPOSSIBLE); + } + + for (size_t k = 0; k < ring.size(); k++) { + out.push_back(ring[k]); + } + } + } +} + +drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip) { + mapbox::geometry::wagyu::wagyu wagyu; + + geom = remove_noop(geom, VT_POLYGON, 0); + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j >= i + 4) { + mapbox::geometry::linear_ring lr; + + for (size_t k = i; k < j; k++) { + lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); + } + + if (lr.size() >= 3) { + wagyu.add_ring(lr); + } + } + + i = j - 1; + } + } + + if (clip) { + long long area = 0xFFFFFFFF; + if (z != 0) { + area = 1LL << (32 - z); + } + long long clip_buffer = buffer * area / 256; + + mapbox::geometry::linear_ring lr; + + lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); + lr.push_back(mapbox::geometry::point(-clip_buffer, area + clip_buffer)); + lr.push_back(mapbox::geometry::point(area + clip_buffer, area + clip_buffer)); + lr.push_back(mapbox::geometry::point(area + clip_buffer, -clip_buffer)); + lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); + + wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); + } + + mapbox::geometry::multi_polygon result; + try { + wagyu.execute(mapbox::geometry::wagyu::clip_type_union, result, mapbox::geometry::wagyu::fill_type_positive, mapbox::geometry::wagyu::fill_type_positive); + } catch (std::runtime_error &e) { + FILE *f = fopen("/tmp/wagyu.log", "w"); + fprintf(f, "%s\n", e.what()); + fprintf(stderr, "%s\n", e.what()); + fprintf(f, "["); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j >= i + 4) { + mapbox::geometry::linear_ring lr; + + if (i != 0) { + fprintf(f, ","); + } + fprintf(f, "["); + + for (size_t k = i; k < j; k++) { + lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); + if (k != i) { + fprintf(f, ","); + } + fprintf(f, "[%lld,%lld]", geom[k].x, geom[k].y); + } + + fprintf(f, "]"); + + if (lr.size() >= 3) { + } + } + + i = j - 1; + } + } + + fprintf(f, "]"); + fprintf(f, "\n\n\n\n\n"); + + fclose(f); + fprintf(stderr, "Internal error: Polygon cleaning failed. Log in /tmp/wagyu.log\n"); + exit(EXIT_IMPOSSIBLE); + } + + drawvec ret; + decode_clipped(result, ret); + return ret; +} + +void to_tile_scale(drawvec &geom, int z, int detail) { + if (32 - detail - z < 0) { + for (size_t i = 0; i < geom.size(); i++) { + geom[i].x = std::round((double) geom[i].x * (1LL << (-(32 - detail - z)))); + geom[i].y = std::round((double) geom[i].y * (1LL << (-(32 - detail - z)))); + } + } else { + for (size_t i = 0; i < geom.size(); i++) { + geom[i].x = std::round((double) geom[i].x / (1LL << (32 - detail - z))); + geom[i].y = std::round((double) geom[i].y / (1LL << (32 - detail - z))); + } + } +} + +drawvec from_tile_scale(drawvec const &geom, int z, int detail) { + drawvec out; + for (size_t i = 0; i < geom.size(); i++) { + draw d = geom[i]; + d.x *= (1LL << (32 - detail - z)); + d.y *= (1LL << (32 - detail - z)); + out.push_back(d); + } + return out; +} + +drawvec remove_noop(drawvec geom, int type, int shift) { + // first pass: remove empty linetos + + long long x = 0, y = 0; + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_LINETO && (long long) std::round((double) geom[i].x / (1LL << shift)) == x && (long long) std::round((double) geom[i].y / (1LL << shift)) == y) { + continue; + } + + if (geom[i].op == VT_CLOSEPATH) { + out.push_back(geom[i]); + } else { /* moveto or lineto */ + out.push_back(geom[i]); + x = std::round((double) geom[i].x / (1LL << shift)); + y = std::round((double) geom[i].y / (1LL << shift)); + } + } + + // second pass: remove unused movetos + + if (type != VT_POINT) { + geom = out; + out.resize(0); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + if (i + 1 >= geom.size()) { + continue; + } + + if (geom[i + 1].op == VT_MOVETO) { + continue; + } + + if (geom[i + 1].op == VT_CLOSEPATH) { + fprintf(stderr, "Shouldn't happen\n"); + i++; // also remove unused closepath + continue; + } + } + + out.push_back(geom[i]); + } + } + + // second pass: remove empty movetos + + if (type == VT_LINE) { + geom = out; + out.resize(0); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + if (i > 0 && geom[i - 1].op == VT_LINETO && (long long) std::round((double) geom[i - 1].x / (1LL << shift)) == (long long) std::round((double) geom[i].x / (1LL << shift)) && (long long) std::round((double) geom[i - 1].y / (1LL << shift)) == (long long) std::round((double) geom[i].y / (1LL << shift))) { + continue; + } + } + + out.push_back(geom[i]); + } + } + + return out; +} + +double get_area_scaled(const drawvec &geom, size_t i, size_t j) { + const double max_exact_double = (double) ((1LL << 53) - 1); + + // keep scaling the geometry down until we can calculate its area without overflow + for (long long scale = 2; scale < (1LL << 30); scale *= 2) { + long long bx = geom[i].x; + long long by = geom[i].y; + bool again = false; + + // https://en.wikipedia.org/wiki/Shoelace_formula + double area = 0; + for (size_t k = i; k < j; k++) { + area += (double) ((geom[k].x - bx) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].y - by) / scale); + if (std::fabs(area) >= max_exact_double) { + again = true; + break; + } + area -= (double) ((geom[k].y - by) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].x - bx) / scale); + if (std::fabs(area) >= max_exact_double) { + again = true; + break; + } + } + + if (again) { + continue; + } else { + area /= 2; + return area * scale * scale; + } + } + + fprintf(stderr, "get_area_scaled: can't happen\n"); + exit(EXIT_IMPOSSIBLE); +} + +double get_area(const drawvec &geom, size_t i, size_t j) { + const double max_exact_double = (double) ((1LL << 53) - 1); + + // Coordinates in `geom` are 40-bit integers, so there is no good way + // to multiply them without possible precision loss. Since they probably + // do not use the full precision, shift them nearer to the origin so + // their product is more likely to be exactly representable as a double. + // + // (In practice they are actually 34-bit integers: 32 bits for the + // Mercator world plane, plus another two bits so features can stick + // off either the left or right side. But that is still too many bits + // for the product to fit either in a 64-bit long long or in a + // double where the largest exact integer is 2^53.) + // + // If the intermediate calculation still exceeds 2^53, start trying to + // recalculate the area by scaling down the geometry. This will not + // produce as precise an area, but it will still be close, and the + // sign will be correct, which is more important, since the sign + // determines the winding order of the rings. We can then use that + // sign with this generally more precise area calculation. + + long long bx = geom[i].x; + long long by = geom[i].y; + + // https://en.wikipedia.org/wiki/Shoelace_formula + double area = 0; + bool overflow = false; + for (size_t k = i; k < j; k++) { + area += (double) (geom[k].x - bx) * (double) (geom[i + ((k - i + 1) % (j - i))].y - by); + if (std::fabs(area) >= max_exact_double) { + overflow = true; + } + area -= (double) (geom[k].y - by) * (double) (geom[i + ((k - i + 1) % (j - i))].x - bx); + if (std::fabs(area) >= max_exact_double) { + overflow = true; + } + } + area /= 2; + + if (overflow) { + double scaled_area = get_area_scaled(geom, i, j); + if ((area < 0 && scaled_area > 0) || (area > 0 && scaled_area < 0)) { + area = -area; + } + } + + return area; +} + +double get_mp_area(drawvec &geom) { + double ret = 0; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + ret += get_area(geom, i, j); + i = j - 1; + } + } + + return ret; +} + +drawvec close_poly(drawvec &geom) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j - 1 > i) { + if (geom[j - 1].x != geom[i].x || geom[j - 1].y != geom[i].y) { + fprintf(stderr, "Internal error: polygon not closed\n"); + } + } + + for (size_t n = i; n < j - 1; n++) { + out.push_back(geom[n]); + } + out.push_back(draw(VT_CLOSEPATH, 0, 0)); + + i = j - 1; + } + } + + return out; +} diff --git a/geometry.cpp b/geometry.cpp index 082a55eda..a080012a4 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -12,8 +12,6 @@ #include #include #include -#include -#include #include #include "geometry.hpp" #include "projection.hpp" @@ -23,8 +21,6 @@ #include "errors.hpp" #include "projection.hpp" -static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax); - drawvec decode_geometry(char **meta, int z, unsigned tx, unsigned ty, long long *bbox, unsigned initial_x, unsigned initial_y) { drawvec out; @@ -83,331 +79,6 @@ drawvec decode_geometry(char **meta, int z, unsigned tx, unsigned ty, long long return out; } -void to_tile_scale(drawvec &geom, int z, int detail) { - for (size_t i = 0; i < geom.size(); i++) { - geom[i].x = std::round((double) geom[i].x / (1LL << (32 - detail - z))); - geom[i].y = std::round((double) geom[i].y / (1LL << (32 - detail - z))); - } -} - -drawvec from_tile_scale(drawvec const &geom, int z, int detail) { - drawvec out; - for (size_t i = 0; i < geom.size(); i++) { - draw d = geom[i]; - d.x *= (1LL << (32 - detail - z)); - d.y *= (1LL << (32 - detail - z)); - out.push_back(d); - } - return out; -} - -drawvec remove_noop(drawvec geom, int type, int shift) { - // first pass: remove empty linetos - - long long x = 0, y = 0; - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_LINETO && std::round((double) geom[i].x / (1LL << shift)) == x && std::round((double) geom[i].y / (1LL << shift)) == y) { - continue; - } - - if (geom[i].op == VT_CLOSEPATH) { - out.push_back(geom[i]); - } else { /* moveto or lineto */ - out.push_back(geom[i]); - x = std::round((double) geom[i].x / (1LL << shift)); - y = std::round((double) geom[i].y / (1LL << shift)); - } - } - - // second pass: remove unused movetos - - if (type != VT_POINT) { - geom = out; - out.resize(0); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - if (i + 1 >= geom.size()) { - continue; - } - - if (geom[i + 1].op == VT_MOVETO) { - continue; - } - - if (geom[i + 1].op == VT_CLOSEPATH) { - fprintf(stderr, "Shouldn't happen\n"); - i++; // also remove unused closepath - continue; - } - } - - out.push_back(geom[i]); - } - } - - // second pass: remove empty movetos - - if (type == VT_LINE) { - geom = out; - out.resize(0); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - if (i > 0 && geom[i - 1].op == VT_LINETO && std::round((double) geom[i - 1].x / (1LL << shift)) == std::round((double) geom[i].x / (1LL << shift)) && std::round((double) geom[i - 1].y / (1LL << shift)) == std::round((double) geom[i].y / (1LL << shift))) { - continue; - } - } - - out.push_back(geom[i]); - } - } - - return out; -} - -double get_area_scaled(const drawvec &geom, size_t i, size_t j) { - const double max_exact_double = (double) ((1LL << 53) - 1); - - // keep scaling the geometry down until we can calculate its area without overflow - for (long long scale = 2; scale < (1LL << 30); scale *= 2) { - long long bx = geom[i].x; - long long by = geom[i].y; - bool again = false; - - // https://en.wikipedia.org/wiki/Shoelace_formula - double area = 0; - for (size_t k = i; k < j; k++) { - area += (double) ((geom[k].x - bx) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].y - by) / scale); - if (std::fabs(area) >= max_exact_double) { - again = true; - break; - } - area -= (double) ((geom[k].y - by) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].x - bx) / scale); - if (std::fabs(area) >= max_exact_double) { - again = true; - break; - } - } - - if (again) { - continue; - } else { - area /= 2; - return area * scale * scale; - } - } - - fprintf(stderr, "get_area_scaled: can't happen\n"); - exit(EXIT_IMPOSSIBLE); -} - -double get_area(const drawvec &geom, size_t i, size_t j) { - const double max_exact_double = (double) ((1LL << 53) - 1); - - // Coordinates in `geom` are 40-bit integers, so there is no good way - // to multiply them without possible precision loss. Since they probably - // do not use the full precision, shift them nearer to the origin so - // their product is more likely to be exactly representable as a double. - // - // (In practice they are actually 34-bit integers: 32 bits for the - // Mercator world plane, plus another two bits so features can stick - // off either the left or right side. But that is still too many bits - // for the product to fit either in a 64-bit long long or in a - // double where the largest exact integer is 2^53.) - // - // If the intermediate calculation still exceeds 2^53, start trying to - // recalculate the area by scaling down the geometry. This will not - // produce as precise an area, but it will still be close, and the - // sign will be correct, which is more important, since the sign - // determines the winding order of the rings. We can then use that - // sign with this generally more precise area calculation. - - long long bx = geom[i].x; - long long by = geom[i].y; - - // https://en.wikipedia.org/wiki/Shoelace_formula - double area = 0; - bool overflow = false; - for (size_t k = i; k < j; k++) { - area += (double) (geom[k].x - bx) * (double) (geom[i + ((k - i + 1) % (j - i))].y - by); - if (std::fabs(area) >= max_exact_double) { - overflow = true; - } - area -= (double) (geom[k].y - by) * (double) (geom[i + ((k - i + 1) % (j - i))].x - bx); - if (std::fabs(area) >= max_exact_double) { - overflow = true; - } - } - area /= 2; - - if (overflow) { - double scaled_area = get_area_scaled(geom, i, j); - if ((area < 0 && scaled_area > 0) || (area > 0 && scaled_area < 0)) { - area = -area; - } - } - - return area; -} - -double get_mp_area(drawvec &geom) { - double ret = 0; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - ret += get_area(geom, i, j); - i = j - 1; - } - } - - return ret; -} - -static void decode_clipped(mapbox::geometry::multi_polygon &t, drawvec &out) { - out.clear(); - - for (size_t i = 0; i < t.size(); i++) { - for (size_t j = 0; j < t[i].size(); j++) { - drawvec ring; - - for (size_t k = 0; k < t[i][j].size(); k++) { - ring.push_back(draw((k == 0) ? VT_MOVETO : VT_LINETO, t[i][j][k].x, t[i][j][k].y)); - } - - if (ring.size() > 0 && ring[ring.size() - 1] != ring[0]) { - fprintf(stderr, "Had to close ring\n"); - ring.push_back(draw(VT_LINETO, ring[0].x, ring[0].y)); - } - - double area = get_area(ring, 0, ring.size()); - - if ((j == 0 && area < 0) || (j != 0 && area > 0)) { - fprintf(stderr, "Ring area has wrong sign: %f for %zu\n", area, j); - exit(EXIT_IMPOSSIBLE); - } - - for (size_t k = 0; k < ring.size(); k++) { - out.push_back(ring[k]); - } - } - } -} - -drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip) { - mapbox::geometry::wagyu::wagyu wagyu; - - geom = remove_noop(geom, VT_POLYGON, 0); - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j >= i + 4) { - mapbox::geometry::linear_ring lr; - - for (size_t k = i; k < j; k++) { - lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); - } - - if (lr.size() >= 3) { - wagyu.add_ring(lr); - } - } - - i = j - 1; - } - } - - if (clip) { - long long area = 0xFFFFFFFF; - if (z != 0) { - area = 1LL << (32 - z); - } - long long clip_buffer = buffer * area / 256; - - mapbox::geometry::linear_ring lr; - - lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); - lr.push_back(mapbox::geometry::point(-clip_buffer, area + clip_buffer)); - lr.push_back(mapbox::geometry::point(area + clip_buffer, area + clip_buffer)); - lr.push_back(mapbox::geometry::point(area + clip_buffer, -clip_buffer)); - lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); - - wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); - } - - mapbox::geometry::multi_polygon result; - try { - wagyu.execute(mapbox::geometry::wagyu::clip_type_union, result, mapbox::geometry::wagyu::fill_type_positive, mapbox::geometry::wagyu::fill_type_positive); - } catch (std::runtime_error &e) { - FILE *f = fopen("/tmp/wagyu.log", "w"); - fprintf(f, "%s\n", e.what()); - fprintf(stderr, "%s\n", e.what()); - fprintf(f, "["); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j >= i + 4) { - mapbox::geometry::linear_ring lr; - - if (i != 0) { - fprintf(f, ","); - } - fprintf(f, "["); - - for (size_t k = i; k < j; k++) { - lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); - if (k != i) { - fprintf(f, ","); - } - fprintf(f, "[%lld,%lld]", geom[k].x, geom[k].y); - } - - fprintf(f, "]"); - - if (lr.size() >= 3) { - } - } - - i = j - 1; - } - } - - fprintf(f, "]"); - fprintf(f, "\n\n\n\n\n"); - - fclose(f); - fprintf(stderr, "Internal error: Polygon cleaning failed. Log in /tmp/wagyu.log\n"); - exit(EXIT_IMPOSSIBLE); - } - - drawvec ret; - decode_clipped(result, ret); - return ret; -} - /* pnpoly: Copyright (c) 1970-2003, Wm. Randolph Franklin @@ -515,90 +186,6 @@ void check_polygon(drawvec &geom) { } } -drawvec close_poly(drawvec &geom) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j - 1 > i) { - if (geom[j - 1].x != geom[i].x || geom[j - 1].y != geom[i].y) { - fprintf(stderr, "Internal error: polygon not closed\n"); - } - } - - for (size_t n = i; n < j - 1; n++) { - out.push_back(geom[n]); - } - out.push_back(draw(VT_CLOSEPATH, 0, 0)); - - i = j - 1; - } - } - - return out; -} - -drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - mapbox::geometry::point min(minx, miny); - mapbox::geometry::point max(maxx, maxy); - mapbox::geometry::box bbox(min, max); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - mapbox::geometry::linear_ring ring; - for (size_t k = i; k < j; k++) { - ring.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); - } - - mapbox::geometry::linear_ring lr = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox); - - if (lr.size() > 0) { - for (size_t k = 0; k < lr.size(); k++) { - if (k == 0) { - out.push_back(draw(VT_MOVETO, lr[k].x, lr[k].y)); - } else { - out.push_back(draw(VT_LINETO, lr[k].x, lr[k].y)); - } - } - - if (lr.size() > 0 && lr[0] != lr[lr.size() - 1]) { - out.push_back(draw(VT_LINETO, lr[0].x, lr[0].y)); - } - } - - i = j - 1; - } else { - fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); - exit(EXIT_IMPOSSIBLE); - } - } - - return out; -} - -drawvec simple_clip_poly(drawvec &geom, int z, int buffer) { - long long area = 1LL << (32 - z); - long long clip_buffer = buffer * area / 256; - - return simple_clip_poly(geom, -clip_buffer, -clip_buffer, area + clip_buffer, area + clip_buffer); -} - drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area, serial_feature *this_feature, serial_feature *tiny_feature) { drawvec out; const double pixel = (1LL << (32 - detail - z)) * (double) tiny_polygon_size; @@ -716,28 +303,6 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double return out; } -drawvec clip_point(drawvec &geom, int z, long long buffer) { - long long min = 0; - long long area = 1LL << (32 - z); - - min -= buffer * area / 256; - area += buffer * area / 256; - - return clip_point(geom, min, min, area, area); -} - -drawvec clip_point(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].x >= minx && geom[i].y >= miny && geom[i].x <= maxx && geom[i].y <= maxy) { - out.push_back(geom[i]); - } - } - - return out; -} - int quick_check(long long *bbox, int z, long long buffer) { long long min = 0; long long area = 1LL << (32 - z); @@ -771,45 +336,6 @@ bool point_within_tile(long long x, long long y, int z) { return x >= 0 && y >= 0 && x < area && y < area; } -drawvec clip_lines(drawvec &geom, int z, long long buffer) { - long long min = 0; - long long area = 1LL << (32 - z); - min -= buffer * area / 256; - area += buffer * area / 256; - - return clip_lines(geom, min, min, area, area); -} - -drawvec clip_lines(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (i > 0 && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO) && geom[i].op == VT_LINETO) { - double x1 = geom[i - 1].x; - double y1 = geom[i - 1].y; - - double x2 = geom[i - 0].x; - double y2 = geom[i - 0].y; - - int c = clip(&x1, &y1, &x2, &y2, minx, miny, maxx, maxy); - - if (c > 1) { // clipped - out.push_back(draw(VT_MOVETO, x1, y1)); - out.push_back(draw(VT_LINETO, x2, y2)); - out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); - } else if (c == 1) { // unchanged - out.push_back(geom[i]); - } else { // clipped away entirely - out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); - } - } else { - out.push_back(geom[i]); - } - } - - return out; -} - static double square_distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) { long long p2x = segB_x - segA_x; long long p2y = segB_y - segA_y; @@ -1245,89 +771,6 @@ std::vector chop_polygon(std::vector &geoms) { } } -#define INSIDE 0 -#define LEFT 1 -#define RIGHT 2 -#define BOTTOM 4 -#define TOP 8 - -static int computeOutCode(double x, double y, double xmin, double ymin, double xmax, double ymax) { - int code = INSIDE; - - if (x < xmin) { - code |= LEFT; - } else if (x > xmax) { - code |= RIGHT; - } - - if (y < ymin) { - code |= BOTTOM; - } else if (y > ymax) { - code |= TOP; - } - - return code; -} - -static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax) { - int outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); - int outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); - int accept = 0; - int changed = 0; - - while (1) { - if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop - accept = 1; - break; - } else if (outcode0 & outcode1) { // Bitwise AND is not 0. Trivially reject and get out of loop - break; - } else { - // failed both tests, so calculate the line segment to clip - // from an outside point to an intersection with clip edge - double x = *x0, y = *y0; - - // At least one endpoint is outside the clip rectangle; pick it. - int outcodeOut = outcode0 ? outcode0 : outcode1; - - // Now find the intersection point; - // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0) - if (outcodeOut & TOP) { // point is above the clip rectangle - x = *x0 + (*x1 - *x0) * (ymax - *y0) / (*y1 - *y0); - y = ymax; - } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle - x = *x0 + (*x1 - *x0) * (ymin - *y0) / (*y1 - *y0); - y = ymin; - } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle - y = *y0 + (*y1 - *y0) * (xmax - *x0) / (*x1 - *x0); - x = xmax; - } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle - y = *y0 + (*y1 - *y0) * (xmin - *x0) / (*x1 - *x0); - x = xmin; - } - - // Now we move outside point to intersection point to clip - // and get ready for next pass. - if (outcodeOut == outcode0) { - *x0 = x; - *y0 = y; - outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); - changed = 1; - } else { - *x1 = x; - *y1 = y; - outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); - changed = 1; - } - } - } - - if (accept == 0) { - return 0; - } else { - return changed + 1; - } -} - drawvec stairstep(drawvec &geom, int z, int detail) { drawvec out; double scale = 1 << (32 - detail - z); diff --git a/geometry.hpp b/geometry.hpp index 0e49ead0f..50a86be18 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -67,6 +67,7 @@ drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip); drawvec simple_clip_poly(drawvec &geom, int z, int buffer); drawvec close_poly(drawvec &geom); drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area, serial_feature *this_feature, serial_feature *tiny_feature); +int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax); drawvec clip_lines(drawvec &geom, int z, long long buffer); drawvec stairstep(drawvec &geom, int z, int detail); bool point_within_tile(long long x, long long y, int z); diff --git a/man/tippecanoe.1 b/man/tippecanoe.1 index e3dc4ffd1..6903f895a 100644 --- a/man/tippecanoe.1 +++ b/man/tippecanoe.1 @@ -1219,3 +1219,29 @@ Join block geometries to employment attributes: $ tippecanoe\-json\-tool \-c in_wac_S000_JT00_2015.csv tl_2010_18157_tabblock10.sort.json > blocks\-wac.json .fi .RE +.SH tippecanoe\-overzoom +.PP +The \fB\fCtippecanoe\-overzoom\fR utility creates a vector tile from one of its parent tiles, +clipping and scaling the geometry from the parent tile and excluding features that +are clipped away. The idea is that if you create very high resolution tiles +(using \fB\fC\-\-extra\-detail\fR) at a moderate zoom level, you can use \fB\fCtippecanoe\-overzoom\fR +to turn those into moderate detail tiles at high zoom levels, for the benefit of +renderers that cannot internally overzoom high\-resolution tiles without losing +some of the precision. Running: +.PP +.RS +.nf +tippecanoe\-overzoom \-o out.mvt.gz inz/inx/iny outz/outx/outy in.mvt.gz +.fi +.RE +.PP +reads tile \fB\fCinz/inx/iny\fR of \fB\fCin.mvt.gz\fR and produces tile \fB\fCoutz/outx/outy\fR of \fB\fCout.mvt.gz\fR\&. +.SS Options +.RS +.IP \(bu 2 +\fB\fC\-b\fR \fIbuffer\fP: Set the tile buffer in the output tile (default 5) +.IP \(bu 2 +\fB\fC\-d\fR \fIdetail\fP: Set the detail of the output tile (default 12) +.IP \(bu 2 +\fB\fC\-y\fR \fIattribute\fP: Retain the specified \fIattribute\fP in the output features. All attributes that are not named in a \fB\fC\-y\fR option will be removed. +.RE diff --git a/overzoom.cpp b/overzoom.cpp new file mode 100644 index 000000000..b22399cec --- /dev/null +++ b/overzoom.cpp @@ -0,0 +1,217 @@ +#include +#include +#include +#include +#include +#include "errors.hpp" +#include "mvt.hpp" +#include "geometry.hpp" + +extern char *optarg; +extern int optind; + +int detail = 12; // tippecanoe-style: mvt extent == 1 << detail +int buffer = 5; // tippecanoe-style: mvt buffer == extent * buffer / 256; + +std::set keep; + +std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny) { + mvt_tile tile, outtile; + bool was_compressed; + + try { + if (!tile.decode(s, was_compressed)) { + fprintf(stderr, "Couldn't parse tile %d/%u/%u\n", oz, ox, oy); + exit(EXIT_MVT); + } + } catch (std::exception const &e) { + fprintf(stderr, "PBF decoding error in tile %d/%u/%u\n", oz, ox, oy); + exit(EXIT_PROTOBUF); + } + + for (auto const &layer : tile.layers) { + mvt_layer outlayer = mvt_layer(); + + outlayer.name = layer.name; + outlayer.version = layer.version; + outlayer.extent = 1LL << detail; + + for (auto const &feature : layer.features) { + mvt_feature outfeature; + drawvec geom; + int t = feature.type; + + // Convert feature geometry to world coordinates + + long long tilesize = 1LL << (32 - oz); // source tile size in world coordinates + draw ring_closure(0, 0, 0); + + for (auto const &g : feature.geometry) { + if (g.op == mvt_closepath) { + geom.push_back(ring_closure); + } else { + geom.emplace_back(g.op, + g.x * tilesize / layer.extent + ox * tilesize, + g.y * tilesize / layer.extent + oy * tilesize); + + if (g.op == mvt_moveto) { + ring_closure = geom.back(); + ring_closure.op = mvt_lineto; + } + } + } + + // Now offset from world coordinates to output tile coordinates, + // but retain world scale, because that is what tippecanoe clipping expects + + long long outtilesize = 1LL << (32 - nz); // destination tile size in world coordinates + for (auto &g : geom) { + g.x -= nx * outtilesize; + g.y -= ny * outtilesize; + } + + // Clip to output tile + + if (t == VT_LINE) { + geom = clip_lines(geom, nz, buffer); + } else if (t == VT_POLYGON) { + geom = simple_clip_poly(geom, nz, buffer); + } else if (t == VT_POINT) { + geom = clip_point(geom, nz, buffer); + } + + // Scale to output tile extent + + to_tile_scale(geom, nz, detail); + + // Clean geometries + + geom = remove_noop(geom, t, 0); + if (t == VT_POLYGON) { + geom = clean_or_clip_poly(geom, 0, 0, false); + geom = close_poly(geom); + } + + // Add geometry to output feature + + outfeature.type = t; + for (auto const &g : geom) { + outfeature.geometry.emplace_back(g.op, g.x, g.y); + } + + // ID and attributes, if it didn't get clipped away + + if (outfeature.geometry.size() > 0) { + if (feature.has_id) { + outfeature.has_id = true; + outfeature.id = feature.id; + } + + for (size_t i = 0; i + 1 < feature.tags.size(); i += 2) { + if (keep.size() == 0 || keep.find(layer.keys[feature.tags[i]]) != keep.end()) { + outlayer.tag(outfeature, layer.keys[feature.tags[i]], layer.values[feature.tags[i + 1]]); + } + } + + outlayer.features.push_back(outfeature); + } + } + + if (outlayer.features.size() > 0) { + outtile.layers.push_back(outlayer); + } + } + + if (outtile.layers.size() > 0) { + std::string pbf = outtile.encode(); + std::string compressed; + compress(pbf, compressed, true); + + return compressed; + } else { + return ""; + } +} + +void usage(char **argv) { + fprintf(stderr, "Usage: %s -o newtile.pbf.gz tile.pbf.gz oz/ox/oy nz/nx/ny\n", argv[0]); + fprintf(stderr, "to create tile nz/nx/ny from tile oz/ox/oy\n"); + exit(EXIT_FAILURE); +} + +int main(int argc, char **argv) { + int i; + const char *outfile = NULL; + + while ((i = getopt(argc, argv, "y:o:d:b:")) != -1) { + switch (i) { + case 'y': + keep.insert(optarg); + break; + + case 'o': + outfile = optarg; + break; + + case 'd': + detail = atoi(optarg); + break; + + case 'b': + buffer = atoi(optarg); + break; + + default: + usage(argv); + } + } + + if (argc - optind != 3) { + usage(argv); + } + + if (outfile == NULL) { + usage(argv); + } + + const char *infile = argv[optind + 0]; + + int oz, ox, oy; + if (sscanf(argv[optind + 1], "%d/%d/%d", &oz, &ox, &oy) != 3) { + fprintf(stderr, "%s: not in z/x/y form\n", argv[optind + 1]); + usage(argv); + } + + int nz, nx, ny; + if (sscanf(argv[optind + 2], "%d/%d/%d", &nz, &nx, &ny) != 3) { + fprintf(stderr, "%s: not in z/x/y form\n", argv[optind + 2]); + usage(argv); + } + + std::string tile; + char buf[1000]; + int len; + + FILE *f = fopen(infile, "rb"); + if (f == NULL) { + perror(infile); + exit(EXIT_FAILURE); + } + + while ((len = fread(buf, sizeof(char), 1000, f)) > 0) { + tile.append(std::string(buf, len)); + } + fclose(f); + + f = fopen(outfile, "wb"); + if (f == NULL) { + perror(outfile); + exit(EXIT_FAILURE); + } + + std::string out = overzoom(tile, oz, ox, oy, nz, nx, ny); + fwrite(out.c_str(), sizeof(char), out.size(), f); + fclose(f); + + return 0; +} diff --git a/tests/pbf/11-327-791.pbf b/tests/pbf/11-327-791.pbf new file mode 100644 index 000000000..f286f98c8 Binary files /dev/null and b/tests/pbf/11-327-791.pbf differ diff --git a/tests/pbf/13-1310-3166-8-30.pbf.json b/tests/pbf/13-1310-3166-8-30.pbf.json new file mode 100644 index 000000000..6c9d94c37 --- /dev/null +++ b/tests/pbf/13-1310-3166-8-30.pbf.json @@ -0,0 +1,17 @@ +{ "type": "FeatureCollection", "properties": { "zoom": 13, "x": 1310, "y": 3166 }, "features": [ +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_land", "version": 2, "extent": 256 }, "features": [ +{ "type": "Feature", "properties": { "scalerank": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -122.388897, 37.792151 ], [ -122.387867, 37.790930 ], [ -122.387009, 37.778720 ], [ -122.382545, 37.774514 ], [ -122.382545, 37.766101 ], [ -122.382889, 37.763794 ], [ -122.382545, 37.763116 ], [ -122.382545, 37.749272 ], [ -122.436790, 37.749272 ], [ -122.436790, 37.792151 ], [ -122.388897, 37.792151 ] ] ] } } +] } +, +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_populated_places", "version": 2, "extent": 256 }, "features": [ +{ "type": "Feature", "properties": { "NAME": "San Francisco" }, "geometry": { "type": "Point", "coordinates": [ -122.399540, 37.784283 ] } } +] } +, +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_roads", "version": 2, "extent": 256 }, "features": [ +{ "type": "Feature", "properties": { "scalerank": 4, "name": "101" }, "geometry": { "type": "LineString", "coordinates": [ [ -122.407265, 37.767865 ], [ -122.412586, 37.781841 ], [ -122.419453, 37.788760 ], [ -122.432156, 37.792151 ] ] } } +, +{ "type": "Feature", "properties": { "scalerank": 3, "name": "80" }, "geometry": { "type": "LineString", "coordinates": [ [ -122.407265, 37.767865 ], [ -122.388210, 37.785232 ], [ -122.382545, 37.789981 ] ] } } +, +{ "type": "Feature", "properties": { "scalerank": 4, "name": "280" }, "geometry": { "type": "LineString", "coordinates": [ [ -122.406406, 37.749272 ], [ -122.405548, 37.752258 ], [ -122.407265, 37.767865 ] ] } } +] } +] } diff --git a/tests/pbf/13-1310-3166.pbf.json b/tests/pbf/13-1310-3166.pbf.json new file mode 100644 index 000000000..bcfbc2869 --- /dev/null +++ b/tests/pbf/13-1310-3166.pbf.json @@ -0,0 +1,17 @@ +{ "type": "FeatureCollection", "properties": { "zoom": 13, "x": 1310, "y": 3166 }, "features": [ +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_land", "version": 2, "extent": 4096 }, "features": [ +{ "type": "Feature", "properties": { "featurecla": "Land", "scalerank": 0, "min_zoom": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -122.387760, 37.788760 ], [ -122.387094, 37.778754 ], [ -122.386837, 37.778517 ], [ -122.386837, 37.752665 ], [ -122.432499, 37.752665 ], [ -122.432499, 37.788760 ], [ -122.387760, 37.788760 ] ] ] } } +] } +, +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_populated_places", "version": 2, "extent": 4096 }, "features": [ +{ "type": "Feature", "properties": { "SCALERANK": 1, "NATSCALE": 300, "LABELRANK": 1, "FEATURECLA": "Populated place", "NAME": "San Francisco", "NAMEALT": "San Francisco-Oakland", "NAMEASCII": "San Francisco", "ADM0CAP": 0, "WORLDCITY": 1, "MEGACITY": 1, "SOV0NAME": "United States", "SOV_A3": "USA", "ADM0NAME": "United States of America", "ADM0_A3": "USA", "ADM1NAME": "California", "ISO_A2": "US", "LATITUDE": 37.769196, "LONGITUDE": -122.417169, "POP_MAX": 3450000, "POP_MIN": 732072, "POP_OTHER": 27400, "RANK_MAX": 12, "RANK_MIN": 11, "MEGANAME": "San Francisco-Oakland", "LS_NAME": "San Francisco1", "MAX_POP10": 988636, "MAX_POP20": 1130999, "MAX_POP50": 1371285, "MAX_POP300": 4561697, "MAX_POP310": 4561697, "MAX_NATSCA": 300, "MIN_AREAKM": 218, "MAX_AREAKM": 1748, "MIN_AREAMI": 84, "MAX_AREAMI": 675, "MIN_PERKM": 126, "MAX_PERKM": 755, "MIN_PERMI": 78, "MAX_PERMI": 469, "MIN_BBXMIN": -122.516667, "MAX_BBXMIN": -122.516667, "MIN_BBXMAX": -122.358333, "MAX_BBXMAX": -121.733333, "MIN_BBYMIN": 37.191667, "MAX_BBYMIN": 37.575, "MIN_BBYMAX": 37.816667, "MAX_BBYMAX": 38.041667, "MEAN_BBXC": -122.301354, "MEAN_BBYC": 37.622288, "TIMEZONE": "America/Los_Angeles", "UN_FID": 570, "POP1950": 1855, "POP1955": 2021, "POP1960": 2200, "POP1965": 2361, "POP1970": 2529, "POP1975": 2590, "POP1980": 2656, "POP1985": 2805, "POP1990": 2961, "POP1995": 3095, "POP2000": 3236, "POP2005": 3387, "POP2010": 3450, "POP2015": 3544, "POP2020": 3684, "POP2025": 3803, "POP2050": 3898, "MIN_ZOOM": 2.7, "WIKIDATAID": "Q62", "WOF_ID": 85922583, "CAPALT": 0, "NAME_EN": "San Francisco", "NAME_DE": "San Francisco", "NAME_ES": "San Francisco", "NAME_FR": "San Francisco", "NAME_PT": "São Francisco", "NAME_RU": "Сан-Франциско", "NAME_ZH": "旧金山", "NAME_AR": "سان فرانسيسكو", "NAME_BN": "সান ফ্রান্সিস্কো", "NAME_EL": "Σαν Φρανσίσκο", "NAME_HI": "सैन फ्रांसिस्को", "NAME_HU": "San Francisco", "NAME_ID": "San Francisco", "NAME_IT": "San Francisco", "NAME_JA": "サンフランシスコ", "NAME_KO": "샌프란시스코", "NAME_NL": "San Francisco", "NAME_PL": "San Francisco", "NAME_SV": "San Francisco", "NAME_TR": "San Francisco", "NAME_VI": "San Francisco", "NE_ID": 1159151479, "NAME_FA": "سان فرانسیسکو", "NAME_HE": "סן פרנסיסקו", "NAME_UK": "Сан-Франциско", "NAME_UR": "سان فرانسسکو", "NAME_ZHT": "舊金山", "GEONAMESID": 5391959 }, "geometry": { "type": "Point", "coordinates": [ -122.399583, 37.784249 ] } } +] } +, +{ "type": "FeatureCollection", "properties": { "layer": "ne_10m_roads", "version": 2, "extent": 4096 }, "features": [ +{ "type": "Feature", "properties": { "scalerank": 4, "featurecla": "Road", "type": "Major Highway", "sov_a3": "USA", "edited": "New in version 2.0.0", "name": "101", "question": 0, "length_km": 11, "toll": 0, "ne_part": "ne_1d4_original", "labelrank": 0, "ignore": 0, "add": 0, "rwdb_rd_id": 0, "orig_fid": 0, "uident": 118705, "continent": "North America", "expressway": 1, "level": "Federal", "min_zoom": 4, "min_label": 7 }, "geometry": { "type": "LineString", "coordinates": [ [ -122.407351, 37.767933 ], [ -122.412586, 37.781841 ], [ -122.419496, 37.788760 ] ] } } +, +{ "type": "Feature", "properties": { "scalerank": 3, "featurecla": "Road", "type": "Major Highway", "sov_a3": "USA", "edited": "New in version 2.0.0", "name": "80", "question": 0, "length_km": 14, "toll": 1, "ne_part": "ne_1d4_original", "labelrank": 0, "ignore": 0, "add": 0, "rwdb_rd_id": 0, "orig_fid": 0, "uident": 311505, "continent": "North America", "expressway": 1, "level": "Interstate", "min_zoom": 3, "min_label": 6 }, "geometry": { "type": "LineString", "coordinates": [ [ -122.407351, 37.767933 ], [ -122.388253, 37.785300 ], [ -122.386837, 37.786453 ] ] } } +, +{ "type": "Feature", "properties": { "scalerank": 4, "featurecla": "Road", "type": "Major Highway", "sov_a3": "USA", "edited": "New in version 2.0.0", "name": "280", "question": 0, "length_km": 52, "toll": 0, "ne_part": "ne_1d4_original", "labelrank": 0, "ignore": 0, "add": 0, "rwdb_rd_id": 0, "orig_fid": 0, "uident": 124805, "continent": "North America", "expressway": 1, "level": "Interstate", "min_zoom": 4, "min_label": 7 }, "geometry": { "type": "LineString", "coordinates": [ [ -122.405677, 37.752665 ], [ -122.407351, 37.767933 ] ] } } +] } +] } diff --git a/tile-join.cpp b/tile-join.cpp index 2984786a2..d393e2be0 100644 --- a/tile-join.cpp +++ b/tile-join.cpp @@ -129,6 +129,10 @@ void handle(std::string message, int z, unsigned x, unsigned y, std::map outlayer.extent) { + // this always scales up the existing layer instead of scaling down + // the layer that is being added, because the assumption is that + // scaling up is safe while scaling down requires geometry cleaning. + for (size_t i = 0; i < outlayer.features.size(); i++) { for (size_t j = 0; j < outlayer.features[i].geometry.size(); j++) { outlayer.features[i].geometry[j].x = outlayer.features[i].geometry[j].x * layer.extent / outlayer.extent; diff --git a/version.hpp b/version.hpp index 647182167..2c3fb7a8f 100644 --- a/version.hpp +++ b/version.hpp @@ -1,6 +1,6 @@ #ifndef VERSION_HPP #define VERSION_HPP -#define VERSION "v2.28.1" +#define VERSION "v2.29.0" #endif