diff --git a/geojson.c b/geojson.c index 041492c..c933657 100644 --- a/geojson.c +++ b/geojson.c @@ -65,21 +65,21 @@ int CPUS; int TEMP_FILES; void init_cpus() { - CPUS = sysconf(_SC_NPROCESSORS_ONLN); - if (CPUS < 1) { - CPUS = 1; - } + CPUS = sysconf(_SC_NPROCESSORS_ONLN); + if (CPUS < 1) { + CPUS = 1; + } - TEMP_FILES = 64; - struct rlimit rl; - if (getrlimit(RLIMIT_NOFILE, &rl) != 0) { - perror("getrlimit"); - } else { - TEMP_FILES = rl.rlim_cur / 3; - if (TEMP_FILES > CPUS * 4) { - TEMP_FILES = CPUS * 4; - } - } + TEMP_FILES = 64; + struct rlimit rl; + if (getrlimit(RLIMIT_NOFILE, &rl) != 0) { + perror("getrlimit"); + } else { + TEMP_FILES = rl.rlim_cur / 3; + if (TEMP_FILES > CPUS * 4) { + TEMP_FILES = CPUS * 4; + } + } } size_t fwrite_check(const void *ptr, size_t size, size_t nitems, FILE *stream, const char *fname) { diff --git a/geometry.cc b/geometry.cc index 586bfa4..36ba691 100644 --- a/geometry.cc +++ b/geometry.cc @@ -357,41 +357,53 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double double area = 0; for (unsigned k = i; k < j; k++) { - area += geom[k].x * geom[i + ((k - i + 1) % (j - i))].y; - area -= geom[k].y * geom[i + ((k - i + 1) % (j - i))].x; + area += (long double) geom[k].x * (long double) geom[i + ((k - i + 1) % (j - i))].y; + area -= (long double) geom[k].y * (long double) geom[i + ((k - i + 1) % (j - i))].x; } area = area / 2; - if (fabs(area) <= pixel * pixel || (area < 0 && !included_last_outer)) { - // printf("area is only %f vs %lld so using square\n", area, pixel * pixel); + // XXX There is an ambiguity here: If the area of a ring is 0 and it is followed by holes, + // we don't know whether the area-0 ring was a hole too or whether it was the outer ring + // that these subsequent holes are somehow being subtracted from. I hope that if a polygon + // was simplified down to nothing, its holes also became nothing. - *accum_area += area; - if (*accum_area > pixel * pixel) { - // XXX use centroid; + if (area != 0) { + // These are pixel coordinates, so area > 0 for the outer ring. + // If the outer ring of a polygon was reduced to a pixel, its + // inner rings must just have their area de-accumulated rather + // than being drawn since we don't really know where they are. - out.push_back(draw(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y - pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y + pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y + pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); + if (fabs(area) <= pixel * pixel || (area < 0 && !included_last_outer)) { + // printf("area is only %f vs %lld so using square\n", area, pixel * pixel); - *accum_area -= pixel * pixel; - } + *accum_area += area; + if (area > 0 && *accum_area > pixel * pixel) { + // XXX use centroid; - if (area >= 0) { - included_last_outer = false; - } - } else { - // printf("area is %f so keeping instead of %lld\n", area, pixel * pixel); + out.push_back(draw(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y - pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y + pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y + pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); - for (unsigned k = i; k <= j && k < geom.size(); k++) { - out.push_back(geom[k]); - } + *accum_area -= pixel * pixel; + } - *reduced = false; + if (area > 0) { + included_last_outer = false; + } + } else { + // printf("area is %f so keeping instead of %lld\n", area, pixel * pixel); - if (area >= 0) { - included_last_outer = true; + for (unsigned k = i; k <= j && k < geom.size(); k++) { + out.push_back(geom[k]); + } + + *reduced = false; + + if (area > 0) { + included_last_outer = true; + } } } diff --git a/projection.c b/projection.c index 839370e..ccd4f00 100644 --- a/projection.c +++ b/projection.c @@ -3,6 +3,23 @@ // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames void latlon2tile(double lat, double lon, int zoom, long long *x, long long *y) { + // Must limit latitude somewhere to prevent overflow. + // 89.9 degrees latitude is 0.621 worlds beyond the edge of the flat earth, + // hopefully far enough out that there are few expectations about the shape. + if (lat < -89.9) { + lat = -89.9; + } + if (lat > 89.9) { + lat = 89.9; + } + + if (lon < -360) { + lon = -360; + } + if (lon > 360) { + lon = 360; + } + double lat_rad = lat * M_PI / 180; unsigned long long n = 1LL << zoom; diff --git a/tile.cc b/tile.cc index 6daa00a..b21e6b0 100644 --- a/tile.cc +++ b/tile.cc @@ -550,11 +550,17 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi // shifted by 360 degrees, and then make sure both copies get clipped down to size. unsigned n = geom.size(); - for (unsigned i = 0; i < n; i++) { - geom.push_back(draw(geom[i].op, geom[i].x - (1LL << 32), geom[i].y)); + + if (bbox[0] < 0) { + for (unsigned i = 0; i < n; i++) { + geom.push_back(draw(geom[i].op, geom[i].x + (1LL << 32), geom[i].y)); + } } - for (unsigned i = 0; i < n; i++) { - geom.push_back(draw(geom[i].op, geom[i].x + (1LL << 32), geom[i].y)); + + if (bbox[2] > 1LL << 32) { + for (unsigned i = 0; i < n; i++) { + geom.push_back(draw(geom[i].op, geom[i].x - (1LL << 32), geom[i].y)); + } } bbox[0] = 0;