Skip to content

Commit f76b76b

Browse files
ansismourner
andcommitted
fix albers edge cases and lambert center (#11229)
* fix albers edge cases and lambert center The albers implementation was not solid. Switch to one based on d3-geo (thanks!) which is much more robust. Lambert was not using the center longitude. * bring back albers wrap correction Co-authored-by: Vladimir Agafonkin <[email protected]>
1 parent ea73618 commit f76b76b

File tree

5 files changed

+32
-29
lines changed

5 files changed

+32
-29
lines changed

src/geo/projection/albers.js

Lines changed: 29 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -13,47 +13,50 @@ export default {
1313

1414
conic: true,
1515

16+
// based on https://github.com/d3/d3-geo-projection, MIT-licensed
17+
1618
initializeConstants() {
1719
if (this.constants && vec2.exactEquals(this.parallels, this.constants.parallels)) {
1820
return;
1921
}
2022

21-
const p1 = degToRad(this.parallels[0]);
22-
const p2 = degToRad(this.parallels[1]);
23-
const sinp1 = Math.sin(p1);
24-
const cosp1 = Math.cos(p1);
25-
const cosp12 = cosp1 * cosp1;
26-
const r = 0.5;
27-
const n = 0.5 * (sinp1 + Math.sin(p2));
28-
const c = cosp12 + 2 * n * sinp1;
29-
const b = r / n * Math.sqrt(c);
30-
31-
this.constants = {n, b, c, parallels: this.parallels};
23+
const sy0 = Math.sin(degToRad(this.parallels[0]));
24+
const n = (sy0 + Math.sin(degToRad(this.parallels[1]))) / 2;
25+
const c = 1 + sy0 * (2 * n - sy0);
26+
const r0 = Math.sqrt(c) / n;
27+
28+
this.constants = {n, c, r0, parallels: this.parallels};
3229
},
30+
3331
project(lng: number, lat: number) {
3432
this.initializeConstants();
3533

36-
const {n, b, c} = this.constants;
37-
const theta = n * degToRad(lng - this.center[0]);
38-
const a = 0.5 / n * Math.sqrt(c - 2 * n * Math.sin(degToRad(lat)));
39-
const x = a * Math.sin(theta);
40-
const y = b - a * Math.cos(theta);
34+
const lambda = degToRad(lng - this.center[0]);
35+
const phi = degToRad(lat);
4136

42-
return {x: 1 + 0.5 * x, y: 1 - 0.5 * y};
37+
const {n, c, r0} = this.constants;
38+
const r = Math.sqrt(c - 2 * n * Math.sin(phi)) / n;
39+
const x = r * Math.sin(lambda * n);
40+
const y = r * Math.cos(lambda * n) - r0;
41+
return {x, y};
4342
},
43+
4444
unproject(x: number, y: number) {
4545
this.initializeConstants();
46+
const {n, c, r0} = this.constants;
4647

47-
const {n, b, c} = this.constants;
48-
const x_ = (x - 1) * 2;
49-
const y_ = (y - 1) * -2;
50-
const y2 = -(y_ - b);
48+
const r0y = r0 + y;
49+
let l = Math.atan2(x, Math.abs(r0y)) * Math.sign(r0y);
50+
if (r0y * n < 0) {
51+
l -= Math.PI * Math.sign(x) * Math.sign(r0y);
52+
}
5153
const dt = degToRad(this.center[0]) * n;
52-
const theta = wrap(Math.atan2(x_, y2), -Math.PI - dt, Math.PI - dt);
53-
const lng = radToDeg(theta / n) + this.center[0];
54-
const a = x_ / Math.sin(theta);
55-
const s = clamp((Math.pow(a / 0.5 * n, 2) - c) / (-2 * n), -1, 1);
56-
const lat = clamp(radToDeg(Math.asin(s)), -MAX_MERCATOR_LATITUDE, MAX_MERCATOR_LATITUDE);
54+
l = wrap(l, -Math.PI - dt, Math.PI - dt);
55+
56+
const lng = radToDeg(l / n) + this.center[0];
57+
const phi = Math.asin(clamp((c - (x * x + r0y * r0y) * n * n) / (2 * n), -1, 1));
58+
const lat = clamp(radToDeg(phi), -MAX_MERCATOR_LATITUDE, MAX_MERCATOR_LATITUDE);
59+
5760
return new LngLat(lng, lat);
5861
}
5962
};

src/geo/projection/cylindrical_equal_area.js

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ import {clamp, degToRad, radToDeg} from '../../util/util.js';
44
import {MAX_MERCATOR_LATITUDE} from '../mercator_coordinate.js';
55

66
export default function(phi: number) {
7-
const cosPhi = Math.cos(degToRad(phi));
7+
const cosPhi = Math.max(0.01, Math.cos(degToRad(phi)));
88
// scale coordinates between 0 and 1 to avoid constraint issues
99
const scale = 1 / (2 * Math.max(Math.PI * cosPhi, 1 / cosPhi));
1010

src/geo/projection/lambert.js

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ export default {
3838

3939
// based on https://github.com/d3/d3-geo, MIT-licensed
4040
lat = degToRad(lat);
41-
lng = degToRad(lng);
41+
lng = degToRad(lng - this.center[0]);
4242

4343
const epsilon = 1e-6;
4444
const {n, f} = this.constants;
@@ -73,7 +73,7 @@ export default {
7373

7474
if (fy * n < 0) l -= Math.PI * Math.sign(x) * signFy;
7575

76-
const lng = clamp(radToDeg(l / n), -180, 180);
76+
const lng = clamp(radToDeg(l / n) + this.center[0], -180, 180);
7777
const phi = 2 * Math.atan(Math.pow(f / r, 1 / n)) - halfPi;
7878
const lat = clamp(radToDeg(phi), -MAX_MERCATOR_LATITUDE, MAX_MERCATOR_LATITUDE);
7979

349 Bytes
Loading
-824 Bytes
Loading

0 commit comments

Comments
 (0)