369 lines
11 KiB
JavaScript
369 lines
11 KiB
JavaScript
|
// QSC projection rewritten from the original PROJ4
|
||
|
|
// https://github.com/OSGeo/proj.4/blob/master/src/PJ_qsc.c
|
||
|
|
|
||
|
|
import {EPSLN, TWO_PI, SPI, HALF_PI, FORTPI} from '../constants/values';
|
||
|
|
|
||
|
|
/* constants */
|
||
|
|
var FACE_ENUM = {
|
||
|
|
FRONT: 1,
|
||
|
|
RIGHT: 2,
|
||
|
|
BACK: 3,
|
||
|
|
LEFT: 4,
|
||
|
|
TOP: 5,
|
||
|
|
BOTTOM: 6
|
||
|
|
};
|
||
|
|
|
||
|
|
var AREA_ENUM = {
|
||
|
|
AREA_0: 1,
|
||
|
|
AREA_1: 2,
|
||
|
|
AREA_2: 3,
|
||
|
|
AREA_3: 4
|
||
|
|
};
|
||
|
|
|
||
|
|
export function init() {
|
||
|
|
|
||
|
|
this.x0 = this.x0 || 0;
|
||
|
|
this.y0 = this.y0 || 0;
|
||
|
|
this.lat0 = this.lat0 || 0;
|
||
|
|
this.long0 = this.long0 || 0;
|
||
|
|
this.lat_ts = this.lat_ts || 0;
|
||
|
|
this.title = this.title || "Quadrilateralized Spherical Cube";
|
||
|
|
|
||
|
|
/* Determine the cube face from the center of projection. */
|
||
|
|
if (this.lat0 >= HALF_PI - FORTPI / 2.0) {
|
||
|
|
this.face = FACE_ENUM.TOP;
|
||
|
|
} else if (this.lat0 <= -(HALF_PI - FORTPI / 2.0)) {
|
||
|
|
this.face = FACE_ENUM.BOTTOM;
|
||
|
|
} else if (Math.abs(this.long0) <= FORTPI) {
|
||
|
|
this.face = FACE_ENUM.FRONT;
|
||
|
|
} else if (Math.abs(this.long0) <= HALF_PI + FORTPI) {
|
||
|
|
this.face = this.long0 > 0.0 ? FACE_ENUM.RIGHT : FACE_ENUM.LEFT;
|
||
|
|
} else {
|
||
|
|
this.face = FACE_ENUM.BACK;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Fill in useful values for the ellipsoid <-> sphere shift
|
||
|
|
* described in [LK12]. */
|
||
|
|
if (this.es !== 0) {
|
||
|
|
this.one_minus_f = 1 - (this.a - this.b) / this.a;
|
||
|
|
this.one_minus_f_squared = this.one_minus_f * this.one_minus_f;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// QSC forward equations--mapping lat,long to x,y
|
||
|
|
// -----------------------------------------------------------------
|
||
|
|
export function forward(p) {
|
||
|
|
var xy = {x: 0, y: 0};
|
||
|
|
var lat, lon;
|
||
|
|
var theta, phi;
|
||
|
|
var t, mu;
|
||
|
|
/* nu; */
|
||
|
|
var area = {value: 0};
|
||
|
|
|
||
|
|
// move lon according to projection's lon
|
||
|
|
p.x -= this.long0;
|
||
|
|
|
||
|
|
/* Convert the geodetic latitude to a geocentric latitude.
|
||
|
|
* This corresponds to the shift from the ellipsoid to the sphere
|
||
|
|
* described in [LK12]. */
|
||
|
|
if (this.es !== 0) {//if (P->es != 0) {
|
||
|
|
lat = Math.atan(this.one_minus_f_squared * Math.tan(p.y));
|
||
|
|
} else {
|
||
|
|
lat = p.y;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Convert the input lat, lon into theta, phi as used by QSC.
|
||
|
|
* This depends on the cube face and the area on it.
|
||
|
|
* For the top and bottom face, we can compute theta and phi
|
||
|
|
* directly from phi, lam. For the other faces, we must use
|
||
|
|
* unit sphere cartesian coordinates as an intermediate step. */
|
||
|
|
lon = p.x; //lon = lp.lam;
|
||
|
|
if (this.face === FACE_ENUM.TOP) {
|
||
|
|
phi = HALF_PI - lat;
|
||
|
|
if (lon >= FORTPI && lon <= HALF_PI + FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
theta = lon - HALF_PI;
|
||
|
|
} else if (lon > HALF_PI + FORTPI || lon <= -(HALF_PI + FORTPI)) {
|
||
|
|
area.value = AREA_ENUM.AREA_1;
|
||
|
|
theta = (lon > 0.0 ? lon - SPI : lon + SPI);
|
||
|
|
} else if (lon > -(HALF_PI + FORTPI) && lon <= -FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_2;
|
||
|
|
theta = lon + HALF_PI;
|
||
|
|
} else {
|
||
|
|
area.value = AREA_ENUM.AREA_3;
|
||
|
|
theta = lon;
|
||
|
|
}
|
||
|
|
} else if (this.face === FACE_ENUM.BOTTOM) {
|
||
|
|
phi = HALF_PI + lat;
|
||
|
|
if (lon >= FORTPI && lon <= HALF_PI + FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
theta = -lon + HALF_PI;
|
||
|
|
} else if (lon < FORTPI && lon >= -FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_1;
|
||
|
|
theta = -lon;
|
||
|
|
} else if (lon < -FORTPI && lon >= -(HALF_PI + FORTPI)) {
|
||
|
|
area.value = AREA_ENUM.AREA_2;
|
||
|
|
theta = -lon - HALF_PI;
|
||
|
|
} else {
|
||
|
|
area.value = AREA_ENUM.AREA_3;
|
||
|
|
theta = (lon > 0.0 ? -lon + SPI : -lon - SPI);
|
||
|
|
}
|
||
|
|
} else {
|
||
|
|
var q, r, s;
|
||
|
|
var sinlat, coslat;
|
||
|
|
var sinlon, coslon;
|
||
|
|
|
||
|
|
if (this.face === FACE_ENUM.RIGHT) {
|
||
|
|
lon = qsc_shift_lon_origin(lon, +HALF_PI);
|
||
|
|
} else if (this.face === FACE_ENUM.BACK) {
|
||
|
|
lon = qsc_shift_lon_origin(lon, +SPI);
|
||
|
|
} else if (this.face === FACE_ENUM.LEFT) {
|
||
|
|
lon = qsc_shift_lon_origin(lon, -HALF_PI);
|
||
|
|
}
|
||
|
|
sinlat = Math.sin(lat);
|
||
|
|
coslat = Math.cos(lat);
|
||
|
|
sinlon = Math.sin(lon);
|
||
|
|
coslon = Math.cos(lon);
|
||
|
|
q = coslat * coslon;
|
||
|
|
r = coslat * sinlon;
|
||
|
|
s = sinlat;
|
||
|
|
|
||
|
|
if (this.face === FACE_ENUM.FRONT) {
|
||
|
|
phi = Math.acos(q);
|
||
|
|
theta = qsc_fwd_equat_face_theta(phi, s, r, area);
|
||
|
|
} else if (this.face === FACE_ENUM.RIGHT) {
|
||
|
|
phi = Math.acos(r);
|
||
|
|
theta = qsc_fwd_equat_face_theta(phi, s, -q, area);
|
||
|
|
} else if (this.face === FACE_ENUM.BACK) {
|
||
|
|
phi = Math.acos(-q);
|
||
|
|
theta = qsc_fwd_equat_face_theta(phi, s, -r, area);
|
||
|
|
} else if (this.face === FACE_ENUM.LEFT) {
|
||
|
|
phi = Math.acos(-r);
|
||
|
|
theta = qsc_fwd_equat_face_theta(phi, s, q, area);
|
||
|
|
} else {
|
||
|
|
/* Impossible */
|
||
|
|
phi = theta = 0;
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Compute mu and nu for the area of definition.
|
||
|
|
* For mu, see Eq. (3-21) in [OL76], but note the typos:
|
||
|
|
* compare with Eq. (3-14). For nu, see Eq. (3-38). */
|
||
|
|
mu = Math.atan((12 / SPI) * (theta + Math.acos(Math.sin(theta) * Math.cos(FORTPI)) - HALF_PI));
|
||
|
|
t = Math.sqrt((1 - Math.cos(phi)) / (Math.cos(mu) * Math.cos(mu)) / (1 - Math.cos(Math.atan(1 / Math.cos(theta)))));
|
||
|
|
|
||
|
|
/* Apply the result to the real area. */
|
||
|
|
if (area.value === AREA_ENUM.AREA_1) {
|
||
|
|
mu += HALF_PI;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_2) {
|
||
|
|
mu += SPI;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_3) {
|
||
|
|
mu += 1.5 * SPI;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Now compute x, y from mu and nu */
|
||
|
|
xy.x = t * Math.cos(mu);
|
||
|
|
xy.y = t * Math.sin(mu);
|
||
|
|
xy.x = xy.x * this.a + this.x0;
|
||
|
|
xy.y = xy.y * this.a + this.y0;
|
||
|
|
|
||
|
|
p.x = xy.x;
|
||
|
|
p.y = xy.y;
|
||
|
|
return p;
|
||
|
|
}
|
||
|
|
|
||
|
|
// QSC inverse equations--mapping x,y to lat/long
|
||
|
|
// -----------------------------------------------------------------
|
||
|
|
export function inverse(p) {
|
||
|
|
var lp = {lam: 0, phi: 0};
|
||
|
|
var mu, nu, cosmu, tannu;
|
||
|
|
var tantheta, theta, cosphi, phi;
|
||
|
|
var t;
|
||
|
|
var area = {value: 0};
|
||
|
|
|
||
|
|
/* de-offset */
|
||
|
|
p.x = (p.x - this.x0) / this.a;
|
||
|
|
p.y = (p.y - this.y0) / this.a;
|
||
|
|
|
||
|
|
/* Convert the input x, y to the mu and nu angles as used by QSC.
|
||
|
|
* This depends on the area of the cube face. */
|
||
|
|
nu = Math.atan(Math.sqrt(p.x * p.x + p.y * p.y));
|
||
|
|
mu = Math.atan2(p.y, p.x);
|
||
|
|
if (p.x >= 0.0 && p.x >= Math.abs(p.y)) {
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
} else if (p.y >= 0.0 && p.y >= Math.abs(p.x)) {
|
||
|
|
area.value = AREA_ENUM.AREA_1;
|
||
|
|
mu -= HALF_PI;
|
||
|
|
} else if (p.x < 0.0 && -p.x >= Math.abs(p.y)) {
|
||
|
|
area.value = AREA_ENUM.AREA_2;
|
||
|
|
mu = (mu < 0.0 ? mu + SPI : mu - SPI);
|
||
|
|
} else {
|
||
|
|
area.value = AREA_ENUM.AREA_3;
|
||
|
|
mu += HALF_PI;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Compute phi and theta for the area of definition.
|
||
|
|
* The inverse projection is not described in the original paper, but some
|
||
|
|
* good hints can be found here (as of 2011-12-14):
|
||
|
|
* http://fits.gsfc.nasa.gov/fitsbits/saf.93/saf.9302
|
||
|
|
* (search for "Message-Id: <9302181759.AA25477 at fits.cv.nrao.edu>") */
|
||
|
|
t = (SPI / 12) * Math.tan(mu);
|
||
|
|
tantheta = Math.sin(t) / (Math.cos(t) - (1 / Math.sqrt(2)));
|
||
|
|
theta = Math.atan(tantheta);
|
||
|
|
cosmu = Math.cos(mu);
|
||
|
|
tannu = Math.tan(nu);
|
||
|
|
cosphi = 1 - cosmu * cosmu * tannu * tannu * (1 - Math.cos(Math.atan(1 / Math.cos(theta))));
|
||
|
|
if (cosphi < -1) {
|
||
|
|
cosphi = -1;
|
||
|
|
} else if (cosphi > +1) {
|
||
|
|
cosphi = +1;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Apply the result to the real area on the cube face.
|
||
|
|
* For the top and bottom face, we can compute phi and lam directly.
|
||
|
|
* For the other faces, we must use unit sphere cartesian coordinates
|
||
|
|
* as an intermediate step. */
|
||
|
|
if (this.face === FACE_ENUM.TOP) {
|
||
|
|
phi = Math.acos(cosphi);
|
||
|
|
lp.phi = HALF_PI - phi;
|
||
|
|
if (area.value === AREA_ENUM.AREA_0) {
|
||
|
|
lp.lam = theta + HALF_PI;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_1) {
|
||
|
|
lp.lam = (theta < 0.0 ? theta + SPI : theta - SPI);
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_2) {
|
||
|
|
lp.lam = theta - HALF_PI;
|
||
|
|
} else /* area.value == AREA_ENUM.AREA_3 */ {
|
||
|
|
lp.lam = theta;
|
||
|
|
}
|
||
|
|
} else if (this.face === FACE_ENUM.BOTTOM) {
|
||
|
|
phi = Math.acos(cosphi);
|
||
|
|
lp.phi = phi - HALF_PI;
|
||
|
|
if (area.value === AREA_ENUM.AREA_0) {
|
||
|
|
lp.lam = -theta + HALF_PI;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_1) {
|
||
|
|
lp.lam = -theta;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_2) {
|
||
|
|
lp.lam = -theta - HALF_PI;
|
||
|
|
} else /* area.value == AREA_ENUM.AREA_3 */ {
|
||
|
|
lp.lam = (theta < 0.0 ? -theta - SPI : -theta + SPI);
|
||
|
|
}
|
||
|
|
} else {
|
||
|
|
/* Compute phi and lam via cartesian unit sphere coordinates. */
|
||
|
|
var q, r, s;
|
||
|
|
q = cosphi;
|
||
|
|
t = q * q;
|
||
|
|
if (t >= 1) {
|
||
|
|
s = 0;
|
||
|
|
} else {
|
||
|
|
s = Math.sqrt(1 - t) * Math.sin(theta);
|
||
|
|
}
|
||
|
|
t += s * s;
|
||
|
|
if (t >= 1) {
|
||
|
|
r = 0;
|
||
|
|
} else {
|
||
|
|
r = Math.sqrt(1 - t);
|
||
|
|
}
|
||
|
|
/* Rotate q,r,s into the correct area. */
|
||
|
|
if (area.value === AREA_ENUM.AREA_1) {
|
||
|
|
t = r;
|
||
|
|
r = -s;
|
||
|
|
s = t;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_2) {
|
||
|
|
r = -r;
|
||
|
|
s = -s;
|
||
|
|
} else if (area.value === AREA_ENUM.AREA_3) {
|
||
|
|
t = r;
|
||
|
|
r = s;
|
||
|
|
s = -t;
|
||
|
|
}
|
||
|
|
/* Rotate q,r,s into the correct cube face. */
|
||
|
|
if (this.face === FACE_ENUM.RIGHT) {
|
||
|
|
t = q;
|
||
|
|
q = -r;
|
||
|
|
r = t;
|
||
|
|
} else if (this.face === FACE_ENUM.BACK) {
|
||
|
|
q = -q;
|
||
|
|
r = -r;
|
||
|
|
} else if (this.face === FACE_ENUM.LEFT) {
|
||
|
|
t = q;
|
||
|
|
q = r;
|
||
|
|
r = -t;
|
||
|
|
}
|
||
|
|
/* Now compute phi and lam from the unit sphere coordinates. */
|
||
|
|
lp.phi = Math.acos(-s) - HALF_PI;
|
||
|
|
lp.lam = Math.atan2(r, q);
|
||
|
|
if (this.face === FACE_ENUM.RIGHT) {
|
||
|
|
lp.lam = qsc_shift_lon_origin(lp.lam, -HALF_PI);
|
||
|
|
} else if (this.face === FACE_ENUM.BACK) {
|
||
|
|
lp.lam = qsc_shift_lon_origin(lp.lam, -SPI);
|
||
|
|
} else if (this.face === FACE_ENUM.LEFT) {
|
||
|
|
lp.lam = qsc_shift_lon_origin(lp.lam, +HALF_PI);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Apply the shift from the sphere to the ellipsoid as described
|
||
|
|
* in [LK12]. */
|
||
|
|
if (this.es !== 0) {
|
||
|
|
var invert_sign;
|
||
|
|
var tanphi, xa;
|
||
|
|
invert_sign = (lp.phi < 0 ? 1 : 0);
|
||
|
|
tanphi = Math.tan(lp.phi);
|
||
|
|
xa = this.b / Math.sqrt(tanphi * tanphi + this.one_minus_f_squared);
|
||
|
|
lp.phi = Math.atan(Math.sqrt(this.a * this.a - xa * xa) / (this.one_minus_f * xa));
|
||
|
|
if (invert_sign) {
|
||
|
|
lp.phi = -lp.phi;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
lp.lam += this.long0;
|
||
|
|
p.x = lp.lam;
|
||
|
|
p.y = lp.phi;
|
||
|
|
return p;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Helper function for forward projection: compute the theta angle
|
||
|
|
* and determine the area number. */
|
||
|
|
function qsc_fwd_equat_face_theta(phi, y, x, area) {
|
||
|
|
var theta;
|
||
|
|
if (phi < EPSLN) {
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
theta = 0.0;
|
||
|
|
} else {
|
||
|
|
theta = Math.atan2(y, x);
|
||
|
|
if (Math.abs(theta) <= FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_0;
|
||
|
|
} else if (theta > FORTPI && theta <= HALF_PI + FORTPI) {
|
||
|
|
area.value = AREA_ENUM.AREA_1;
|
||
|
|
theta -= HALF_PI;
|
||
|
|
} else if (theta > HALF_PI + FORTPI || theta <= -(HALF_PI + FORTPI)) {
|
||
|
|
area.value = AREA_ENUM.AREA_2;
|
||
|
|
theta = (theta >= 0.0 ? theta - SPI : theta + SPI);
|
||
|
|
} else {
|
||
|
|
area.value = AREA_ENUM.AREA_3;
|
||
|
|
theta += HALF_PI;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
return theta;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Helper function: shift the longitude. */
|
||
|
|
function qsc_shift_lon_origin(lon, offset) {
|
||
|
|
var slon = lon + offset;
|
||
|
|
if (slon < -SPI) {
|
||
|
|
slon += TWO_PI;
|
||
|
|
} else if (slon > +SPI) {
|
||
|
|
slon -= TWO_PI;
|
||
|
|
}
|
||
|
|
return slon;
|
||
|
|
}
|
||
|
|
|
||
|
|
export var names = ["Quadrilateralized Spherical Cube", "Quadrilateralized_Spherical_Cube", "qsc"];
|
||
|
|
export default {
|
||
|
|
init: init,
|
||
|
|
forward: forward,
|
||
|
|
inverse: inverse,
|
||
|
|
names: names
|
||
|
|
};
|
||
|
|
|