Merge pull request #122 from mapbox/dateline2

More fixes for large polygons and clipping at the antimeridian

geometry.cc

@@ -357,16 +357,27 @@ 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 += (long double) geom[k].x * (long double) 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].y * (long double) geom[i + ((k - i + 1) % (j - i))].x;
 			}
 			area = area / 2;
 
+			// 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.
+
+			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.
+
 				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 += area;
-				if (*accum_area > pixel * pixel) {
+					if (area > 0 && *accum_area > pixel * pixel) {
 						// XXX use centroid;
 
 						out.push_back(draw(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2));
@@ -378,7 +389,7 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double
 						*accum_area -= pixel * pixel;
 					}
 
-				if (area >= 0) {
+					if (area > 0) {
 						included_last_outer = false;
 					}
 				} else {
@@ -390,10 +401,11 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double
 
 					*reduced = false;
 
-				if (area >= 0) {
+					if (area > 0) {
 						included_last_outer = true;
 					}
 				}
+			}
 
 			i = j - 1;
 		} else {

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;
 

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();
+
+					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));
+						}
+					}
+
+					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));
 						}
-					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;