Refactor raw projections, handle projection options (#10913)

* refactor raw projections, handle projection options

* add unprojection to winkel tripel

* fix flow

* Pin chrome to version 91 (#10887)

* Pin to chrome version 91

* Pin chrome version for test-browser

* fix lint

* remove to superfluous sin calls

Co-authored-by: Arindam Bose <arindam.bose@mapbox.com>

LICENSE.txt

@@ -53,8 +53,10 @@ EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 --------------------------------------------------------------------------------
 
 Contains a portion of d3-color https://github.com/d3/d3-color
+Contains a portion of d3-geo https://github.com/d3/d3-geo
+Contains a portion of d3-geo-projection https://github.com/d3/d3-geo-projection
 
-Copyright 2010-2016 Mike Bostock
+Copyright 2010-2021 Mike Bostock
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification,

debug/projections.html

@@ -88,15 +88,13 @@
         alaska: 4,
         winkel: 1.2,
         mercator: 1.0,
-        sinusoidal: 1.0,
         wgs84: 1.0
     };
     const centers = {
         albers: [-122.414, 37.776],
         alaska: [-154, 63],
         winkel: [0, 0],
         mercator: [0, 0],
-        sinusoidal: [0, 0],
         wgs84: [0, 0]
     };
     const projection = el.options[el.selectedIndex].value;

src/data/load_geometry.js

@@ -4,7 +4,7 @@ import {warnOnce, clamp} from '../util/util.js';
 
 import EXTENT from './extent.js';
 import {lngFromMercatorX, latFromMercatorY} from '../geo/mercator_coordinate.js';
-import projections from '../geo/projection/index.js';
+import getProjection from '../geo/projection/index.js';
 import tileTransform from '../geo/projection/tile_transform.js';
 import Point from '@mapbox/point-geometry';
 import type {CanonicalTileID} from '../source/tile_id.js';
@@ -19,8 +19,8 @@ const MIN = -MAX - 1;
 
 let projection;
 
-export function setProjection(projectionName: string) {
-    projection = projections[projectionName];
+export function setProjection(config: {name: string} | string) {
+    projection = getProjection(config);
 }
 
 function clampPoint(point: Point) {
@@ -52,7 +52,7 @@ export default function loadGeometry(feature: VectorTileFeature, canonical?: Can
     const scale = EXTENT / featureExtent;
     let cs, z2;
     if (canonical) {
-        cs = tileTransform(canonical, projection.project);
+        cs = tileTransform(canonical, projection);
         z2 = Math.pow(2, canonical.z);
     }
 

src/geo/projection/albers.js

@@ -1,74 +1,48 @@
 // @flow
 import LngLat from '../lng_lat.js';
 
-const albersConstants = {
-    refLng: -96,
-    p1r: 29.5,
-    p2r: 45.5
-};
-
-const alaskaConstants = {
-    refLng: -154,
-    p1r: 55,
-    p2r: 65
-};
-
-function project(lng, lat, constants) {
-    const p1 = constants.p1r / 180 * Math.PI;
-    const p2 = constants.p2r / 180 * Math.PI;
-    const n = 0.5 * (Math.sin(p1) + Math.sin(p2));
-    const theta = n * ((lng - constants.refLng) / 180 * Math.PI);
-    const c = Math.pow(Math.cos(p1), 2) + 2 * n * Math.sin(p1);
-    const r = 0.5;
-    const a = r / n * Math.sqrt(c - 2 * n * Math.sin(lat / 180 * Math.PI));
-    const b = r / n * Math.sqrt(c - 2 * n * Math.sin(0 / 180 * Math.PI));
-    const x = a * Math.sin(theta);
-    const y = b - a * Math.cos(theta);
-    const ret = {x: 0.5 + 0.5 * x, y: 0.5 + 0.5 * -y};
-    ret.x += 0.5;
-    ret.y += 0.5;
-
-    return ret;
-}
-
-function unproject(x, y, constants) {
-    const p1 = constants.p1r / 180 * Math.PI;
-    const p2 = constants.p2r / 180 * Math.PI;
-    const n = 0.5 * (Math.sin(p1) + Math.sin(p2));
-    const c = Math.pow(Math.cos(p1), 2) + 2 * n * Math.sin(p1);
-    const r = 0.5;
-    const b = r / n * Math.sqrt(c - 2 * n * Math.sin(0 / 180 * Math.PI));
-    const x__ = x - 0.5;
-    const y__ = y - 0.5;
-    const x_ = (x__ - 0.5) * 2;
-    const y_ = (y__ - 0.5) * -2;
-    const y2 = -(y_ - b);
-    const theta = Math.atan2(x_, y2);
-    const lng = (theta / n * 180 / Math.PI) + constants.refLng;
-    const a = x_ / Math.sin(theta);
-    const lat = Math.asin((Math.pow(a / r * n, 2) - c) / (-2 * n)) * 180 / Math.PI;
-
-    return new LngLat(lng, lat);
-}
-
 export const albers = {
     name: 'albers',
     range: [3.5, 7],
+
+    center: [-96, 37.5],
+    parallels: [29.5, 45.5],
+
     project(lng: number, lat: number) {
-        return project(lng, lat, albersConstants);
+        const p1 = this.parallels[0] / 180 * Math.PI;
+        const p2 = this.parallels[1] / 180 * Math.PI;
+        const n = 0.5 * (Math.sin(p1) + Math.sin(p2));
+        const theta = n * ((lng - this.center[0]) / 180 * Math.PI);
+        const c = Math.pow(Math.cos(p1), 2) + 2 * n * Math.sin(p1);
+        const r = 0.5;
+        const a = r / n * Math.sqrt(c - 2 * n * Math.sin(lat / 180 * Math.PI));
+        const b = r / n * Math.sqrt(c - 2 * n * Math.sin(0 / 180 * Math.PI));
+        const x = a * Math.sin(theta);
+        const y = b - a * Math.cos(theta);
+        return {x: 1 + 0.5 * x, y: 1 - 0.5 * y};
     },
-    unproject: (x: number, y: number) => {
-        return unproject(x, y, albersConstants);
+    unproject(x: number, y: number) {
+        const p1 = this.parallels[0] / 180 * Math.PI;
+        const p2 = this.parallels[1] / 180 * Math.PI;
+        const n = 0.5 * (Math.sin(p1) + Math.sin(p2));
+        const c = Math.pow(Math.cos(p1), 2) + 2 * n * Math.sin(p1);
+        const r = 0.5;
+        const b = r / n * Math.sqrt(c - 2 * n * Math.sin(0 / 180 * Math.PI));
+        const x_ = (x - 1) * 2;
+        const y_ = (y - 1) * -2;
+        const y2 = -(y_ - b);
+        const theta = Math.atan2(x_, y2);
+        const lng = (theta / n * 180 / Math.PI) + this.center[0];
+        const a = x_ / Math.sin(theta);
+        const lat = Math.asin((Math.pow(a / r * n, 2) - c) / (-2 * n)) * 180 / Math.PI;
+        return new LngLat(lng, lat);
     }
 };
 
 export const alaska = {
+    ...albers,
     name: 'alaska',
     range: [4, 7],
-    project(lng: number, lat: number) {
-        return project(lng, lat, alaskaConstants);
-    },
-    unproject: (x: number, y: number) => {
-        return unproject(x, y, alaskaConstants);
-    }
+    center: [-154, 50],
+    parallels: [55, 65]
 };

src/geo/projection/index.js

@@ -1,23 +1,29 @@
 // @flow
 import {albers, alaska} from './albers.js';
 import mercator from './mercator.js';
-import sinusoidal from './sinusoidal.js';
 import wgs84 from './wgs84.js';
 import winkel from './winkelTripel.js';
 import LngLat from '../lng_lat.js';
 
 export type Projection = {
     name: string,
-    range: Array<number>,
-    project: (lng: number, lat: number, options?: Object) => {x: number, y: number},
+    center: [number, number],
+    range?: Array<number>,
+    project: (lng: number, lat: number) => {x: number, y: number},
     unproject: (x: number, y: number) => LngLat
 };
 
-export default {
+const projections = {
     albers,
     alaska,
     mercator,
-    sinusoidal,
     wgs84,
     winkel
 };
+
+export default function getProjection(config: {name: string} | string) {
+    if (typeof config === 'string') {
+        return projections[config];
+    }
+    return {...projections[config.name], ...config};
+}

src/geo/projection/mercator.js

@@ -1,13 +1,18 @@
 // @flow
-import MercatorCoordinate, {mercatorXfromLng, mercatorYfromLat} from '../mercator_coordinate.js';
+import {mercatorXfromLng, mercatorYfromLat, lngFromMercatorX, latFromMercatorY} from '../mercator_coordinate.js';
+import LngLat from '../lng_lat.js';
 
 export default {
     name: 'mercator',
-    range: [],
+    center: [0, 0],
     project(lng: number, lat: number) {
         const x = mercatorXfromLng(lng);
         const y = mercatorYfromLat(lat);
         return {x, y};
     },
-    unproject: (x: number, y: number) => new MercatorCoordinate(x, y).toLngLat()
+    unproject(x: number, y: number) {
+        return new LngLat(
+            lngFromMercatorX(x),
+            latFromMercatorY(y));
+    }
 };

src/geo/projection/tile_transform.js

@@ -1,8 +1,9 @@
 // @flow
 import {lngFromMercatorX, latFromMercatorY} from '../mercator_coordinate.js';
 import {number as interpolate} from '../../style-spec/util/interpolate.js';
+import type {Projection} from './index.js';
 
-export default function tileTransform(id: Object, project: function) {
+export default function tileTransform(id: Object, projection: Projection) {
     const s = Math.pow(2, -id.z);
     const x1 = (id.x) * s;
     const x2 = (id.x + 1) * s;
@@ -22,10 +23,10 @@ export default function tileTransform(id: Object, project: function) {
         const lng = lngFromMercatorX(interpolate(x1, x2, i / n));
         const lat = latFromMercatorY(interpolate(y1, y2, i / n));
 
-        const p0 = project(lng, lat1);
-        const p1 = project(lng, lat2);
-        const p2 = project(lng1, lat);
-        const p3 = project(lng2, lat);
+        const p0 = projection.project(lng, lat1);
+        const p1 = projection.project(lng, lat2);
+        const p2 = projection.project(lng1, lat);
+        const p3 = projection.project(lng2, lat);
 
         minX = Math.min(minX, p0.x, p1.x, p2.x, p3.x);
         maxX = Math.max(maxX, p0.x, p1.x, p2.x, p3.x);

src/geo/projection/wgs84.js

@@ -1,12 +1,17 @@
 // @flow
+import LngLat from '../lng_lat.js';
 
 export default {
     name: 'wgs84',
-    range: [],
+    center: [0, 0],
     project(lng: number, lat: number) {
         const x = 0.5 + lng / 360;
         const y = 0.5 - lat / 360;
         return {x, y};
     },
-    unproject: () => {}
+    unproject(x: number, y: number) {
+        return new LngLat(
+            (x - 0.5) * 360,
+            (0.5 - y) * 360);
+    }
 };

src/geo/projection/winkelTripel.js

@@ -1,8 +1,60 @@
 // @flow
-// leaving Winkel Tripel methods as noop for now while we decide on implementation issues
+import LngLat from '../lng_lat.js';
+
 export default {
     name: 'winkel',
+    center: [0, 0],
     range: [3.5, 7],
-    project: () => {},
-    unproject: () => {}
+
+    project(lng: number, lat: number) {
+        lat = lat / 180 * Math.PI;
+        lng = lng / 180 * Math.PI;
+        const phi1 = Math.acos(2 / Math.PI);
+        const alpha = Math.acos(Math.cos(lat) * Math.cos(lng / 2));
+        const x = 0.5 * (lng * Math.cos(phi1) + (2 * Math.cos(lat) * Math.sin(lng / 2)) / (Math.sin(alpha) / alpha)) || 0;
+        const y = 0.5 * (lat + Math.sin(lat) / (Math.sin(alpha) / alpha)) || 0;
+        return {
+            x: (x / Math.PI + 0.5) * 0.5,
+            y: 1 - (y / Math.PI + 0.5) * 0.5
+        };
+    },
+
+    unproject(x: number, y: number) {
+        // based on https://github.com/d3/d3-geo-projection, MIT-licensed
+        x = (2 * x - 0.5) * Math.PI;
+        y = (2 * (1 - y) - 0.5) * Math.PI;
+        let lambda = x;
+        let phi = y;
+        let i = 25;
+        const epsilon = 1e-6;
+        let dlambda = 0, dphi = 0;
+        do {
+            const cosphi = Math.cos(phi),
+                sinphi = Math.sin(phi),
+                sinphi2 = 2 * sinphi * cosphi,
+                sin2phi = sinphi * sinphi,
+                cos2phi = cosphi * cosphi,
+                coslambda2 = Math.cos(lambda / 2),
+                sinlambda2 = Math.sin(lambda / 2),
+                sinlambda = 2 * coslambda2 * sinlambda2,
+                sin2lambda2 = sinlambda2 * sinlambda2,
+                C = 1 - cos2phi * coslambda2 * coslambda2,
+                F = C ? 1 / C : 0,
+                E = C ? Math.acos(cosphi * coslambda2) * Math.sqrt(1 / C) : 0,
+                fx = 0.5 * (2 * E * cosphi * sinlambda2 + lambda * 2 / Math.PI) - x,
+                fy = 0.5 * (E * sinphi + phi) - y,
+                dxdlambda = 0.5 * F * (cos2phi * sin2lambda2 + E * cosphi * coslambda2 * sin2phi) + 1 / Math.PI,
+                dxdphi = F * (sinlambda * sinphi2 / 4 - E * sinphi * sinlambda2),
+                dydlambda = 0.125 * F * (sinphi2 * sinlambda2 - E * sinphi * cos2phi * sinlambda),
+                dydphi = 0.5 * F * (sin2phi * coslambda2 + E * sin2lambda2 * cosphi) + 0.5,
+                denominator = dxdphi * dydlambda - dydphi * dxdlambda;
+
+            dlambda = (fy * dxdphi - fx * dydphi) / denominator;
+            dphi = (fx * dydlambda - fy * dxdlambda) / denominator;
+            lambda -= dlambda;
+            phi -= dphi;
+        } while ((Math.abs(dlambda) > epsilon || Math.abs(dphi) > epsilon) && --i > 0);
+
+        return new LngLat(lambda * 180 / Math.PI, phi * 180 / Math.PI);
+    }
 };

src/geo/transform.js

@@ -3,7 +3,7 @@
 import LngLat from './lng_lat.js';
 import LngLatBounds from './lng_lat_bounds.js';
 import MercatorCoordinate, {mercatorXfromLng, mercatorYfromLat, mercatorZfromAltitude, latFromMercatorY} from './mercator_coordinate.js';
-import projections from './projection/index.js';
+import getProjection from './projection/index.js';
 import tileTransform from '../geo/projection/tile_transform.js';
 import Point from '@mapbox/point-geometry';
 import {wrap, clamp, radToDeg, degToRad, getAABBPointSquareDist, furthestTileCorner} from '../util/util.js';
@@ -93,6 +93,7 @@ class Transform {
     cameraElevationReference: ElevationReference;
     fogCullDistSq: ?number;
     _averageElevation: number;
+    projectionOptions: {name: string} | string;
     projection: Projection;
     _elevation: ?Elevation;
     _fov: number;
@@ -129,7 +130,8 @@ class Transform {
         this.setMaxBounds();
 
         if (!projection) projection = 'mercator';
-        this.projection = projections[projection];
+        this.projectionOptions = projection;
+        this.projection = getProjection(projection);
 
         this.width = 0;
         this.height = 0;
@@ -678,7 +680,7 @@ class Transform {
         const minRange = options.isTerrainDEM ? -maxRange : this._elevation ? this._elevation.getMinElevationBelowMSL() : 0;
 
         const aabbForTile = (z, x, y, wrap, min, max) => {
-            const tt = tileTransform({z, x, y}, this.projection.project);
+            const tt = tileTransform({z, x, y}, this.projection);
             const tx = tt.x / tt.scale;
             const ty = tt.y / tt.scale;
             const tx2 = tt.x2 / tt.scale;
@@ -1265,7 +1267,7 @@ class Transform {
             scaledX = unwrappedX * scale;
             scaledY = canonical.y * scale;
         } else {
-            const cs = tileTransform(canonical, this.projection.project);
+            const cs = tileTransform(canonical, this.projection);
             scale = 1;
             scaledX = cs.x;
             scaledY = cs.y;