Add minimal support for alternate input projections

CHANGELOG.md

@@ -1,3 +1,7 @@
+## 1.11.9
+
+* Add minimal support for alternate input projections (EPSG:3857).
+
 ## 1.11.8
 
 * Add an option to calculate the density of features as a feature attribute

README.md

@@ -90,6 +90,7 @@ Options
  * -D _detail_ or --low-detail=_detail_: Detail at lower zoom levels (default 12, for tile resolution of 4096)
  * -m _detail_ or --minimum-detail=_detail_: Minimum detail that it will try if tiles are too big at regular detail (default 7)
  * -b _pixels_ or --buffer=_pixels_: Buffer size where features are duplicated from adjacent tiles. Units are "screen pixels"--1/256th of the tile width or height. (default 5)
+ * -s _projection_ or --projection=_projection_: Specify the projection of the input data. Currently supported are EPSG:4326 (WGS84, the default) and EPSG:3857 (Web Mercator).
 
 All internal math is done in terms of a 32-bit tile coordinate system, so 1/(2^32) of the size of Earth,
 or about 1cm, is the smallest distinguishable distance. If _maxzoom_ + _detail_ > 32, no additional

decode.cpp

@@ -164,7 +164,7 @@ void handle(std::string message, int z, unsigned x, unsigned y, int describe) {
 					long long wy = scale * y + (scale / extent) * py;
 
 					double lat, lon;
-					tile2latlon(wx, wy, 32, &lat, &lon);
+					tile2lonlat(wx, wy, 32, &lon, &lat);
 
 					ops.push_back(lonlat(op, lon, lat, px, py));
 				} else {

geojson.cpp

@@ -90,7 +90,7 @@ long long parse_geometry(int t, json_object *j, long long *bbox, long long *fpos
 			long long x, y;
 			double lon = j->array[0]->number;
 			double lat = j->array[1]->number;
-			latlon2tile(lat, lon, 32, &x, &y);
+			projection(lon, lat, 32, &x, &y);
 
 			if (j->length > 2) {
 				static int warned = 0;

main.cpp

@@ -58,6 +58,7 @@ int geometry_scale = 0;
 
 int prevent[256];
 int additional[256];
+void (*projection)(double ix, double iy, int zoom, long long *ox, long long *oy) = lonlat2tile;
 
 struct source {
 	std::string layer;
@@ -1634,8 +1635,8 @@ int read_input(std::vector<source> &sources, char *fname, const char *layername,
 
 	double minlat = 0, minlon = 0, maxlat = 0, maxlon = 0, midlat = 0, midlon = 0;
 
-	tile2latlon(midx, midy, maxzoom, &maxlat, &minlon);
+	tile2lonlat(midx, midy, maxzoom, &maxlon, &minlat);
-	tile2latlon(midx + 1, midy + 1, maxzoom, &minlat, &maxlon);
+	tile2lonlat(midx + 1, midy + 1, maxzoom, &minlon, &maxlat);
 
 	midlat = (maxlat + minlat) / 2;
 	midlon = (maxlon + minlon) / 2;
@@ -1674,8 +1675,8 @@ int read_input(std::vector<source> &sources, char *fname, const char *layername,
 		file_bbox[3] = (1LL << 32) - 1;
 	}
 
-	tile2latlon(file_bbox[0], file_bbox[1], 32, &maxlat, &minlon);
+	tile2lonlat(file_bbox[0], file_bbox[1], 32, &minlon, &maxlat);
-	tile2latlon(file_bbox[2], file_bbox[3], 32, &minlat, &maxlon);
+	tile2lonlat(file_bbox[2], file_bbox[3], 32, &maxlon, &minlat);
 
 	if (midlat < minlat) {
 		midlat = minlat;
@@ -1705,6 +1706,17 @@ static bool has_name(struct option *long_options, int *pl) {
 	return false;
 }
 
+struct projection {
+	const char *name;
+	void (*project)(double ix, double iy, int zoom, long long *ox, long long *oy);
+};
+
+struct projection projections[] = {
+	{"EPSG:4326", lonlat2tile},
+	{"EPSG:3857", epsg3857totile},
+	{NULL, NULL},
+};
+
 int main(int argc, char **argv) {
 #ifdef MTRACE
 	mtrace();
@@ -1763,6 +1775,7 @@ int main(int argc, char **argv) {
 		{"gamma", required_argument, 0, 'g'},
 		{"prevent", required_argument, 0, 'p'},
 		{"additional", required_argument, 0, 'a'},
+		{"projection", required_argument, 0, 's'},
 
 		{"exclude-all", no_argument, 0, 'X'},
 		{"force", no_argument, 0, 'f'},
@@ -1811,7 +1824,7 @@ int main(int argc, char **argv) {
 		}
 	}
 
-	while ((i = getopt_long(argc, argv, "n:l:z:Z:B:d:D:m:o:x:y:r:b:t:g:p:a:XfFqvPL:A:", long_options, NULL)) != -1) {
+	while ((i = getopt_long(argc, argv, "n:l:z:Z:B:d:D:m:o:x:y:r:b:t:g:p:a:XfFqvPL:A:s:", long_options, NULL)) != -1) {
 		switch (i) {
 		case 0:
 			break;
@@ -1978,6 +1991,20 @@ int main(int argc, char **argv) {
 			read_parallel = 1;
 			break;
 
+		case 's': {
+			struct projection *p;
+			for (p = projections; p->name != NULL; p++) {
+				if (strcmp(p->name, optarg) == 0) {
+					projection = p->project;
+					break;
+				}
+			}
+			if (p == NULL) {
+				fprintf(stderr, "Unknown projection (-s): %s\n", optarg);
+				exit(EXIT_FAILURE);
+			}
+		} break;
+
 		default: {
 			int width = 7 + strlen(argv[0]);
 			fprintf(stderr, "Usage: %s", argv[0]);

main.hpp

@@ -13,3 +13,5 @@ extern int quiet;
 
 extern int CPUS;
 extern int TEMP_FILES;
+
+extern void (*projection)(double ix, double iy, int zoom, long long *ox, long long *oy);

man/tippecanoe.1

@@ -107,6 +107,8 @@ compensate for the larger marker, or \-Bf\fInumber\fP to allow at most \fInumber
 \-m \fIdetail\fP or \-\-minimum\-detail=\fIdetail\fP: Minimum detail that it will try if tiles are too big at regular detail (default 7)
 .IP \(bu 2
 \-b \fIpixels\fP or \-\-buffer=\fIpixels\fP: Buffer size where features are duplicated from adjacent tiles. Units are "screen pixels"\-\-1/256th of the tile width or height. (default 5)
+.IP \(bu 2
+\-s \fIprojection\fP or \-\-projection=\fIprojection\fP: Specify the projection of the input data. Currently supported are EPSG:4326 (WGS84, the default) and EPSG:3857 (Web Mercator).
 .RE
 .PP
 All internal math is done in terms of a 32\-bit tile coordinate system, so 1/(2 of the size of Earth,

projection.cpp

@@ -1,8 +1,9 @@
+#include <stdio.h>
 #include <math.h>
 #include "projection.hpp"
 
 // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
-void latlon2tile(double lat, double lon, int zoom, long long *x, long long *y) {
+void lonlat2tile(double lon, double lat, 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.
@@ -31,12 +32,22 @@ void latlon2tile(double lat, double lon, int zoom, long long *x, long long *y) {
 }
 
 // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
-void tile2latlon(long long x, long long y, int zoom, double *lat, double *lon) {
+void tile2lonlat(long long x, long long y, int zoom, double *lon, double *lat) {
 	unsigned long long n = 1LL << zoom;
 	*lon = 360.0 * x / n - 180.0;
 	*lat = atan(sinh(M_PI * (1 - 2.0 * y / n))) * 180.0 / M_PI;
 }
 
+void epsg3857totile(double ix, double iy, int zoom, long long *x, long long *y) {
+	*x = ix * (1LL << 31) / 6378137.0 / M_PI + (1LL << 31);
+	*y = ((1LL << 32) - 1) - (iy * (1LL << 31) / 6378137.0 / M_PI + (1LL << 31));
+
+	if (zoom != 0) {
+		*x >>= (32 - zoom);
+		*y >>= (32 - zoom);
+	}
+}
+
 unsigned long long encode(unsigned int wx, unsigned int wy) {
 	unsigned long long out = 0;
 

projection.hpp

@@ -1,4 +1,5 @@
-void latlon2tile(double lat, double lon, int zoom, long long *x, long long *y);
+void lonlat2tile(double lon, double lat, int zoom, long long *x, long long *y);
-void tile2latlon(long long x, long long y, int zoom, double *lat, double *lon);
+void epsg3857totile(double ix, double iy, int zoom, long long *x, long long *y);
+void tile2lonlat(long long x, long long y, int zoom, double *lon, double *lat);
 unsigned long long encode(unsigned int wx, unsigned int wy);
 void decode(unsigned long long index, unsigned *wx, unsigned *wy);

version.hpp

@@ -1 +1 @@
-#define VERSION "tippecanoe v1.11.8\n"
+#define VERSION "tippecanoe v1.11.9\n"