-
Notifications
You must be signed in to change notification settings - Fork 5
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
3 changed files
with
38 additions
and
26 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,40 +1,49 @@ | ||
pj_add(pj_cupola, 'cupola', 'Cupola', '\n\tMisc., Sph., NoInv.'); | ||
|
||
// Source: https://github.com/OSGeo/PROJ/issues/2706 | ||
// Source: https://www.tandfonline.com/eprint/EE7Y8RK4GXA4ITWUTQPY/full?target=10.1080/23729333.2020.1862962 | ||
// See also: http://www.at-a-lanta.nl/weia/cupola.html | ||
|
||
function pj_cupola(P) { | ||
// parameters for Cupola | ||
var c1 = 0.5253; // part of the equator on intermediate sphere, default = 1 | ||
var c2 = 0.7264; // sin of angle of polarline; default = 1 | ||
var c3 = 0.4188; // height of the equator, can be negative, default = 0 | ||
var c4 = 22.00; // phi of centre projection in degrees, default = 0 | ||
var c5 = 0.9701; // stretch in plane, default = 1 | ||
var c6 = 11.023; // central meridian 11.023 degrees, also defines border of the map | ||
var c7 = 180.00; // degrees to the right of c6, default = 180 | ||
var r2 = 1.61885660611815; | ||
var w123 = 0.258701109423297; | ||
var c4r = 0.383972435438752; | ||
var pc = 0.559562341761853; | ||
var de = 0.5253; // part of the equator on intermediate sphere, default = 1 | ||
var dp = 0.7264; // sin of angle of polar line; default = 1 | ||
var ri = 1 / Math.sqrt(de * dp); | ||
// height of the equator, can be negative, default = 0 | ||
var he = 0.4188; | ||
// phi of projection center | ||
var phi0 = 22 * DEG_TO_RAD; | ||
|
||
// NOTE: Proj applies the +lon_0 (central meridian) parameter before passing | ||
// coordinates to the projection function. | ||
// TODO: Use 11.023 as default central meridian | ||
// // var lam0 = 11.023 * DEG_TO_RAD; | ||
|
||
var se = 0.9701; // stretch in plane, default = 1 | ||
// center of projection on intermediate sphere | ||
var pc = calcP(phi0); | ||
var qc = calcQ(0); | ||
var spc = sin(pc); | ||
var cpc = cos(pc); | ||
var c6r = DEG_TO_RAD * (c6-c7+180); // c6 in [rad], v of center of projection | ||
var qc = c1*c6r; | ||
|
||
P.es = 0; | ||
P.fwd = s_fwd; | ||
|
||
function calcP(phi) { | ||
return asin(dp * sin(phi) + he * sqrt(de * dp)); | ||
} | ||
|
||
function calcQ(lam) { | ||
return de * lam; | ||
} | ||
|
||
function s_fwd(lp, xy) { | ||
// p, q: radians on intermediate sphere | ||
var p = asin(c2 * sin(lp.phi) + w123); | ||
var p = calcP(lp.phi); | ||
var q = calcQ(lp.lam); | ||
var sp = sin(p); | ||
var cp = cos(p); | ||
var q = c1 * lp.lam; | ||
var sqqc = sin(q-qc); | ||
var cqqc = cos(q-qc); | ||
// k is Snyder's constant, taken as product with r2: | ||
var r2k = r2 * sqrt(2/(1+spc*sp+cpc*cp*cqqc)); | ||
xy.x = r2k*cp*sqqc*c5; | ||
xy.y = r2k*(cpc*sp-spc*cp*cqqc)/c5; | ||
var sqqc = sin(q - qc); | ||
var cqqc = cos(q - qc); | ||
var K = sqrt(2 / (1 + sin(pc) * sp + cpc * cp * cqqc)); | ||
xy.x = ri * K * cp * sqqc * se; | ||
xy.y = ri * K * (cpc * sp - spc * cp * cqqc) / se; | ||
} | ||
} |