diff --git a/Makefile b/Makefile index ec14c4db..7b482892 100644 --- a/Makefile +++ b/Makefile @@ -77,7 +77,7 @@ tippecanoe-json-tool: jsontool.o jsonpull/jsonpull.o csv.o text.o geojson-loop.o unit: unit.o text.o sort.o mvt.o $(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread -tippecanoe-overzoom: overzoom.o mvt.o clip.o evaluator.o jsonpull/jsonpull.o text.o attribute.o +tippecanoe-overzoom: overzoom.o mvt.o clip.o evaluator.o jsonpull/jsonpull.o text.o attribute.o read_json.o projection.o $(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread -include $(wildcard *.d) diff --git a/clip.cpp b/clip.cpp index 3538da70..15d65776 100644 --- a/clip.cpp +++ b/clip.cpp @@ -1397,3 +1397,151 @@ std::string overzoom(std::vector const &tiles, int nz, int nx, int return ""; } } + +drawvec fix_polygon(const drawvec &geom, bool use_winding, bool reverse_winding) { + int outer = 1; + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_CLOSEPATH) { + outer = 1; + } else if (geom[i].op == VT_MOVETO) { + // Find the end of the ring + + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + // A polygon ring must contain at least three points + // (and really should contain four). If this one does + // not have any, avoid a division by zero trying to + // calculate the centroid below. + if (j - i < 1) { + i = j - 1; + outer = 0; + continue; + } + + // Make a temporary copy of the ring. + // Close it if it isn't closed. + + drawvec ring; + for (size_t a = i; a < j; a++) { + ring.push_back(geom[a]); + } + if (j - i != 0 && (ring[0].x != ring[j - i - 1].x || ring[0].y != ring[j - i - 1].y)) { + ring.push_back(ring[0]); + } + + // A polygon ring at this point should contain at least four points. + // Flesh it out with some vertex copies if it doesn't. + + while (ring.size() < 4) { + ring.push_back(ring[0]); + } + + // Reverse ring if winding order doesn't match + // inner/outer expectation + + bool reverse_ring = false; + if (use_winding) { + // GeoJSON winding is reversed from vector winding + reverse_ring = true; + } else if (reverse_winding) { + // GeoJSON winding is reversed from vector winding + reverse_ring = false; + } else { + double area = get_area(ring, 0, ring.size()); + if ((area > 0) != outer) { + reverse_ring = true; + } + } + + if (reverse_ring) { + drawvec tmp; + for (int a = ring.size() - 1; a >= 0; a--) { + tmp.push_back(ring[a]); + } + ring = tmp; + } + + // Now we are rotating the ring to make the first/last point + // one that would be unlikely to be simplified away. + + // calculate centroid + // a + 1 < size() because point 0 is duplicated at the end + long long xtotal = 0; + long long ytotal = 0; + long long count = 0; + for (size_t a = 0; a + 1 < ring.size(); a++) { + xtotal += ring[a].x; + ytotal += ring[a].y; + count++; + } + xtotal /= count; + ytotal /= count; + + // figure out which point is furthest from the centroid + long long dist2 = 0; + long long furthest = 0; + for (size_t a = 0; a + 1 < ring.size(); a++) { + // division by 16 because these are z0 coordinates and we need to avoid overflow + long long xd = (ring[a].x - xtotal) / 16; + long long yd = (ring[a].y - ytotal) / 16; + long long d2 = xd * xd + yd * yd; + if (d2 > dist2 || (d2 == dist2 && ring[a] < ring[furthest])) { + dist2 = d2; + furthest = a; + } + } + + // then figure out which point is furthest from *that*, + // which will hopefully be a good origin point since it should be + // at a far edge of the shape. + long long dist2b = 0; + long long furthestb = 0; + for (size_t a = 0; a + 1 < ring.size(); a++) { + // division by 16 because these are z0 coordinates and we need to avoid overflow + long long xd = (ring[a].x - ring[furthest].x) / 16; + long long yd = (ring[a].y - ring[furthest].y) / 16; + long long d2 = xd * xd + yd * yd; + if (d2 > dist2b || (d2 == dist2b && ring[a] < ring[furthestb])) { + dist2b = d2; + furthestb = a; + } + } + + // rotate ring so the furthest point is the duplicated one. + // the idea is that simplification will then be more efficient, + // never wasting the start and end points, which are always retained, + // on a point that has little impact on the shape. + + // Copy ring into output, fixing the moveto/lineto ops if necessary because of + // reversal or closing + + for (size_t a = 0; a < ring.size(); a++) { + size_t a2 = (a + furthestb) % (ring.size() - 1); + + if (a == 0) { + out.push_back(draw(VT_MOVETO, ring[a2].x, ring[a2].y)); + } else { + out.push_back(draw(VT_LINETO, ring[a2].x, ring[a2].y)); + } + } + + // Next ring or polygon begins on the non-lineto that ended this one + // and is not an outer ring unless there is a terminator first + + i = j - 1; + outer = 0; + } else { + fprintf(stderr, "Internal error: polygon ring begins with %d, not moveto\n", geom[i].op); + exit(EXIT_IMPOSSIBLE); + } + } + + return out; +} diff --git a/geometry.cpp b/geometry.cpp index a0330308..7246f6da 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -381,154 +381,6 @@ drawvec reorder_lines(const drawvec &geom) { return geom; } -drawvec fix_polygon(const drawvec &geom) { - int outer = 1; - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_CLOSEPATH) { - outer = 1; - } else if (geom[i].op == VT_MOVETO) { - // Find the end of the ring - - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - // A polygon ring must contain at least three points - // (and really should contain four). If this one does - // not have any, avoid a division by zero trying to - // calculate the centroid below. - if (j - i < 1) { - i = j - 1; - outer = 0; - continue; - } - - // Make a temporary copy of the ring. - // Close it if it isn't closed. - - drawvec ring; - for (size_t a = i; a < j; a++) { - ring.push_back(geom[a]); - } - if (j - i != 0 && (ring[0].x != ring[j - i - 1].x || ring[0].y != ring[j - i - 1].y)) { - ring.push_back(ring[0]); - } - - // A polygon ring at this point should contain at least four points. - // Flesh it out with some vertex copies if it doesn't. - - while (ring.size() < 4) { - ring.push_back(ring[0]); - } - - // Reverse ring if winding order doesn't match - // inner/outer expectation - - bool reverse_ring = false; - if (prevent[P_USE_SOURCE_POLYGON_WINDING]) { - // GeoJSON winding is reversed from vector winding - reverse_ring = true; - } else if (prevent[P_REVERSE_SOURCE_POLYGON_WINDING]) { - // GeoJSON winding is reversed from vector winding - reverse_ring = false; - } else { - double area = get_area(ring, 0, ring.size()); - if ((area > 0) != outer) { - reverse_ring = true; - } - } - - if (reverse_ring) { - drawvec tmp; - for (int a = ring.size() - 1; a >= 0; a--) { - tmp.push_back(ring[a]); - } - ring = tmp; - } - - // Now we are rotating the ring to make the first/last point - // one that would be unlikely to be simplified away. - - // calculate centroid - // a + 1 < size() because point 0 is duplicated at the end - long long xtotal = 0; - long long ytotal = 0; - long long count = 0; - for (size_t a = 0; a + 1 < ring.size(); a++) { - xtotal += ring[a].x; - ytotal += ring[a].y; - count++; - } - xtotal /= count; - ytotal /= count; - - // figure out which point is furthest from the centroid - long long dist2 = 0; - long long furthest = 0; - for (size_t a = 0; a + 1 < ring.size(); a++) { - // division by 16 because these are z0 coordinates and we need to avoid overflow - long long xd = (ring[a].x - xtotal) / 16; - long long yd = (ring[a].y - ytotal) / 16; - long long d2 = xd * xd + yd * yd; - if (d2 > dist2 || (d2 == dist2 && ring[a] < ring[furthest])) { - dist2 = d2; - furthest = a; - } - } - - // then figure out which point is furthest from *that*, - // which will hopefully be a good origin point since it should be - // at a far edge of the shape. - long long dist2b = 0; - long long furthestb = 0; - for (size_t a = 0; a + 1 < ring.size(); a++) { - // division by 16 because these are z0 coordinates and we need to avoid overflow - long long xd = (ring[a].x - ring[furthest].x) / 16; - long long yd = (ring[a].y - ring[furthest].y) / 16; - long long d2 = xd * xd + yd * yd; - if (d2 > dist2b || (d2 == dist2b && ring[a] < ring[furthestb])) { - dist2b = d2; - furthestb = a; - } - } - - // rotate ring so the furthest point is the duplicated one. - // the idea is that simplification will then be more efficient, - // never wasting the start and end points, which are always retained, - // on a point that has little impact on the shape. - - // Copy ring into output, fixing the moveto/lineto ops if necessary because of - // reversal or closing - - for (size_t a = 0; a < ring.size(); a++) { - size_t a2 = (a + furthestb) % (ring.size() - 1); - - if (a == 0) { - out.push_back(draw(VT_MOVETO, ring[a2].x, ring[a2].y)); - } else { - out.push_back(draw(VT_LINETO, ring[a2].x, ring[a2].y)); - } - } - - // Next ring or polygon begins on the non-lineto that ended this one - // and is not an outer ring unless there is a terminator first - - i = j - 1; - outer = 0; - } else { - fprintf(stderr, "Internal error: polygon ring begins with %d, not moveto\n", geom[i].op); - exit(EXIT_IMPOSSIBLE); - } - } - - return out; -} - #if 0 std::vector chop_polygon(std::vector &geoms) { while (1) { diff --git a/geometry.hpp b/geometry.hpp index d85ddcac..ec6baf2d 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -83,7 +83,7 @@ int quick_check(const long long *bbox, int z, long long buffer); void douglas_peucker(drawvec &geom, int start, int n, double e, size_t kept, size_t retain, bool prevent_simplify_shared_nodes); drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool mark_tile_bounds, double simplification, size_t retain, drawvec const &shared_nodes, struct node *shared_nodes_map, size_t nodepos, std::string const &shared_nodes_bloom); drawvec reorder_lines(const drawvec &geom); -drawvec fix_polygon(const drawvec &geom); +drawvec fix_polygon(const drawvec &geom, bool use_winding, bool reverse_winding); std::vector chop_polygon(std::vector &geoms); void check_polygon(drawvec &geom); double get_area(const drawvec &geom, size_t i, size_t j); diff --git a/plugin.cpp b/plugin.cpp index 742222eb..cca70ed5 100644 --- a/plugin.cpp +++ b/plugin.cpp @@ -74,182 +74,6 @@ void *run_writer(void *a) { return NULL; } -// XXX deduplicate -static std::vector to_feature(drawvec &geom) { - std::vector out; - - for (size_t i = 0; i < geom.size(); i++) { - out.push_back(mvt_geometry(geom[i].op, geom[i].x, geom[i].y)); - } - - return out; -} - -std::vector parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent) { - std::map ret; - std::shared_ptr tile_stringpool = std::make_shared(); - - json_pull *jp = json_begin_file(fp); - while (1) { - json_object *j = json_read(jp); - if (j == NULL) { - if (jp->error != NULL) { - fprintf(stderr, "Filter output:%d: %s: ", jp->line, jp->error); - if (jp->root != NULL) { - json_context(jp->root); - } else { - fprintf(stderr, "\n"); - } - exit(EXIT_JSON); - } - - json_free(jp->root); - break; - } - - json_object *type = json_hash_get(j, "type"); - if (type == NULL || type->type != JSON_STRING) { - continue; - } - if (strcmp(type->value.string.string, "Feature") != 0) { - continue; - } - - json_object *geometry = json_hash_get(j, "geometry"); - if (geometry == NULL) { - fprintf(stderr, "Filter output:%d: filtered feature with no geometry: ", jp->line); - json_context(j); - json_free(j); - exit(EXIT_JSON); - } - - json_object *properties = json_hash_get(j, "properties"); - if (properties == NULL || (properties->type != JSON_HASH && properties->type != JSON_NULL)) { - fprintf(stderr, "Filter output:%d: feature without properties hash: ", jp->line); - json_context(j); - json_free(j); - exit(EXIT_JSON); - } - - json_object *geometry_type = json_hash_get(geometry, "type"); - if (geometry_type == NULL) { - fprintf(stderr, "Filter output:%d: null geometry (additional not reported): ", jp->line); - json_context(j); - exit(EXIT_JSON); - } - - if (geometry_type->type != JSON_STRING) { - fprintf(stderr, "Filter output:%d: geometry type is not a string: ", jp->line); - json_context(j); - exit(EXIT_JSON); - } - - json_object *coordinates = json_hash_get(geometry, "coordinates"); - if (coordinates == NULL || coordinates->type != JSON_ARRAY) { - fprintf(stderr, "Filter output:%d: feature without coordinates array: ", jp->line); - json_context(j); - exit(EXIT_JSON); - } - - int t; - for (t = 0; t < GEOM_TYPES; t++) { - if (strcmp(geometry_type->value.string.string, geometry_names[t]) == 0) { - break; - } - } - if (t >= GEOM_TYPES) { - fprintf(stderr, "Filter output:%d: Can't handle geometry type %s: ", jp->line, geometry_type->value.string.string); - json_context(j); - exit(EXIT_JSON); - } - - std::string layername = "unknown"; - json_object *tippecanoe = json_hash_get(j, "tippecanoe"); - json_object *layer = NULL; - if (tippecanoe != NULL) { - layer = json_hash_get(tippecanoe, "layer"); - if (layer != NULL && layer->type == JSON_STRING) { - layername = std::string(layer->value.string.string); - } - } - - if (ret.count(layername) == 0) { - mvt_layer l; - l.name = layername; - l.version = 2; - l.extent = extent; - - ret.insert(std::pair(layername, l)); - } - auto l = ret.find(layername); - - drawvec dv; - parse_geometry(t, coordinates, dv, VT_MOVETO, "Filter output", jp->line, j); - if (mb_geometry[t] == VT_POLYGON) { - dv = fix_polygon(dv); - } - - // Scale and offset geometry from global to tile - for (size_t i = 0; i < dv.size(); i++) { - long long scale = 1LL << (32 - z); - dv[i].x = std::round((dv[i].x - scale * x) * extent / (double) scale); - dv[i].y = std::round((dv[i].y - scale * y) * extent / (double) scale); - } - - if (mb_geometry[t] == VT_POLYGON) { - // we can try scaling up because these are tile coordinates - dv = clean_or_clip_poly(dv, 0, 0, false, true); - if (dv.size() < 3) { - dv.clear(); - } - } - dv = remove_noop(dv, mb_geometry[t], 0); - if (mb_geometry[t] == VT_POLYGON) { - dv = close_poly(dv); - } - - if (dv.size() > 0) { - mvt_feature feature; - feature.type = mb_geometry[t]; - feature.geometry = to_feature(dv); - - json_object *id = json_hash_get(j, "id"); - if (id != NULL && id->type == JSON_NUMBER) { - feature.id = id->value.number.number; - if (id->value.number.large_unsigned > 0) { - feature.id = id->value.number.large_unsigned; - } - feature.has_id = true; - } - - for (size_t i = 0; i < properties->value.object.length; i++) { - serial_val sv = stringify_value(properties->value.object.values[i], "Filter output", jp->line, j); - - // Nulls can be excluded here because this is the postfilter - // and it is nearly time to create the vector representation - - if (sv.type != mvt_null) { - mvt_value v = stringified_to_mvt_value(sv.type, sv.s.c_str(), tile_stringpool); - l->second.tag(feature, std::string(properties->value.object.keys[i]->value.string.string), v); - } - } - - l->second.features.push_back(feature); - } - - json_free(j); - } - - json_end(jp); - - std::vector final; - for (auto a : ret) { - final.push_back(a.second); - } - - return final; -} - // Reads from the postfilter std::vector parse_layers(int fd, int z, unsigned x, unsigned y, std::vector> *layermaps, size_t tiling_seg, std::vector> *layer_unmaps, int extent) { FILE *f = fdopen(fd, "r"); @@ -401,7 +225,7 @@ serial_feature parse_feature(json_pull *jp, int z, unsigned x, unsigned y, std:: drawvec dv; parse_geometry(t, coordinates, dv, VT_MOVETO, "Filter output", jp->line, j); if (mb_geometry[t] == VT_POLYGON) { - dv = fix_polygon(dv); + dv = fix_polygon(dv, false, false); } // Scale and offset geometry from global to tile diff --git a/read_json.cpp b/read_json.cpp index 45e13658..2d7bf7f1 100644 --- a/read_json.cpp +++ b/read_json.cpp @@ -166,3 +166,179 @@ serial_val stringify_value(json_object *value, const char *reading, int line, js return sv; } + +// XXX deduplicate +static std::vector to_feature(drawvec &geom) { + std::vector out; + + for (size_t i = 0; i < geom.size(); i++) { + out.push_back(mvt_geometry(geom[i].op, geom[i].x, geom[i].y)); + } + + return out; +} + +std::vector parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent) { + std::map ret; + std::shared_ptr tile_stringpool = std::make_shared(); + + json_pull *jp = json_begin_file(fp); + while (1) { + json_object *j = json_read(jp); + if (j == NULL) { + if (jp->error != NULL) { + fprintf(stderr, "Filter output:%d: %s: ", jp->line, jp->error); + if (jp->root != NULL) { + json_context(jp->root); + } else { + fprintf(stderr, "\n"); + } + exit(EXIT_JSON); + } + + json_free(jp->root); + break; + } + + json_object *type = json_hash_get(j, "type"); + if (type == NULL || type->type != JSON_STRING) { + continue; + } + if (strcmp(type->value.string.string, "Feature") != 0) { + continue; + } + + json_object *geometry = json_hash_get(j, "geometry"); + if (geometry == NULL) { + fprintf(stderr, "Filter output:%d: filtered feature with no geometry: ", jp->line); + json_context(j); + json_free(j); + exit(EXIT_JSON); + } + + json_object *properties = json_hash_get(j, "properties"); + if (properties == NULL || (properties->type != JSON_HASH && properties->type != JSON_NULL)) { + fprintf(stderr, "Filter output:%d: feature without properties hash: ", jp->line); + json_context(j); + json_free(j); + exit(EXIT_JSON); + } + + json_object *geometry_type = json_hash_get(geometry, "type"); + if (geometry_type == NULL) { + fprintf(stderr, "Filter output:%d: null geometry (additional not reported): ", jp->line); + json_context(j); + exit(EXIT_JSON); + } + + if (geometry_type->type != JSON_STRING) { + fprintf(stderr, "Filter output:%d: geometry type is not a string: ", jp->line); + json_context(j); + exit(EXIT_JSON); + } + + json_object *coordinates = json_hash_get(geometry, "coordinates"); + if (coordinates == NULL || coordinates->type != JSON_ARRAY) { + fprintf(stderr, "Filter output:%d: feature without coordinates array: ", jp->line); + json_context(j); + exit(EXIT_JSON); + } + + int t; + for (t = 0; t < GEOM_TYPES; t++) { + if (strcmp(geometry_type->value.string.string, geometry_names[t]) == 0) { + break; + } + } + if (t >= GEOM_TYPES) { + fprintf(stderr, "Filter output:%d: Can't handle geometry type %s: ", jp->line, geometry_type->value.string.string); + json_context(j); + exit(EXIT_JSON); + } + + std::string layername = "unknown"; + json_object *tippecanoe = json_hash_get(j, "tippecanoe"); + json_object *layer = NULL; + if (tippecanoe != NULL) { + layer = json_hash_get(tippecanoe, "layer"); + if (layer != NULL && layer->type == JSON_STRING) { + layername = std::string(layer->value.string.string); + } + } + + if (ret.count(layername) == 0) { + mvt_layer l; + l.name = layername; + l.version = 2; + l.extent = extent; + + ret.insert(std::pair(layername, l)); + } + auto l = ret.find(layername); + + drawvec dv; + parse_geometry(t, coordinates, dv, VT_MOVETO, "Filter output", jp->line, j); + if (mb_geometry[t] == VT_POLYGON) { + dv = fix_polygon(dv, false, false); + } + + // Scale and offset geometry from global to tile + for (size_t i = 0; i < dv.size(); i++) { + long long scale = 1LL << (32 - z); + dv[i].x = std::round((dv[i].x - scale * x) * extent / (double) scale); + dv[i].y = std::round((dv[i].y - scale * y) * extent / (double) scale); + } + + if (mb_geometry[t] == VT_POLYGON) { + // we can try scaling up because these are tile coordinates + dv = clean_or_clip_poly(dv, 0, 0, false, true); + if (dv.size() < 3) { + dv.clear(); + } + } + dv = remove_noop(dv, mb_geometry[t], 0); + if (mb_geometry[t] == VT_POLYGON) { + dv = close_poly(dv); + } + + if (dv.size() > 0) { + mvt_feature feature; + feature.type = mb_geometry[t]; + feature.geometry = to_feature(dv); + + json_object *id = json_hash_get(j, "id"); + if (id != NULL && id->type == JSON_NUMBER) { + feature.id = id->value.number.number; + if (id->value.number.large_unsigned > 0) { + feature.id = id->value.number.large_unsigned; + } + feature.has_id = true; + } + + for (size_t i = 0; i < properties->value.object.length; i++) { + serial_val sv = stringify_value(properties->value.object.values[i], "Filter output", jp->line, j); + + // Nulls can be excluded here because this is the postfilter + // and it is nearly time to create the vector representation + + if (sv.type != mvt_null) { + mvt_value v = stringified_to_mvt_value(sv.type, sv.s.c_str(), tile_stringpool); + l->second.tag(feature, std::string(properties->value.object.keys[i]->value.string.string), v); + } + } + + l->second.features.push_back(feature); + } + + json_free(j); + } + + json_end(jp); + + std::vector final; + for (auto a : ret) { + final.push_back(a.second); + } + + return final; +} diff --git a/read_json.hpp b/read_json.hpp index 1d9374ab..e2bd2184 100644 --- a/read_json.hpp +++ b/read_json.hpp @@ -12,5 +12,6 @@ extern int mb_geometry[GEOM_TYPES]; void json_context(json_object *j); void parse_geometry(int t, json_object *j, drawvec &out, int op, const char *fname, int line, json_object *feature); +std::vector parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent); serial_val stringify_value(json_object *value, const char *reading, int line, json_object *feature); diff --git a/serial.cpp b/serial.cpp index b043c102..6cc0b683 100644 --- a/serial.cpp +++ b/serial.cpp @@ -460,7 +460,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf, std:: // This has to happen after scaling so that the wraparound detection has happened first. // Otherwise the inner/outer calculation will be confused by bad geometries. if (sf.t == VT_POLYGON) { - scaled_geometry = fix_polygon(scaled_geometry); + scaled_geometry = fix_polygon(scaled_geometry, prevent[P_USE_SOURCE_POLYGON_WINDING], prevent[P_REVERSE_SOURCE_POLYGON_WINDING]); } for (auto &c : clipbboxes) {