/*
Copyright 2026
JXG.Geography contributors
This file is part of JXG.Geography, an extension for JSXGraph.
JXG.Geography is free software dual licensed under the GNU LGPL or MIT
License, at your choice — the same terms JSXGraph itself is offered under.
* https://www.gnu.org/licenses/lgpl-3.0.html
* https://opensource.org/licenses/MIT
Geographic data: Natural Earth (naturalearthdata.com), public domain.
*/
/**
* JXG.Geography — core namespace.
*
* Everything here is board- and element-independent: projections, spherical
* measures, the rotation and clipping chain, the shared clip cache, the
* dataset validator, and the two helpers that elements must not reimplement
* (maskedLayer, invert).
*
* Draft 0.1 — see SPEC.md.
*
*
* CONVENTIONS
*
* Positions are [lon, lat] in degrees, longitude first, latitude second.
* Angles are degrees everywhere except where a name says otherwise: only
* `samplesFor(spanRad, ...)` takes radians, and it says so.
*
* Rings are open: the closing point is not repeated. All 116 rings of the
* bundled dataset are stored that way and `validate` insists on it, because a
* repeated point is a zero-length segment that every winding measure has to
* step over.
*
* An outer ring runs counter-clockwise seen from outside the sphere, a hole
* runs clockwise. `orientationOf` reports which, and the bundled dataset holds
* to it: 114 outer rings positive, 2 holes negative, none the other way.
*
* Tolerances are in *device* pixels, not CSS pixels — see `devicePixels()`.
* On a retina screen the same tolerance therefore asks for twice the samples.
*
* Nothing here touches a board or an element. Every function can be called on
* its own, and the test suites do exactly that.
*
*
* KNOWN LIMITS
*
* Orientation cannot always be read, and the two ways it fails are different.
*
* A ring that circles a pole is *flagged*: `readableOrientation` returns false
* and the caller's stated winding is used instead. A 90-degree cap centred off
* the equator is the everyday case — it encircles a pole, its integral comes
* out at exactly zero, and this is why a night cap states its winding rather
* than having it derived.
*
* A ring larger than a hemisphere that circles *no* pole is **not** flagged
* and is answered wrongly: a 120-degree cap about 0/0 encircles the antipode,
* its integral measures the complementary region, and `orientationOf` reports
* -1 for a ring that runs counter-clockwise. That is a real gap, not a
* convention. No coastline ring is remotely that large — the biggest in the
* bundled dataset spans 2 steradians — but a synthetic one can be, and then
* the answer is silently wrong.
*
* The subsolar point ignores the equation of time and can be 4.1 degrees of
* longitude out, about 457 km at the equator.
*
* The dataset is Natural Earth 1:110 million. 29 countries are too small to
* have an outline at that scale and stand in as circles, which a map meant to
* be measured from should switch off.
* @alias JXG.Geography.RAD
* @alias JXG.Geography.version
*/
/**
* The shared namespace: projections, spherical geometry, datasets. Both the
* map and the globe are built on it, and the drawable elements measure with
* it.
*
* @namespace JXG.Geography
*/
(function (root) {
'use strict';
root.JXG = root.JXG || {};
var JXG = root.JXG;
JXG.Geography = JXG.Geography || {};
var G = JXG.Geography;
G.version = '0.1';
// ---------------------------------------------------------------- basics
var RAD = Math.PI / 180, DEG = 180 / Math.PI, TAU = 2 * Math.PI;
var R_EARTH = 6371.0088; // mean radius, km
G.RAD = RAD; G.DEG = DEG; G.EARTH_RADIUS = R_EARTH;
/**
* Clamp a number into a closed interval.
*
* @param {number} v The value.
* @param {number} a Lower bound.
* @param {number} b Upper bound.
* @returns {number} v, or the bound it exceeded.
* @alias JXG.Geography.clamp
*/
function clamp(v, a, b) {
if (v < a) { return a; }
if (v > b) { return b; }
return v;
}
G.clamp = clamp;
var warned = {};
/**
* Report a problem once per message.
*
* Repeating the same complaint on every frame is what makes a console
* useless, so each distinct message is passed on only the first time.
*
* @param {string} msg What went wrong, in a form the caller can act on.
* @returns {void}
*
* @alias JXG.Geography.warn
*/
function warn(msg) {
if (warned[msg]) { return; }
warned[msg] = true;
if (typeof JXG.warn === 'function') {
JXG.warn('JXG.Geography: ' + msg);
} else if (root.console) {
root.console.warn('JXG.Geography: ' + msg);
}
}
G.warn = warn;
/**
* Forget which messages have already been reported.
*
* Only of use to a test that wants to see a warning it has provoked twice.
*
* @returns {void}
* @alias JXG.Geography.resetWarnings
*/
G.resetWarnings = function () { warned = {}; };
// ------------------------------------------------------- spherical maths
/**
* A position as a unit vector.
*
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @returns {Array<number>} Unit vector [x, y, z]; x towards 0°/0°,
* z towards the north pole.
* @alias JXG.Geography.toVector
*/
function toVector(lon, lat) {
var l = lon * RAD, p = lat * RAD, c = Math.cos(p);
return [c * Math.cos(l), c * Math.sin(l), Math.sin(p)];
}
/**
* Inverse of toVector. The vector is normalised first: callers that hand in
* a sum of directions — a centroid, say — would otherwise get a latitude
* read off an unnormalised z, which is wrong by however long the vector is.
*
* @param {Array<number>} v A vector [x, y, z]; need not be a unit vector.
* @returns {Array<number>} The position as [lon, lat] in degrees.
*/
function toGeo(v) {
var L = Math.hypot(v[0], v[1], v[2]);
if (L < 1e-15) { return [0, 0]; }
return [Math.atan2(v[1], v[0]) * DEG,
Math.asin(clamp(v[2] / L, -1, 1)) * DEG];
}
function dot(a, b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; }
/**
* Angular distance between two positions.
*
* @param {Array<number>} a Position as [lon, lat] in degrees.
* @param {Array<number>} b Position as [lon, lat] in degrees.
* @returns {number} Angular distance in degrees, 0 to 180.
* @alias JXG.Geography.angle
*/
function angle(a, b) {
var u = toVector(a[0], a[1]), v = toVector(b[0], b[1]);
var c = [u[1] * v[2] - u[2] * v[1],
u[2] * v[0] - u[0] * v[2],
u[0] * v[1] - u[1] * v[0]];
// atan2 of |u x v| against u.v stays accurate for small and large angles,
// unlike acos(u.v) which loses precision near 0.
return Math.atan2(Math.hypot(c[0], c[1], c[2]), dot(u, v)) * DEG;
}
/**
* Great-circle distance in km.
*
* @param {Array<number>} a Position as [lon, lat] in degrees.
* @param {Array<number>} b Position as [lon, lat] in degrees.
* @returns {number} Great-circle distance in kilometres.
*/
function distance(a, b) { return angle(a, b) * RAD * R_EARTH; }
/**
* Initial course from a to b, degrees clockwise from north.
*
* @param {Array<number>} a Position as [lon, lat] in degrees.
* @param {Array<number>} b Position as [lon, lat] in degrees.
* @returns {number} Initial bearing in degrees, clockwise from north.
*/
function bearing(a, b) {
var p1 = a[1] * RAD, p2 = b[1] * RAD, dl = (b[0] - a[0]) * RAD;
var y = Math.sin(dl) * Math.cos(p2);
var x = Math.cos(p1) * Math.sin(p2) - Math.sin(p1) * Math.cos(p2) * Math.cos(dl);
return (Math.atan2(y, x) * DEG + 360) % 360;
}
/**
* Position at distance r (degrees) from a point on course b (degrees).
*
* Computed in vectors rather than from the spherical law of cosines: the
* closed form degenerates to atan2(0, 0) exactly at the poles, which
* collapses a polar small circle into a single point. Here the tangent
* frame is built explicitly, and where it is undefined -- only at the two
* poles -- the reference direction is fixed to match the limit of the
* closed form, so nothing shifts away from the poles.
*
* @param {number} lon Start longitude in degrees.
* @param {number} lat Start latitude in degrees.
* @param {number} r Angular distance in degrees.
* @param {number} b Initial bearing in degrees, clockwise from north.
* @returns {Array<number>} The arrival point as [lon, lat] in degrees.
* @alias JXG.Geography.destination
*/
function destination(lon, lat, r, b) {
var c = toVector(lon, lat), n, e, ln, rr = r * RAD, bb = b * RAD, cr, sr;
n = [-c[0] * c[2], -c[1] * c[2], 1 - c[2] * c[2]]; // north, tangent
ln = Math.hypot(n[0], n[1], n[2]);
if (ln < 1e-12) { n = c[2] > 0 ? [-1, 0, 0] : [1, 0, 0]; ln = 1; }
n = [n[0] / ln, n[1] / ln, n[2] / ln];
e = [n[1] * c[2] - n[2] * c[1], // east = north x centre
n[2] * c[0] - n[0] * c[2],
n[0] * c[1] - n[1] * c[0]];
cr = Math.cos(rr); sr = Math.sin(rr);
return toGeo([cr * c[0] + sr * (Math.cos(bb) * n[0] + Math.sin(bb) * e[0]),
cr * c[1] + sr * (Math.cos(bb) * n[1] + Math.sin(bb) * e[1]),
cr * c[2] + sr * (Math.cos(bb) * n[2] + Math.sin(bb) * e[2])]);
}
/**
* Great-circle interpolation. t = 0 gives a, t = 1 gives b.
*
* @param {Array<number>} a Position as [lon, lat] in degrees.
* @param {Array<number>} b Position as [lon, lat] in degrees.
* @param {number} t 0 at a, 1 at b.
* @returns {Array<number>} The point on the great circle, as [lon, lat].
*/
function interpolate(a, b, t) {
var u = toVector(a[0], a[1]), v = toVector(b[0], b[1]);
var w = clamp(dot(u, v), -1, 1), o = Math.acos(w);
if (o < 1e-12) { return [a[0], a[1]]; }
var s = Math.sin(o), f1 = Math.sin((1 - t) * o) / s, f2 = Math.sin(t * o) / s;
return toGeo([f1 * u[0] + f2 * v[0], f1 * u[1] + f2 * v[1], f1 * u[2] + f2 * v[2]]);
}
/**
* The two quantities every orientation question needs, in one pass.
*
* `turn` is the total change in longitude with each step normalised to
* ±180, so ±360 means the ring circles a pole. `s` is the trapezoid
* integral of dLambda · mean(sin phi): the signed area of the band between
* the ring and the equator.
*
* Four functions used to walk the ring separately for these, with the same
* eight lines copied each time. That duplication is why the sign convention
* came out wrong in three different places before it came out right.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {{turn: number, s: number}} Degrees of longitude, and steradians.
*/
function ringSums(ring) {
var n = ring.length, s = 0, turn = 0, i, j, dl;
if (n < 3) { return { turn: 0, s: 0 }; }
for (i = 0; i < n; i++) {
j = (i + 1) % n;
dl = ring[j][0] - ring[i][0];
while (dl > 180) { dl -= 360; }
while (dl < -180) { dl += 360; }
turn += dl;
s += dl * RAD * (Math.sin(ring[i][1] * RAD) + Math.sin(ring[j][1] * RAD)) / 2;
}
return { turn: turn, s: s };
}
/**
* Can this ring's orientation be read at all?
*
* Not for one that circles a pole — the winding says which pole, not which
* way round — and not for one covering half the sphere, where the integral
* comes out at zero. A 90-degree cap is both.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {boolean} True when `orientationOf` can answer for this ring.
* @alias JXG.Geography.readableOrientation
*/
function readableOrientation(ring) {
var q = ringSums(ring);
return Math.abs(q.turn) <= 180 && Math.abs(q.s) > 1e-6;
}
/**
* Which way a ring runs: +1 counter-clockwise seen from outside, -1 the
* other way.
*
* A ring that circles a pole needs a different reading from one that does
* not. The trapezoid integral measures the complement in that case and its
* sign flips with it, so the winding decides instead — and the winding is
* unambiguous exactly when the integral is not.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {number} +1 or -1. Meaningless where `readableOrientation` is
* false; such rings keep the orientation they were given.
* @alias JXG.Geography.orientationOf
*/
function orientationOf(ring) {
var q = ringSums(ring);
if (Math.abs(q.turn) > 180) { return q.turn > 0 ? 1 : -1; }
return -q.s >= 0 ? 1 : -1;
}
/**
* Signed spherical area: positive for a ring wound counter-clockwise seen
* from outside, negative for a hole.
*
* `area` returns the magnitude, which is what a measurement wants. Every
* orientation question needs the sign, and taking it from the planar
* shoelace fails for exactly the rings that matter — those near a pole or
* across the antimeridian.
*
* No complement correction, unlike `area`: that correction gives the
* magnitude a pole-circling ring encloses but flips the sign with it, and
* the sign is the whole point. A cap larger than a hemisphere still runs
* counter-clockwise and is still an outer ring. The integral runs the other
* way round from the winding convention, so it is negated.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {number} Steradians, signed. A pole-circling ring reports the
* band it bounds rather than the cap.
* @alias JXG.Geography.signedArea
*/
function signedArea(ring) {
return -ringSums(ring).s;
}
/**
* Area of a spherical polygon in steradians.
*
* Sum of dLambda * mean(sin phi) integrates the band between the ring and
* the equator. A ring that encircles a pole therefore measures that band,
* not the cap, and the cap is its complement in the hemisphere. Densified
* edges make the trapezoid rule accurate enough; the difference to true
* great-circle edges vanishes below the segment length.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {number} Steradians, always positive. Multiply by R squared for
* an area on a sphere of that radius.
* @alias JXG.Geography.area
*/
function area(ring) {
var q = ringSums(ring), a2 = Math.abs(q.s);
if (Math.abs(q.turn) > 180) { a2 = 2 * Math.PI - a2; }
return a2;
}
/**
* Circle of constant angular radius r (degrees) around a position.
*
* Traversed with decreasing bearing, so the ring runs counter-clockwise as
* seen from outside the sphere — the same winding as every ring in a
* dataset. Increasing bearing (N-E-S-W) would be clockwise, and the
* terminator rejoin picks its arc direction from that winding: a reversed
* ring comes back as the complement of the intended area.
*
* @param {number} lon Centre longitude in degrees.
* @param {number} lat Centre latitude in degrees.
* @param {number} r Angular radius in degrees.
* @param {number} [n] Number of samples.
* @returns {Array<Array<number>>} An open ring, counter-clockwise about
* the centre seen from outside.
*/
function smallCircle(lon, lat, r, n) {
var out = [], i;
n = n || 64;
for (i = 0; i <= n; i++) { out.push(destination(lon, lat, r, -360 * i / n)); }
return out;
}
/**
* How many segments an arc needs so its chord stays within `tol` pixels.
*
* A chord across an angle d on a circle of screen radius R misses the arc
* by the sagitta R(1 - cos(d/2)). Setting that below the tolerance and
* solving for the number of segments gives
*
* n >= (span / 2) * sqrt(R / (2 * tol))
*
* so the count follows from how large the thing is drawn, not from a number
* chosen because it looked about right.
*
* @param {number} spanRad Angular span of the arc, in radians.
* @param {number} radiusPx Drawn radius in device pixels.
* @param {number} tolPx Allowed chord error in device pixels.
* @param {number} [min] Never return fewer than this.
* @param {number} [max] Never return more than this.
* @returns {number} Number of samples.
* @alias JXG.Geography.samplesFor
*/
function samplesFor(spanRad, radiusPx, tolPx, min, max) {
var lo = min || 8, hi = max || 512;
var tol = tolPx > 0 ? tolPx : 0.25; // device pixels
var R = radiusPx > 0 ? radiusPx : 1;
var n = Math.ceil(Math.abs(spanRad) / 2 * Math.sqrt(R / (2 * tol)));
// Math.max and Math.min pass NaN straight through, so a bad span would
// become a NaN point count and every curve built from it would vanish.
if (!isFinite(n)) { return lo; }
return Math.max(lo, Math.min(hi, n));
}
G.samplesFor = samplesFor;
/**
* Largest step that keeps the chord within tolerance, for a full circle.
*
* @param {number} radiusPx Drawn radius in device pixels.
* @param {number} tolPx Allowed chord error in device pixels.
* @returns {number} Step in degrees.
* @alias JXG.Geography.stepFor
*/
function stepFor(radiusPx, tolPx) {
var tol = tolPx > 0 ? tolPx : 0.25;
var R = radiusPx > 0 ? radiusPx : 1;
return 2 * Math.acos(clamp(1 - tol / R, -1, 1));
}
G.stepFor = stepFor;
/**
* Device pixels per CSS pixel.
*
* Board coordinates are in CSS pixels, but a corner is visible at the
* resolution the screen actually has: half a CSS pixel is a whole device
* pixel on a retina display, which is exactly where the remaining kinks
* came from. Sampling therefore works in device pixels.
* @alias JXG.Geography.devicePixels
*/
function devicePixels() {
var r = (typeof root !== 'undefined' && root.devicePixelRatio) || 1;
return r > 0 ? r : 1;
}
G.devicePixels = devicePixels;
/**
* Sampling count found by measuring the curve instead of assuming a circle.
*
* What is measured here is how far the curve departs from a great circle:
* for an evenly sampled great circle a point lies exactly on the bisector
* of its neighbours, so the measure is zero. It therefore does not replace
* samplesFor — which covers the chord error of the arc itself — but adds to
* it wherever a curve bends more tightly than a circle of the same screen
* size, as a ground track does at its turning latitudes.
*
* `make(n)` must return an array of [lon, lat]. Called a handful of times
* when the curve is built, never while drawing.
*
* @param {function(number): Array} make Builds the curve at n samples.
* @param {number} n0 Starting sample count.
* @param {number} radiusPx Drawn radius in device pixels.
* @param {number} tolPx Allowed chord error in device pixels.
* @param {number} [maxN] Upper bound on samples.
* @returns {Array} The curve, sampled finely enough.
* @alias JXG.Geography.refineSamples
*/
function refineSamples(make, n0, radiusPx, tolPx, maxN) {
var n = Math.max(4, n0), tol = tolPx > 0 ? tolPx : 0.25;
var R = radiusPx > 0 ? radiusPx : 1, cap = maxN || 2048, pts, i, sag;
while (n <= cap) {
pts = make(n);
sag = 0;
for (i = 1; i < pts.length - 1; i++) {
// deviation of the middle point from the chord, as an angle
var a = toVector(pts[i - 1][0], pts[i - 1][1]);
var b = toVector(pts[i][0], pts[i][1]);
var c = toVector(pts[i + 1][0], pts[i + 1][1]);
var m = [(a[0] + c[0]) / 2, (a[1] + c[1]) / 2, (a[2] + c[2]) / 2];
var L = Math.hypot(m[0], m[1], m[2]);
if (L < 1e-12) { continue; }
sag = Math.max(sag, Math.acos(clamp((b[0] * m[0] + b[1] * m[1] + b[2] * m[2]) / L, -1, 1)));
}
if (sag * R <= tol || n * 2 > cap) { return n; }
n *= 2;
}
return n;
}
G.refineSamples = refineSamples;
G.toVector = toVector; G.toGeo = toGeo; G.dot = dot;
G.angle = angle; G.distance = distance; G.bearing = bearing;
G.destination = destination; G.interpolate = interpolate;
G.area = area; G.smallCircle = smallCircle;
// ------------------------------------------------------------- sun / time
var OBLIQUITY = 23.4365; // epoch 2025
G.OBLIQUITY = OBLIQUITY;
/**
* Subsolar point for a date. Ignores the equation of time.
*
* @param {Date} [date] When; defaults to now.
* @returns {Array<number>} The point the sun is overhead, as [lon, lat]
* in degrees. The equation of time is not applied, so the longitude
* can be 4.1 degrees out — 457 km at the equator, around 31 October.
* That is fine for a terminator on a world map and not fine for
* anything that needs the sun's true position.
* @alias JXG.Geography.subsolar
*/
function subsolar(date) {
var d = date || new Date();
var start = Date.UTC(d.getUTCFullYear(), 0, 0);
var day = (Date.UTC(d.getUTCFullYear(), d.getUTCMonth(), d.getUTCDate()) - start) / 86400000;
var g = TAU / 365.24 * (day + 10);
var dec = -OBLIQUITY * Math.cos(g + 2 * 0.0167 * Math.sin(TAU / 365.24 * (day - 2)));
var utc = d.getUTCHours() + d.getUTCMinutes() / 60 + d.getUTCSeconds() / 3600;
var lon = 180 - utc * 15;
while (lon > 180) { lon -= 360; }
while (lon < -180) { lon += 360; }
return [lon, dec];
}
function terminator(date, n) {
var s = subsolar(date);
return smallCircle(s[0], s[1], 90, n || 240);
}
G.subsolar = subsolar; G.terminator = terminator;
// -------------------------------------------------------------- rotation
/**
* Rotation that brings (lon0, lat0) to the map centre (0, 0).
* For lat0 = 0 this collapses to a shift of the central meridian, which
* every projection supports; lat0 != 0 needs an oblique-capable one.
*
* @param {number} lon0 Longitude of the new centre, in degrees.
* @param {number} lat0 Latitude of the new centre, in degrees.
* @returns {function(number, number): Array<number>} lon, lat in
* degrees to the rotated [lon, lat], also in degrees.
* @alias JXG.Geography.rotator
*/
function rotator(lon0, lat0) {
var cl, sl, cp, sp;
if (lat0 === 0) {
return function (lon, lat) {
var l = lon - lon0;
while (l > 180) { l -= 360; }
while (l < -180) { l += 360; }
return [l, lat];
};
}
cl = Math.cos(-lon0 * RAD); sl = Math.sin(-lon0 * RAD);
cp = Math.cos(lat0 * RAD); sp = Math.sin(lat0 * RAD);
return function (lon, lat) {
var l = lon * RAD, p = lat * RAD, c = Math.cos(p);
var x = c * Math.cos(l), y = c * Math.sin(l), z = Math.sin(p);
var x1 = x * cl - y * sl, y1 = x * sl + y * cl; // Rz(-lon0)
var x2 = x1 * cp + z * sp, z2 = -x1 * sp + z * cp; // Ry(lat0)
return [Math.atan2(y1, x2) * DEG, Math.asin(clamp(z2, -1, 1)) * DEG];
};
}
function unrotator(lon0, lat0) {
var cl, sl, cp, sp;
if (lat0 === 0) {
return function (lon, lat) {
var l = lon + lon0;
while (l > 180) { l -= 360; }
while (l < -180) { l += 360; }
return [l, lat];
};
}
cl = Math.cos(lon0 * RAD); sl = Math.sin(lon0 * RAD);
cp = Math.cos(-lat0 * RAD); sp = Math.sin(-lat0 * RAD);
return function (lon, lat) {
var l = lon * RAD, p = lat * RAD, c = Math.cos(p);
var x = c * Math.cos(l), y = c * Math.sin(l), z = Math.sin(p);
var x1 = x * cp + z * sp, z1 = -x * sp + z * cp; // Ry(-lat0)
var x2 = x1 * cl - y * sl, y2 = x1 * sl + y * cl; // Rz(lon0)
return [Math.atan2(y2, x2) * DEG, Math.asin(clamp(z1, -1, 1)) * DEG];
};
}
/**
* Rotation that carries the map centre to the **north pole**.
*
* The cylindrical and pseudocylindrical projections put the centre on the
* equator, where the map edge is the meridian λ′ = ±180 — a straight cut.
* An azimuthal projection is built around a point instead: everything is
* expressed as an angular distance ρ from the centre and an azimuth θ. With
* the centre at the pole, ρ = 90° − φ′ and θ = λ′, so the map edge becomes
* the parallel φ′ = 90° − ρmax: a circle, and the reason this needs its own
* clipping.
*
* @param {number} lon0 Longitude of the new centre, in degrees.
* @param {number} lat0 Latitude of the new centre, in degrees.
* @returns {function(number, number): Array<number>} lon, lat in
* degrees to the rotated [lon, lat], also in degrees.
* @alias JXG.Geography.poleRotator
*/
function poleRotator(lon0, lat0) {
var cl = Math.cos(-lon0 * RAD), sl = Math.sin(-lon0 * RAD);
var a = (90 - lat0) * RAD, ca = Math.cos(a), sa = Math.sin(a);
return function (lon, lat) {
var l = lon * RAD, p = lat * RAD, c = Math.cos(p);
var x = c * Math.cos(l), y = c * Math.sin(l), z = Math.sin(p);
var x1 = x * cl - y * sl, y1 = x * sl + y * cl; // Rz(-lon0)
var x2 = x1 * ca - z * sa, z2 = x1 * sa + z * ca; // Ry(90 - lat0)
return [Math.atan2(y1, x2) * DEG, Math.asin(clamp(z2, -1, 1)) * DEG];
};
}
/**
* Inverse of poleRotator.
*
* @param {number} lon0 Longitude the centre was carried from.
* @param {number} lat0 Latitude the centre was carried from.
* @returns {function(number, number): Array<number>} Undoes
* `poleRotator` for the same centre.
*/
function poleUnrotator(lon0, lat0) {
var cl = Math.cos(lon0 * RAD), sl = Math.sin(lon0 * RAD);
var a = -(90 - lat0) * RAD, ca = Math.cos(a), sa = Math.sin(a);
return function (lon, lat) {
var l = lon * RAD, p = lat * RAD, c = Math.cos(p);
var x = c * Math.cos(l), y = c * Math.sin(l), z = Math.sin(p);
var x2 = x * ca - z * sa, z2 = x * sa + z * ca; // Ry(-(90 - lat0))
var x1 = x2 * cl - y * sl, y1 = x2 * sl + y * cl; // Rz(lon0)
return [Math.atan2(y1, x1) * DEG, Math.asin(clamp(z2, -1, 1)) * DEG];
};
}
G.signedArea = signedArea;
G.orientationOf = orientationOf;
G.readableOrientation = readableOrientation;
G.poleRotator = poleRotator; G.poleUnrotator = poleUnrotator;
/**
* More points where the azimuth turns fast.
*
* The source ring is sampled finely enough on the sphere, but an azimuthal
* projection maps the antipode of the centre onto the whole rim: a short arc
* passing near it is stretched around a large part of the circle. Segments
* are therefore split — along the great circle, so the inserted points are
* real positions — until no step turns more than maxTurn degrees of azimuth.
*
* @param {Array<Array<number>>} pts Points in the rotated frame.
* @param {number} maxTurn Largest step in azimuth, in degrees.
* @returns {Array<Array<number>>} The same ring with points inserted
* along the great circle, so they are real positions.
*/
function densifyAzimuth(pts, maxTurn) {
var out = [], i, a, b, d, k, n;
for (i = 0; i < pts.length; i++) {
a = pts[i];
b = pts[(i + 1) % pts.length];
out.push(a);
d = Math.abs(((b[0] - a[0] + 540) % 360) - 180);
if (d <= maxTurn) { continue; }
n = Math.min(64, Math.ceil(d / maxTurn));
for (k = 1; k < n; k++) { out.push(interpolate(a, b, k / n)); }
}
return out;
}
/**
* Filled ring clipped to a disc of angular radius rhoMax about the centre.
*
* Works in pole-rotated coordinates, where the disc is everything above the
* parallel φ′ = 90 − rhoMax. Runs of points inside are kept; where a run
* leaves and the next one enters, the two are joined along that parallel, in
* the direction the ring is wound. That join is the same idea as the
* terminator rejoin on the globe, and it is what makes a coastline running
* off the edge close against the rim instead of across the middle.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @param {function(number, number): Array<number>} rot Carries a
* position into the frame the cut is made in.
* @param {number} rhoMax Angular radius of the disc, in degrees.
* @param {number} [arcStep] Step along the rim when joining, in degrees.
* @returns {Array<Array<Array<number>>>} The pieces inside the disc.
* @alias JXG.Geography.clipRingCircle
*/
function clipRingCircle(ring, rot, rhoMax, arcStep) {
var phiC = 90 - rhoMax, step = arcStep || 2;
var n = ring.length, rr = new Array(n), inside = new Array(n), i;
var anyIn = false, allIn = true;
for (i = 0; i < n; i++) {
rr[i] = rot(ring[i][0], ring[i][1]);
inside[i] = rr[i][1] > phiC;
if (inside[i]) { anyIn = true; } else { allIn = false; }
}
if (allIn) { return [densifyAzimuth(rr, 3)]; }
if (!anyIn) {
// Nothing of the ring is on the map. Either it misses the disc entirely,
// or it surrounds it — an ocean ring around a small polar map. The two
// are told apart by the winding of the ring about the centre.
return enclosesCentre(rr, ring)
? [boundaryRing(phiC, ringSign(ring), step)] : [];
}
/**
* Where the segment a→b crosses the parallel, by latitude. When the
* azimuth between the two is not meaningful — the antipodal case — the
* nearer end's azimuth is used instead of an interpolation that would
* land anywhere on the rim.
*
* @param {Array<number>} a A point in the rotated frame.
* @param {Array<number>} b The next point, likewise.
* @returns {Array<number>} The crossing, on the boundary parallel.
*/
function cross(a, b) {
var t = (phiC - a[1]) / (b[1] - a[1]);
if (!usable(a, b)) {
return [t < 0.5 ? a[0] : b[0], phiC];
}
var la = a[0], lb = b[0];
while (lb - la > 180) { lb -= 360; }
while (la - lb > 180) { lb += 360; }
return [wrapLon(la + (lb - la) * t), phiC];
}
var runs = [], j, idx, E, X;
for (i = 0; i < n; i++) {
if (!(inside[i] && !inside[(i - 1 + n) % n])) { continue; }
idx = [i];
j = i;
while (inside[(j + 1) % n] && (j + 1) % n !== i) { j = (j + 1) % n; idx.push(j); }
E = cross(rr[(i - 1 + n) % n], rr[i]);
X = cross(rr[idx[idx.length - 1]], rr[(idx[idx.length - 1] + 1) % n]);
runs.push({ idx: idx, E: E, X: X });
}
if (!runs.length) { return []; }
// The direction to walk the rim comes from the ring's own orientation,
// and that has to be measured on the sphere: the planar shoelace fails
// for rings near a pole or across the antimeridian, which is exactly the
// set of rings a polar map is made of.
var dir = orientationOf(ring);
var out = [], used = new Array(runs.length), start, cur, guard, r, k, best, nxt, dd;
for (start = 0; start < runs.length; start++) {
if (used[start]) { continue; }
var poly = [];
cur = start;
guard = 0;
while (guard++ <= runs.length) {
used[cur] = true;
r = runs[cur];
poly.push(r.E);
for (k = 0; k < r.idx.length; k++) { poly.push(rr[r.idx[k]]); }
poly.push(r.X);
nxt = start;
best = Infinity;
for (k = 0; k < runs.length; k++) {
dd = ((runs[k].E[0] - r.X[0]) * dir) % 360;
if (dd < 0) { dd += 360; }
if (dd < best) { best = dd; nxt = k; }
}
// walk the rim from the exit to the next entry
var steps = Math.max(1, Math.round(best / step));
for (k = 1; k < steps; k++) {
poly.push([wrapLon(r.X[0] + dir * best * k / steps), phiC]);
}
if (nxt === start) { break; }
cur = nxt;
}
// A run that enters and leaves at the same place leaves a sliver of no
// area behind; it draws nothing and only muddles any orientation check.
if (poly.length > 2 && Math.abs(signedArea(poly)) > 1e-12) {
// The chaining can close a polygon the other way round — it depends on
// which rim arc was taken. The orientation is therefore imposed at the
// end: an outer ring keeps the sign of its spherical area, a hole the
// opposite, and the fill comes out the right way whatever route the
// walk took.
if (orientationOf(poly) !== orientationOf(ring)) { poly.reverse(); }
out.push(densifyAzimuth(poly, 3));
}
}
return out;
}
/**
* The rim itself, as a ring wound the given way.
*
* @param {number} phiC The boundary parallel, in degrees.
* @param {number} sign +1 walks it eastward, -1 westward.
* @param {number} [step] Step along the rim, in degrees.
* @returns {Array<Array<number>>} The rim as an open ring.
*/
function boundaryRing(phiC, sign, step) {
var out = [], i, n = Math.max(24, Math.round(360 / (step || 2)));
for (i = 0; i < n; i++) {
out.push([wrapLon(sign >= 0 ? -180 + 360 * i / n : 180 - 360 * i / n), phiC]);
}
return out;
}
/**
* Does the ring contain the map centre — now the north pole?
*
* Winding alone cannot say: a ring around the antipode winds a full turn
* too, only the other way. The direction that settles it is the ring's own
* orientation, and that has to be measured on the sphere. The planar
* shoelace is unreliable for exactly the rings in question, which reach
* across the antimeridian or around a pole.
*
* @param {Array<Array<number>>} rr The ring in the rotated frame.
* @param {Array<Array<number>>} ring The same ring, unrotated. Its
* orientation is what tells a surrounding ring from a far-side one.
* @returns {boolean} Whether the ring surrounds the centre of the disc.
*/
function enclosesCentre(rr, ring) {
return windingAbout(rr) * orientationOf(ring) > 180;
}
function windingAbout(rr) {
var total = 0, i, d;
for (i = 0; i < rr.length; i++) {
d = rr[(i + 1) % rr.length][0] - rr[i][0];
while (d > 180) { d -= 360; }
while (d < -180) { d += 360; }
total += d;
}
return total;
}
function wrapLon(l) {
while (l > 180) { l -= 360; }
while (l < -180) { l += 360; }
return l;
}
/**
* Is the azimuth between two neighbouring points meaningful?
*
* Near the antipode of the centre it is not: two points a degree apart on
* the ground can sit half a turn apart in θ, and anything derived from
* interpolating between them — a crossing point, a straight join — lands
* somewhere arbitrary on the rim and draws a chord across the disc. A large
* step in θ over a short arc is the giveaway.
*
* @param {Array<number>} a A point as [lon, lat] in the rotated frame.
* @param {Array<number>} b The neighbouring point, likewise.
* @returns {boolean} False when the step in azimuth is too large to be real.
*/
function usable(a, b) {
return !(angle(a, b) < 5 &&
Math.abs(((b[0] - a[0] + 540) % 360) - 180) > 20);
}
/**
* Open polyline clipped to the disc; no rim arcs, just the pieces inside.
*
* @param {Array<Array<number>>} pts Open polyline of [lon, lat] in degrees.
* @param {function(number, number): Array<number>} rot Carries a position
* into the frame the cut is made in.
* @param {number} rhoMax Angular radius of the disc, in degrees.
* @returns {Array<Array<Array<number>>>} The pieces inside the disc.
* @alias JXG.Geography.clipLineCircle
*/
function clipLineCircle(pts, rot, rhoMax) {
var phiC = 90 - rhoMax, out = [], run = [], i, b, t, prev = null;
for (i = 0; i < pts.length; i++) {
var q = rot(pts[i][0], pts[i][1]);
if (prev && run.length && !usable(prev, q)) {
if (run.length > 1) { out.push(run); }
run = [];
}
if (q[1] > phiC) {
if (prev && prev[1] <= phiC && usable(prev, q)) {
t = (phiC - prev[1]) / (q[1] - prev[1]);
b = prev[0];
while (q[0] - b > 180) { b += 360; }
while (b - q[0] > 180) { b -= 360; }
run.push([wrapLon(b + (q[0] - b) * t), phiC]);
}
run.push(q);
} else if (prev && prev[1] > phiC) {
if (usable(prev, q)) {
t = (phiC - prev[1]) / (q[1] - prev[1]);
b = q[0];
while (b - prev[0] > 180) { b -= 360; }
while (prev[0] - b > 180) { b += 360; }
run.push([wrapLon(prev[0] + (b - prev[0]) * t), phiC]);
}
if (run.length > 1) { out.push(run); }
run = [];
}
prev = q;
}
if (run.length > 1) { out.push(run); }
return out;
}
G.clipRingCircle = clipRingCircle;
G.clipLineCircle = clipLineCircle;
G.rotator = rotator; G.unrotator = unrotator;
/**
* Map centre from a camera rotation matrix.
*
* lon = atan2(cam_y, cam_x) collapses at the poles: both components go to
* zero and the arctangent swings over its whole range on the slightest
* pointer move. Row 2 of the matrix is exactly the preimage of the screen
* up-axis, (-sin p cos l, -sin p sin l, cos p), and carries the longitude
* well-conditioned wherever the eye vector loses it. With zero roll both
* routes agree, so the switch at |lat| = 45° introduces no jump.
*
*
* @param {Array<Array<number>>} m A 3×3 rotation matrix.
* @returns {Array<number>} The centre it carries, as [lon, lat].
* @alias JXG.Geography.centreFromMatrix
*/
G.centreFromMatrix = function (m) {
var cam = [m[3][1], m[3][2], m[3][3]], r2 = [m[2][1], m[2][2], m[2][3]];
var lat = Math.asin(clamp(cam[2], -1, 1)) * DEG, lon, sg;
if (Math.hypot(cam[0], cam[1]) >= Math.abs(cam[2])) {
lon = Math.atan2(cam[1], cam[0]) * DEG;
} else {
sg = cam[2] >= 0 ? 1 : -1;
lon = Math.atan2(-sg * r2[1], -sg * r2[0]) * DEG;
}
return [lon, lat];
};
// -------------------------------------------------------------- clipping
/**
* Longitudes made continuous, so a ring can be cut against a meridian.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @returns {Array<Array<number>>} A copy whose longitudes run
* continuously, and may leave the range ±180.
* @alias JXG.Geography.unwrap
*/
function unwrap(ring) {
var out = [[ring[0][0], ring[0][1]]], i, lo, prev;
for (i = 1; i < ring.length; i++) {
lo = ring[i][0];
prev = out[i - 1][0];
while (lo - prev > 180) { lo -= 360; }
while (prev - lo > 180) { lo += 360; }
out.push([lo, ring[i][1]]);
}
return out;
}
/**
* Total longitude change around a ring; ±360 means it encloses a pole.
*
* @param {Array<Array<number>>} unw A ring with unwrapped longitudes.
* @returns {{lo: number, hi: number}} The range of turns it spans, so
* a ring can be cut once per turn.
*/
function winding(unw) {
var n = unw.length, lo = unw[0][0], prev = unw[n - 1][0];
while (lo - prev > 180) { lo -= 360; }
while (prev - lo > 180) { lo += 360; }
return unw[n - 1][0] - unw[0][0] + (lo - unw[n - 1][0]);
}
function clipHalf(poly, keepGreater, x) {
var res = [], n = poly.length, i, a, b, ia, ib, t;
function inside(p) { return keepGreater ? p[0] >= x : p[0] <= x; }
for (i = 0; i < n; i++) {
a = poly[i]; b = poly[(i + 1) % n];
ia = inside(a); ib = inside(b);
if (ia) { res.push(a); }
if (ia !== ib) {
t = (x - a[0]) / (b[0] - a[0]);
res.push([x, a[1] + t * (b[1] - a[1])]);
}
}
return res;
}
/**
* Twice the signed planar area of a ring, in its own coordinates.
*
* A planar measure on planar data. It is *not* a way to read the winding of
* a geographic ring: for one near a pole or across the antimeridian it
* answers confidently and wrongly. Use `orientationOf` for that.
*
* @param {Array<Array<number>>} r Open ring of [x, y].
* @returns {number} Twice the signed area; positive counter-clockwise.
* @alias JXG.Geography.shoelace
*/
function shoelace(r) {
var a = 0, n = r.length, i, j;
for (i = 0; i < n; i++) { j = (i + 1) % n; a += r[i][0] * r[j][1] - r[j][0] * r[i][1]; }
return a / 2;
}
/**
* Edges created by clipping and by pole closure are long and would be drawn
* as chords instead of following the curved map boundary.
*
* @param {Array<Array<number>>} poly Points as [lon, lat] in degrees.
* @param {number} maxStep Longest segment to leave alone, in degrees.
* @returns {Array<Array<number>>} The polygon with points inserted.
*/
function densify(poly, maxStep) {
var out = [], i, k, a, b, n, t;
for (i = 0; i < poly.length; i++) {
a = poly[i]; b = poly[(i + 1) % poly.length];
out.push(a);
n = Math.floor(Math.max(Math.abs(b[0] - a[0]), Math.abs(b[1] - a[1])) / maxStep);
for (k = 1; k <= n; k++) {
t = k / (n + 1);
out.push([a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]);
}
}
return out;
}
/**
* Filled ring: rotate, unwrap, clip into the strip [-180, 180]. Three shifts
* cover every wrapping case, and Sutherland–Hodgman closes each piece along
* the map edge by itself.
*
* @param {Array<Array<number>>} pts Points as [lon, lat].
* @param {number} dx Shift in degrees of longitude.
* @returns {Array<Array<number>>} A shifted copy.
*/
function shiftBy(pts, dx) {
var out = [], i;
for (i = 0; i < pts.length; i++) { out.push([pts[i][0] + dx, pts[i][1]]); }
return out;
}
/**
* Filled ring clipped to a sheet, split where it crosses the antimeridian.
*
* The ring is unwrapped, cut against ±180 in each turn it spans, and closed
* along the pole line where it circles a pole. The orientation is imposed
* afterwards, on the densified result: the winding is read from steps
* normalised to ±180, and a single step longer than that normalises the
* other way round, so measuring before the split and shipping after it gave
* the wrong answer for every ring reaching a pole.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @param {function(number, number): Array<number>} rot Carries a position
* into the frame the cut is made in.
* @param {number} maxStep Longest segment to leave undivided, in degrees.
* @returns {Array<Array<Array<number>>>} The pieces, each an open ring.
* @alias JXG.Geography.clipRing
*/
function clipRing(ring, rot, maxStep) {
var unw = unwrap(ring.map(function (p) { return rot(p[0], p[1]); }));
var w = winding(unw), out = [], k, shifted, poly, pole, a, b;
if (Math.abs(w) > 1) {
// Encloses a rotated pole, so the unwrapped run is open. Close it over
// the pole line: with counter-clockwise rings the interior is on the
// left, hence positive winding means the north pole.
pole = w > 0 ? 90 : -90;
a = unw[unw.length - 1]; b = unw[0];
unw = unw.concat([[a[0], pole], [b[0], pole]]);
}
for (k = -1; k <= 1; k++) {
// shift bound per turn: a closure over the loop variable would read
// whatever k happens to be when it is finally called
shifted = shiftBy(unw, 360 * k);
poly = clipHalf(shifted, true, -180);
if (poly.length < 3) { continue; }
poly = clipHalf(poly, false, 180);
if (poly.length >= 3 && Math.abs(shoelace(poly)) > 1e-4) {
// Densify first, then decide. The winding is read from normalised
// steps, and a single step longer than 180° normalises the other way:
// splitting the pole-line closure into short pieces therefore changes
// the answer. Measuring before the split and shipping after it was
// wrong in exactly the cases that reach a pole.
poly = densify(poly, maxStep);
// Closing a ring against the map edge can reverse it, so the source
// orientation is imposed — but only where it can be read. A ring that
// circles a pole, or one whose integral vanishes because it covers
// half the sphere, has no readable orientation, and guessing at it
// would turn a night cap inside out. Those keep what they came with.
if (readableOrientation(ring) &&
orientationOf(poly) !== orientationOf(ring)) {
poly.reverse();
}
out.push(poly);
}
}
return out;
}
/**
* Open polyline: split where it leaves the map, inserting the crossing.
*
* Longitudes are unwrapped first, the way filled rings already are. Reading
* a jump of more than 180 degrees as an antimeridian crossing fails near a
* rotated pole, where a small step on the ground is a large step in
* longitude: the line then gets a spurious break and runs off to the map
* edge and back — the stray diagonals across an obliquely centred map.
*
* @param {Array<Array<number>>} pts Open polyline of [lon, lat].
* @param {function(number, number): Array<number>} rot Carries a
* position into the frame the cut is made in.
* @returns {Array<Array<Array<number>>>} The pieces.
*/
function clipLine(pts, rot) {
var rr = pts.map(function (p) { return rot(p[0], p[1]); });
if (rr.length < 2) { return []; }
var unw = unwrap(rr), out = [], cur = [unw[0]],
i, a, b, k, lo, hi, edge, t, lat, guard;
for (i = 1; i < unw.length; i++) {
a = cur[cur.length - 1];
b = unw[i];
// Every odd multiple of 180 between the two is a map edge, and a long
// segment can cross several. The bounds are recomputed each time round,
// because the remaining part of the segment starts at the last crossing.
guard = 0;
while (guard++ < 16) {
lo = Math.min(a[0], b[0]);
hi = Math.max(a[0], b[0]);
k = Math.floor((lo - 180) / 360) + 1;
edge = k * 360 + 180;
if (!(edge > lo && edge < hi) || b[0] === a[0]) { break; }
t = (edge - a[0]) / (b[0] - a[0]);
lat = a[1] + (b[1] - a[1]) * t;
cur.push([edge, lat]);
out.push(cur);
// the next piece starts at the same crossing and carries on in the
// unwrapped frame; folding happens once, at the end, per piece
cur = [[edge, lat]];
a = [edge, lat];
}
cur.push(b);
}
if (cur.length > 1) { out.push(cur); }
/**
* Fold a whole piece back into [-180, 180] with one shift. Clamping point
* by point would tear the piece apart again: two neighbours either side of
* the edge would land 360 degrees from each other.
*/
var folded = out.map(function (piece) {
var min = Infinity, max = -Infinity, shift, j;
for (j = 0; j < piece.length; j++) {
if (piece[j][0] < min) { min = piece[j][0]; }
if (piece[j][0] > max) { max = piece[j][0]; }
}
shift = -360 * Math.round((min + max) / 720);
if (min + shift < -180.001) { shift += 360; }
if (max + shift > 180.001) { shift -= 360; }
return piece.map(function (q) { return [q[0] + shift, q[1]]; });
});
// A last pass along each folded piece. Its points now run smoothly, so a
// longitude beyond the frame really is a crossing and can be cut there —
// unlike at the start, where the same test misfired near a rotated pole.
var done = [];
folded.forEach(function (piece) {
var run = [], j, q, prev, tt, la;
for (j = 0; j < piece.length; j++) {
q = piece[j];
if (Math.abs(q[0]) <= 180.0001) {
if (run.length === 0 && j > 0) {
prev = piece[j - 1];
tt = ((prev[0] > 180 ? 180 : -180) - prev[0]) / (q[0] - prev[0]);
la = prev[1] + (q[1] - prev[1]) * tt;
run.push([prev[0] > 180 ? 180 : -180, la]);
}
run.push(q);
} else {
if (run.length) {
prev = run[run.length - 1];
tt = ((q[0] > 180 ? 180 : -180) - prev[0]) / (q[0] - prev[0]);
la = prev[1] + (q[1] - prev[1]) * tt;
run.push([q[0] > 180 ? 180 : -180, la]);
if (run.length > 1) { done.push(run); }
run = [];
}
}
}
if (run.length > 1) { done.push(run); }
});
return done;
}
/**
* Concatenate fragments into one point list with NaN separators.
*
* @param {Array<Array<Array<number>>>} parts Pieces to be joined.
* @param {boolean} close Whether the result is a closed ring.
* @returns {Array<Array<Array<number>>>} The pieces, joined where
* their ends meet.
*/
function join(parts, close) {
var out = [], i, j;
for (i = 0; i < parts.length; i++) {
if (i) { out.push([NaN, NaN]); }
for (j = 0; j < parts[i].length; j++) { out.push(parts[i][j]); }
if (close && parts[i].length > 2) { out.push(parts[i][0]); }
}
return out;
}
/**
* Drop consecutive duplicates, and the closing repeat of the first point.
*
* @param {Array<Array<number>>} pts Points, possibly with repeats.
* @returns {Array<Array<number>>} The same points, repeats dropped.
*/
function dedupe(pts) {
if (!pts || !pts.length) { return []; }
var out = [pts[0]], i, p, q;
for (i = 1; i < pts.length; i++) {
p = out[out.length - 1]; q = pts[i];
if (Math.abs(p[0] - q[0]) > 1e-12 || Math.abs(p[1] - q[1]) > 1e-12 ||
Math.abs(p[2] - q[2]) > 1e-12) { out.push(q); }
}
while (out.length > 3) {
p = out[0]; q = out[out.length - 1];
if (Math.abs(p[0] - q[0]) < 1e-12 && Math.abs(p[1] - q[1]) < 1e-12 &&
Math.abs(p[2] - q[2]) < 1e-12) { out.pop(); } else { break; }
}
return out;
}
/**
* The visible stretches of a ring, each with the exact points where it
* enters and leaves the visible side and their angles on the silhouette.
* Kept apart from the rejoin so that neither half has to be read with the
* other in mind.
*
* @param {number} n Number of points in the ring.
* @param {Array<boolean>} m Visibility per point.
* @param {Array} cross Where visibility changes, per segment.
* @param {Array<number>} psi Angle along the limb, per crossing.
* @returns {Array<Object>} One entry per visible run.
*/
function visibleRuns(n, m, cross, psi) {
var runs = [], i, j, idx, E, X;
for (i = 0; i < n; i++) {
if (!(m[i] && !m[(i - 1 + n) % n])) { continue; }
idx = [i];
j = i;
while (m[(j + 1) % n] && (j + 1) % n !== i) { j = (j + 1) % n; idx.push(j); }
E = cross((i - 1 + n) % n, i);
X = cross(idx[idx.length - 1], (idx[idx.length - 1] + 1) % n);
runs.push({ idx: idx, E: E, X: X, psiE: psi(E), psiX: psi(X) });
}
return runs;
}
/**
* Terminator clipping for a filled ring on the sphere.
*
* Returns the parts of `pts` on the hemisphere `side` faces, each closed
* over arcs of the silhouette circle. Where a ring crosses the terminator
* more than twice the pieces must be chained along that circle: closing
* each piece on its own produces the complement of the intended area.
*
* r1, r2 are the screen axes in world coordinates (rows 1 and 2 of the
* camera rotation), cam is the eye vector (row 3). Points are unit vectors;
* the caller scales and projects.
*
* @param {Array<Array<number>>} pts Open ring of [lon, lat] in degrees.
* @param {Array<number>} r1 First row of the view rotation.
* @param {Array<number>} r2 Second row of the view rotation.
* @param {Array<number>} cam Camera direction as a unit vector.
* @param {number} side +1 keeps the hemisphere facing the camera, -1
* the one behind it. It is the sign the dot product with `cam` must
* have for a point to count as visible.
* @param {number} [arcStep] Step along the limb when rejoining.
* @param {number} [orient=1] Winding of the ring, +1 or -1. Stated by
* the caller rather than derived, because a great circle spans every
* longitude and no measure can read its orientation.
* @returns {Array<Array<Array<number>>>} The visible pieces.
* @alias JXG.Geography.clipSphereRing
*/
function clipSphereRing(pts, r1, r2, cam, side, arcStep, orient) {
// A ring handed in closed carries a zero-length segment, and a freehand
// stroke can repeat a point outright. Both make the crossing computation
// divide by zero, so they are removed before anything else happens.
if (!pts || pts.length < 3) { return []; }
pts = dedupe(pts);
if (pts.length < 3) { return []; }
var n = pts.length, d = new Array(n), m = new Array(n),
anyIn = false, allIn = true, i, out = [], runs, used, start, cur,
guard, r, nxt, best, dd, steps, k, a, ca, sa;
arcStep = arcStep || 3 * RAD;
// Which way the closing arc runs depends on the ring's own winding, not
// only on the hemisphere. Holes are stored clockwise so the nonzero rule
// punches them out; feeding one through the counter-clockwise assumption
// returns the complement — a fragment that swallows the whole disc and
// flips land against ocean for the entire layer.
var dir = side * (orient === undefined ? 1 : orient);
for (i = 0; i < n; i++) {
d[i] = dot(pts[i], cam);
m[i] = d[i] * side > 0;
if (m[i]) { anyIn = true; } else { allIn = false; }
}
if (!anyIn) { return []; }
if (allIn) { return [pts.concat([pts[0]])]; }
// Angular radius of the ring about its own centre: how much of the sphere
// it can possibly cover.
var mid = [0, 0, 0], extent = 0;
for (i = 0; i < n; i++) { mid[0] += pts[i][0]; mid[1] += pts[i][1]; mid[2] += pts[i][2]; }
var midL = Math.hypot(mid[0], mid[1], mid[2]);
if (midL > 1e-12) {
mid = [mid[0] / midL, mid[1] / midL, mid[2] / midL];
for (i = 0; i < n; i++) { extent = Math.max(extent, Math.acos(clamp(dot(pts[i], mid), -1, 1))); }
} else { extent = Math.PI; }
function cross(ia, ib) {
var t = d[ia] / (d[ia] - d[ib]), p = pts[ia], q = pts[ib];
var v = [p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1]), p[2] + t * (q[2] - p[2])];
var L = Math.hypot(v[0], v[1], v[2]);
return [v[0] / L, v[1] / L, v[2] / L];
}
function psi(q) { return Math.atan2(dot(q, r2), dot(q, r1)); }
runs = visibleRuns(n, m, cross, psi);
if (!runs.length) { return []; }
used = new Array(runs.length);
for (start = 0; start < runs.length; start++) {
if (used[start]) { continue; }
var poly = [];
cur = start; guard = 0;
while (guard++ <= runs.length) {
used[cur] = true;
r = runs[cur];
poly.push(r.E);
for (k = 0; k < r.idx.length; k++) { poly.push(pts[r.idx[k]]); }
poly.push(r.X);
nxt = start; best = Infinity;
for (k = 0; k < runs.length; k++) {
dd = ((runs[k].psiE - r.psiX) * dir) % TAU;
if (dd < 0) { dd += TAU; }
if (dd < best) { best = dd; nxt = k; }
}
// A ring reaching `extent` from its own centre can only ever need a
// closing arc of about twice that. Anything longer has wrapped the
// wrong way round and the complement is meant — this is the flash
// that turns land against ocean for a whole layer. Rings covering a
// large part of the sphere are left alone: there the long arc can be
// the right answer.
if (extent < Math.PI / 2 && best > Math.min(Math.PI, 2 * extent + 0.2)) {
best -= TAU;
}
steps = Math.max(1, Math.floor(Math.abs(best) / arcStep));
for (k = 1; k < steps; k++) {
a = r.psiX + dir * best * k / steps;
ca = Math.cos(a); sa = Math.sin(a);
poly.push([ca * r1[0] + sa * r2[0], ca * r1[1] + sa * r2[1], ca * r1[2] + sa * r2[2]]);
}
if (nxt === start || used[nxt]) { break; }
cur = nxt;
}
poly.push(poly[0]);
out.push(poly);
}
return out;
}
G.clipSphereRing = clipSphereRing;
/**
* Winding of a lon/lat ring: +1 counter-clockwise seen from outside the
* sphere, -1 clockwise. Holes are stored clockwise so the nonzero fill rule
* punches them out, and the terminator rejoin needs to know which it has.
*
* This is a planar measure in lon/lat and only agrees with the spherical
* one for rings that do not circle the sphere. A great circle spans every
* longitude, and its shoelace then reports the area under a sine wave
* rather than the side its interior lies on. Callers who know the winding
* by construction — a small circle, for instance — should say so instead
* of asking.
*
* @param {Array<Array<number>>} ll Open ring of [lon, lat] in degrees.
* @returns {number} +1 or -1. Planar, and therefore wrong for rings near
* a pole or across the antimeridian; `orientationOf` is the one to ask.
* @alias JXG.Geography.projections
* @alias JXG.Geography.ringSign
*/
function ringSign(ll) {
return shoelace(ll) < 0 ? -1 : 1;
}
G.ringSign = ringSign;
G.unwrap = unwrap; G.winding = winding; G.densify = densify;
G.clipRing = clipRing; G.clipLine = clipLine; G.join = join;
G.shoelace = shoelace;
// ----------------------------------------------------------- projections
G.projections = {};
/**
* Add a projection under a name.
*
* @param {string} name The key it is looked up by.
* @param {Object} desc Needs `forward(lon, lat)`; may carry `inverse(x, y)`,
* `latMax`, `oblique`, `clip` and `rhoMax`. Without an inverse a click
* cannot be turned back into a position, so `G.invert` falls back to
* searching, which is slower and can fail near the edges.
* @returns {Object} The projection as registered.
* @alias JXG.Geography.registerProjection
*/
G.registerProjection = function (name, desc) {
if (typeof desc.forward !== 'function') {
warn('projection "' + name + '" has no forward()');
return null;
}
desc.id = name;
desc.latMax = desc.latMax === undefined ? 90 : desc.latMax;
desc.oblique = !!desc.oblique;
desc.clip = desc.clip || 'antimeridian';
desc.kind = desc.kind || 'compromise';
G.projections[name] = desc;
return desc;
};
/**
* Look up a projection by name, or pass one through.
*
* @param {(string|Object)} p A registered name, or a projection itself.
* @returns {?Object} The projection, or null with a warning.
* @alias JXG.Geography.projection
*/
G.projection = function (p) {
if (typeof p !== 'string') { return p; }
if (!G.projections[p]) { warn('unknown projection "' + p + '"'); return null; }
return G.projections[p];
};
(function registerBuiltins() {
var E1 = 1.340264, E2 = -0.081106, E3 = 0.000893, E4 = 0.003796,
EM = Math.sqrt(3) / 2, WT1 = Math.acos(2 / Math.PI);
var ROB_X = [1.0000, 0.9986, 0.9954, 0.9900, 0.9822, 0.9730, 0.9600, 0.9427,
0.9216, 0.8962, 0.8679, 0.8350, 0.7986, 0.7597, 0.7186, 0.6732,
0.6213, 0.5722, 0.5322];
var ROB_Y = [0.0000, 0.0620, 0.1240, 0.1860, 0.2480, 0.3100, 0.3720, 0.4340,
0.4958, 0.5571, 0.6176, 0.6769, 0.7346, 0.7903, 0.8435, 0.8936,
0.9394, 0.9761, 1.0000];
function catmull(tab, t) {
var n = tab.length, i = clamp(Math.floor(t), 0, n - 2), u = t - i;
var p0 = tab[Math.max(0, i - 1)], p1 = tab[i],
p2 = tab[i + 1], p3 = tab[Math.min(n - 1, i + 2)];
return 0.5 * (2 * p1 + (p2 - p0) * u + (2 * p0 - 5 * p1 + 4 * p2 - p3) * u * u +
(3 * p1 - 3 * p2 + p3 - p0) * u * u * u);
}
/**
* Azimuthal projections.
*
* These are built around a point, not a line: `forward` receives
* pole-rotated coordinates, in which the map centre sits at the north pole
* and ρ = 90 − φ′ is the angular distance from it. Each one differs only
* in how ρ becomes a radius on the sheet, and each buys one property at
* the cost of the others:
*
* equidistant r = ρ distances from the centre are true
* equal-area r = 2 sin(ρ/2) areas are true, at any radius
* stereographic r = 2 tan(ρ/2) angles are true; the rim runs away
* orthographic r = sin ρ the view from infinitely far off
*
* `clip: 'circle'` tells the map to cut against a circle of `rhoMax`
* rather than against a meridian.
*/
/**
* Register one azimuthal projection.
*
* @param {string} id The key it is looked up by.
* @param {Object} names Display names by language.
* @param {string} kind What it preserves: equidistant, equalarea,
* conformal or perspective.
* @param {function(number): number} radiusOf Angular distance from the
* centre, in radians, to a radius on the sheet. Carries its own
* `invert`, so a click needs no search.
* @param {number} rhoMax How far from the centre the disc reaches, in
* degrees. Short of 180 for all of them: right at the antipode the
* azimuth has no meaning.
* @returns {void}
*/
function azimuthal(id, names, kind, radiusOf, rhoMax) {
G.registerProjection(id, {
name: names, kind: kind, oblique: true, azimuthal: true,
clip: 'circle', rhoMax: rhoMax, latMax: 90,
forward: function (lon, lat) {
// In the pole-rotated frame θ = 180° points north and θ = 90° east,
// so north is up and east is right only with y = −r cos θ. Getting
// this wrong mirrors the map without otherwise breaking anything.
var rho = (90 - lat) * RAD, th = lon * RAD, r = radiusOf(rho);
return [r * Math.sin(th), -r * Math.cos(th)];
},
inverse: function (x, y) {
var r = Math.hypot(x, y);
if (r < 1e-12) { return [0, 90]; }
return [Math.atan2(x, -y) * DEG, 90 - radiusOf.invert(r) * DEG];
}
});
}
function withInverse(f, g) { f.invert = g; return f; }
azimuthal('azimuthalequidistant',
{ de: 'Azimutal abstandstreu', en: 'Azimuthal equidistant' },
'equidistant',
withInverse(function (rho) { return rho; },
function (r) { return r; }), 150);
azimuthal('azimuthalequalarea',
{ de: 'Lambert azimutal', en: 'Lambert azimuthal' },
'equalarea',
withInverse(function (rho) { return 2 * Math.sin(rho / 2); },
function (r) { return 2 * Math.asin(clamp(r / 2, -1, 1)); }),
150);
azimuthal('stereographic',
{ de: 'Stereografisch', en: 'Stereographic' },
'conformal',
withInverse(function (rho) { return 2 * Math.tan(rho / 2); },
function (r) { return 2 * Math.atan(r / 2); }), 150);
azimuthal('orthographic',
{ de: 'Orthografisch', en: 'Orthographic' },
'perspective',
withInverse(function (rho) { return Math.sin(rho); },
function (r) { return Math.asin(clamp(r, -1, 1)); }), 90);
G.registerProjection('equalearth', {
name: { de: 'Equal Earth', en: 'Equal Earth' },
kind: 'equalarea', oblique: true,
forward: function (lon, lat) {
var l = lon * RAD, t = Math.asin(EM * Math.sin(lat * RAD));
var t2 = t * t, t6 = t2 * t2 * t2, t8 = t6 * t2;
var d = 9 * E4 * t8 + 7 * E3 * t6 + 3 * E2 * t2 + E1;
return [2 * Math.sqrt(3) * l * Math.cos(t) / (3 * d),
t * (E1 + E2 * t2 + E3 * t6 + E4 * t8)];
}
});
G.registerProjection('mollweide', {
name: { de: 'Mollweide', en: 'Mollweide' },
kind: 'equalarea', oblique: true,
forward: function (lon, lat) {
var l = lon * RAD, p = lat * RAD, t = p, i, f, fp, d;
if (Math.abs(Math.abs(p) - Math.PI / 2) > 1e-9) {
for (i = 0; i < 24; i++) {
f = 2 * t + Math.sin(2 * t) - Math.PI * Math.sin(p);
fp = 2 + 2 * Math.cos(2 * t);
if (Math.abs(fp) < 1e-12) { break; }
d = f / fp; t -= d;
if (Math.abs(d) < 1e-13) { break; }
}
} else { t = (p < 0 ? -1 : 1) * Math.PI / 2; }
return [2 * Math.SQRT2 / Math.PI * l * Math.cos(t), Math.SQRT2 * Math.sin(t)];
}
});
G.registerProjection('hammer', {
name: { de: 'Hammer', en: 'Hammer' },
kind: 'equalarea', oblique: true,
forward: function (lon, lat) {
var l = lon * RAD, p = lat * RAD;
var z = Math.sqrt(1 + Math.cos(p) * Math.cos(l / 2));
return [2 * Math.SQRT2 * Math.cos(p) * Math.sin(l / 2) / z,
Math.SQRT2 * Math.sin(p) / z];
}
});
G.registerProjection('sinusoidal', {
name: { de: 'Sinusoidal', en: 'Sinusoidal' },
kind: 'equalarea',
forward: function (lon, lat) {
var p = lat * RAD;
return [lon * RAD * Math.cos(p), p];
}
});
G.registerProjection('robinson', {
name: { de: 'Robinson', en: 'Robinson' },
kind: 'compromise',
forward: function (lon, lat) {
var a = Math.abs(lat) / 5;
return [0.8487 * catmull(ROB_X, a) * lon * RAD,
1.3523 * catmull(ROB_Y, a) * (lat < 0 ? -1 : 1)];
}
});
G.registerProjection('winkel', {
name: { de: 'Winkel Tripel', en: 'Winkel tripel' },
kind: 'compromise',
forward: function (lon, lat) {
var l = lon * RAD, p = lat * RAD;
var a = Math.acos(clamp(Math.cos(p) * Math.cos(l / 2), -1, 1));
var s = Math.abs(a) < 1e-10 ? 1 : Math.sin(a) / a;
return [0.5 * (l * Math.cos(WT1) + 2 * Math.cos(p) * Math.sin(l / 2) / s),
0.5 * (p + Math.sin(p) / s)];
}
});
G.registerProjection('mercator', {
name: { de: 'Mercator', en: 'Mercator' },
kind: 'conformal', latMax: 84,
forward: function (lon, lat) {
var p = clamp(lat, -84, 84) * RAD;
return [lon * RAD, Math.log(Math.tan(Math.PI / 4 + p / 2))];
}
});
G.registerProjection('platecarree', {
name: { de: 'Plattkarte', en: 'Plate carrée' },
kind: 'equidistant',
forward: function (lon, lat) { return [lon * RAD, lat * RAD]; }
});
}());
/**
* Generic inverse: coarse grid start, damped Newton with a numerical
* Jacobian and backtracking, then a derivative-free pattern search as a
* fallback. Projections need only supply forward(); a closed-form inverse
* is an optimisation, not a requirement.
*
* `tol` is the residual, in projected units, below which the point counts
* as being on the map. Callers derive it from a pixel tolerance and their
* own scale — a fixed number would mean different things on a small phone
* and a large screen, and different things again for projections whose
* extents range from 5.1 to 6.9.
*
* Returns {lon, lat, ok}. ok is false when the point lies outside the map.
*
* @param {Object} P A registered projection.
* @param {number} x Projected x.
* @param {number} y Projected y.
* @param {Array<number>} [seed] Starting guess as [lon, lat]; ignored
* when the projection has an exact inverse of its own.
* @param {number} [tol] Tolerance in projected units.
* @returns {{ok: boolean, lon: number, lat: number}} ok is false when
* the point lies outside the sheet.
* @alias JXG.Geography.invert
*/
G.invert = function (P, x, y, seed, tol) {
// A projection that knows its own inverse is asked directly; Newton is
// only for the ones given as a formula in one direction.
if (typeof P.inverse === 'function') {
var back = P.inverse(x, y);
var img = P.forward(back[0], back[1]);
return { lon: back[0], lat: back[1],
ok: Math.hypot(img[0] - x, img[1] - y) < (tol > 0 ? tol : 0.012) };
}
var LM = P.latMax - 1e-6, lon, lat, r, i, best, a, b, d, f, fa, fb,
j11, j12, j21, j22, det, dx, dy, dl, dp, step, moved, k, nl, np, nr, st;
function res(A, B) { var q = P.forward(A, B); return Math.hypot(q[0] - x, q[1] - y); }
if (seed) {
lon = clamp(seed[0], -180, 180); lat = clamp(seed[1], -LM, LM);
} else {
best = Infinity;
// The grid deliberately avoids the pole lines: y is stationary in
// latitude there and the Jacobian is singular.
for (a = -180; a <= 180; a += 10) {
for (b = -LM + 1; b <= LM - 1; b += 10) {
d = res(a, b);
if (d < best) { best = d; lon = a; lat = b; }
}
}
}
r = res(lon, lat);
var h = 1e-6;
for (i = 0; i < 80 && r > 1e-12; i++) {
f = P.forward(lon, lat);
fa = P.forward(clamp(lon + h, -180, 180), lat);
fb = P.forward(lon, clamp(lat + h, -LM, LM));
j11 = (fa[0] - f[0]) / h; j12 = (fb[0] - f[0]) / h;
j21 = (fa[1] - f[1]) / h; j22 = (fb[1] - f[1]) / h;
det = j11 * j22 - j12 * j21;
if (Math.abs(det) < 1e-14) { lat *= 0.995; r = res(lon, lat); continue; }
dx = f[0] - x; dy = f[1] - y;
dl = (-dx * j22 + dy * j12) / det;
dp = (dx * j21 - dy * j11) / det;
step = 1; moved = false;
for (k = 0; k < 28; k++) {
nl = clamp(lon + dl * step, -180, 180);
np = clamp(lat + dp * step, -LM, LM);
nr = res(nl, np);
if (nr < r) { lon = nl; lat = np; r = nr; moved = true; break; }
step *= 0.5;
}
if (!moved) { break; }
}
var lim = tol > 0 ? tol : 0.012;
if (r >= lim) {
st = 8;
while (st > 1e-8) {
moved = false;
var offs = [[st, 0], [-st, 0], [0, st], [0, -st],
[st, st], [st, -st], [-st, st], [-st, -st]];
for (k = 0; k < offs.length; k++) {
nl = clamp(lon + offs[k][0], -180, 180);
np = clamp(lat + offs[k][1], -LM, LM);
nr = res(nl, np);
if (nr < r) { lon = nl; lat = np; r = nr; moved = true; }
}
if (!moved) { st *= 0.5; }
}
}
return { lon: lon, lat: lat, ok: r < lim };
};
// ----------------------------------------------- graticule and boundaries
/**
* The outline of the sheet, in the coordinates the projection reads.
*
* For a cylindrical map that is the rectangle bounded by the meridians
* ±180 and the parallels ±latMax. For an azimuthal one it is the rim of the
* disc: a single parallel at φ′ = 90 − rhoMax, which is what the projection
* turns into a circle.
*
* @param {Object} P A registered projection.
* @param {number} n Samples per edge, or around the rim for a disc.
* @returns {Array<Array<number>>} The outline as [lon, lat] in degrees.
* @alias JXG.Geography.frameLL
*/
G.frameLL = function (P, n) {
var out = [], i;
if (P.clip === 'circle') {
var phiC = 90 - Math.min(P.rhoMax, 179.9);
for (i = 0; i <= n; i++) { out.push([-180 + 360 * i / n, phiC]); }
return out;
}
var L = P.latMax;
for (i = 0; i <= n; i++) { out.push([180, -L + 2 * L * i / n]); }
for (i = 0; i <= n; i++) { out.push([180 - 360 * i / n, L]); }
for (i = 0; i <= n; i++) { out.push([-180, L - 2 * L * i / n]); }
for (i = 0; i <= n; i++) { out.push([-180 + 360 * i / n, -L]); }
return out;
};
// Generated independently of the projection so that point counts stay equal
// across a switch between projections with different latMax — the morph
// depends on it. Clamping happens at projection time.
/**
* A meridian, pole to pole.
*
* @param {number} lon Longitude in degrees.
* @param {number} [n] Number of samples.
* @returns {Array<Array<number>>} Points as [lon, lat] in degrees.
* @alias JXG.Geography.meridianLL
*/
G.meridianLL = function (lon, n) {
var out = [], i;
for (i = 0; i <= n; i++) { out.push([lon, -90 + 180 * i / n]); }
return out;
};
/**
* A parallel, right around the sphere.
*
* The ring is open and runs eastward, so it circles the north pole
* counter-clockwise seen from outside — which is what `orientationOf`
* reports for it.
*
* @param {number} lat Latitude in degrees.
* @param {number} [n] Number of samples.
* @returns {Array<Array<number>>} Points as [lon, lat] in degrees.
* @alias JXG.Geography.features
* @alias JXG.Geography.parallelLL
*/
G.parallelLL = function (lat, n) {
var out = [], i;
for (i = 0; i <= n; i++) { out.push([-180 + 360 * i / n, lat]); }
return out;
};
// ------------------------------------------------------------- features
G.features = {
equator: function () { return [G.parallelLL(0, 240)]; },
// both halves together, and each on its own: the northern tropic is the
// Tropic of Cancer, the southern the Tropic of Capricorn, and they are
// rarely wanted at the same time
tropics: function () { return [G.parallelLL(OBLIQUITY, 240), G.parallelLL(-OBLIQUITY, 240)]; },
tropicnorth: function () { return [G.parallelLL(OBLIQUITY, 240)]; },
tropicsouth: function () { return [G.parallelLL(-OBLIQUITY, 240)]; },
polarcircles: function () {
return [G.parallelLL(90 - OBLIQUITY, 240), G.parallelLL(-(90 - OBLIQUITY), 240)];
},
arcticcircle: function () { return [G.parallelLL(90 - OBLIQUITY, 240)]; },
antarcticcircle: function () { return [G.parallelLL(-(90 - OBLIQUITY), 240)]; },
primemeridian: function () { return [G.meridianLL(0, 120)]; },
dateline: function () { return [G.meridianLL(180, 120), G.meridianLL(-180, 120)]; },
timezones: function () {
var o = [], l;
for (l = -180; l < 180; l += 15) { if (l % 30) { o.push(G.meridianLL(l, 120)); } }
return o;
},
terminator: function (date) {
var c = terminator(date, 240);
return [c.concat([c[0]])];
},
subsolar: function (date) {
var s = subsolar(date), c = smallCircle(s[0], s[1], 2.2, 20);
return [c.concat([c[0]])];
}
};
/**
* Nominal time zones: the 15-degree bands, and one of them singled out.
*
* These are the meridian-based zones, not the legal ones — a country's real
* zone follows its border and can be offset by 30 or 45 minutes. What can be
* had from a formula is the nominal band; the legal boundaries need their
* own dataset, and the two disagree by up to several hundred kilometres.
*
* @param {number} offsetHours Offset from UTC, in hours.
* @returns {Array<Array<number>>} The band as a closed ring of
* [lon, lat]. The nominal 15-degree slice, not the legal boundary.
* @alias JXG.Geography.timeZoneBand
*/
G.timeZoneBand = function (offsetHours) {
var centre = offsetHours * 15;
var west = centre - 7.5, east = centre + 7.5, out = [], i;
for (i = 0; i <= 60; i++) { out.push([west, -90 + 180 * i / 60]); }
for (i = 0; i <= 60; i++) { out.push([west + 15 * i / 60, 90]); }
for (i = 0; i <= 60; i++) { out.push([east, 90 - 180 * i / 60]); }
for (i = 0; i <= 60; i++) { out.push([east - 15 * i / 60, -90]); }
return out.map(function (p) {
var lon = p[0];
while (lon > 180) { lon -= 360; }
while (lon < -180) { lon += 360; }
return [lon, p[1]];
});
};
/**
* Local time in a nominal zone, with an optional summer-time rule.
*
* `rule` is 'eu', 'us' or 'none'. Both rules are stated in local time and
* differ by more than their dates: the European change happens at the same
* instant everywhere, the American one at 2 a.m. in each zone separately.
*
* @param {number} offsetHours Standard offset from UTC, in hours.
* @param {Date} [when] The moment; defaults to now.
* @param {string} [rule] Summer-time rule: "eu", "us" or "none".
* @returns {Date} Local time, summer time included where the rule says.
* @alias JXG.Geography.localTime
*/
G.localTime = function (offsetHours, when, rule) {
var d = when || new Date();
var dst = G.summerTime(d, rule, offsetHours);
var t = new Date(d.getTime() + (offsetHours + (dst ? 1 : 0)) * 3600000);
return { hours: t.getUTCHours(), minutes: t.getUTCMinutes(),
summer: dst, offset: offsetHours + (dst ? 1 : 0) };
};
/**
* Whether summer time is in force. Northern-hemisphere rules only.
*
* @param {Date} when The moment to test.
* @param {string} rule "eu", "us" or "none".
* @param {number} offsetHours Standard offset from UTC, in hours.
* @returns {boolean} True when summer time is in force.
* @alias JXG.Geography.summerTime
*/
G.summerTime = function (when, rule, offsetHours) {
if (!rule || rule === 'none') { return false; }
var y = when.getUTCFullYear(), start, end;
function lastSunday(month) { // month is 0-based
var d = new Date(Date.UTC(y, month + 1, 0));
d.setUTCDate(d.getUTCDate() - d.getUTCDay());
return d;
}
function nthSunday(month, n) {
var d = new Date(Date.UTC(y, month, 1));
d.setUTCDate(1 + ((7 - d.getUTCDay()) % 7) + (n - 1) * 7);
return d;
}
if (rule === 'eu') {
// 01:00 UTC on the last Sunday in March and in October, everywhere at once
start = lastSunday(2); start.setUTCHours(1);
end = lastSunday(9); end.setUTCHours(1);
} else {
// 02:00 local on the second Sunday in March and the first in November
start = nthSunday(2, 2); start.setUTCHours(2 - (offsetHours || 0));
end = nthSunday(10, 1); end.setUTCHours(2 - (offsetHours || 0));
}
return when >= start && when < end;
};
/**
* Tissot indicatrices: circles of equal angular radius on a lattice.
*
* @param {Object} [opts] `radius` in degrees, `lonStep` and `latStep`
* in degrees, `n` samples per circle.
* @returns {Array<Array<Array<number>>>} One ring per indicatrix.
* @alias JXG.Geography.tissot
*/
G.tissot = function (opts) {
opts = opts || {};
var step = opts.step || 30, r = opts.radius || 6, out = [], lat, lon, c;
var latMax = opts.latMax === undefined ? 60 : opts.latMax;
for (lat = -latMax; lat <= latMax; lat += step) {
// starts at -180 so the antimeridian carries a circle too; the old
// offset of 1.5 steps left a visible gap there
for (lon = -180; lon < 180; lon += step * 1.5) {
c = smallCircle(lon, lat, r, 40);
out.push(c.concat([c[0]]));
}
}
return out;
};
// ---------------------------------------------------------- dataset access
G.datasets = {};
/**
* Add a dataset under a name, so an element can ask for it by that name.
*
* @param {string} name The key it is looked up by.
* @param {Object} ds The dataset. It is not validated here; run
* `G.validate` when building one.
* @returns {Object} The dataset as registered.
* @alias JXG.Geography.registerDataset
*/
G.registerDataset = function (name, ds) { ds.id = ds.id || name; G.datasets[name] = ds; return ds; };
/**
* Bounding box of a group of rings, as [[lonW, latS], [lonE, latN]].
*
* Longitudes are unwrapped against the first point, so a shape straddling
* the antimeridian yields a west edge east of its east edge — Fiji comes out
* as 175 to 185 rather than as the whole world. Callers that need a plain
* pair can fold the east edge back themselves; a viewport must not.
*/
/* ---------------------------------------------------------------------
* Feature layers
*
* A dataset holds groups of things that differ only in what they are called:
* countries, states, lakes, rivers, cities. Rather than a function per kind,
* they are all *layers* — a list of features with the same three questions
* asked of them: where is it, what is it called, and how important is it.
*
* feature = { id, name, parent?, rank?, bbox, polys | line | point }
*
* `parent` carries the hierarchy: a state's parent is its country. `rank`
* carries importance, so a view can draw only what its scale justifies —
* the one thing that keeps a detailed dataset usable.
* ------------------------------------------------------------------- */
/**
* Lat-band buckets, so a lookup need not walk every feature.
*
* @param {Array<Object>} features Features carrying a `bbox`.
* @param {number} [size=10] Height of a band, in degrees.
* @returns {{size: number, buckets: Object}} Feature ids per band.
* @alias JXG.Geography.bandIndex
*/
G.bandIndex = function (features, size) {
var n = size || 10, buckets = {}, i, f, lo, hi, b;
for (i = 0; i < features.length; i++) {
f = features[i];
if (!f.bbox) { continue; }
lo = Math.floor((f.bbox[1] + 90) / n);
hi = Math.floor((f.bbox[3] + 90) / n);
for (b = lo; b <= hi; b++) {
(buckets[String(b)] = buckets[String(b)] || []).push(f.id);
}
}
return { size: n, buckets: buckets };
};
/**
* The existing groups, presented as layers without copying their geometry.
*
* `countries` and the rest keep the shape they have had since the first
* build. This is an adapter, not a migration: both readings work on the same
* arrays, so nothing has to be converted and nothing can drift apart.
*
* @param {Object} ds A dataset, as registered with `registerDataset`.
* @param {string} name A group in the dataset, such as "countries".
* @returns {?Object} A layer, or null with a warning if there is no such
* group. The features are the dataset's own arrays, not copies.
* @alias JXG.Geography.asLayer
*/
G.asLayer = function (ds, name) {
var src = ds[name];
if (!src) { warn('no group "' + name + '" in this dataset'); return null; }
var kind = name === 'capitals' ? 'points' : 'polys';
var features = src.map ? src : [];
return { name: name, kind: kind, features: features,
index: (name === 'countries' && ds.index) ? ds.index.bands
: G.bandIndex(features) };
};
/**
* Which feature of this layer contains the point, or null.
*
* @param {Object} layer A feature layer.
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @returns {?Object} The feature containing the position, or null. A
* position in a hole counts as outside.
* @alias JXG.Geography.featureAt
*/
G.featureAt = function (layer, lon, lat) {
var ids = null, i, j, k, f, poly, byId;
if (layer.index && layer.index.buckets) {
ids = layer.index.buckets[String(Math.floor((lat + 90) / layer.index.size))] || [];
byId = {};
for (i = 0; i < layer.features.length; i++) {
byId[layer.features[i].id] = layer.features[i];
}
}
var list = ids ? ids.map(function (id) { return byId[id]; }) : layer.features;
for (i = 0; i < list.length; i++) {
f = list[i];
if (!f || !f.polys) { continue; }
if (f.bbox && (lon < f.bbox[0] || lon > f.bbox[2] ||
lat < f.bbox[1] || lat > f.bbox[3])) { continue; }
for (j = 0; j < f.polys.length; j++) {
poly = f.polys[j];
if (!inRing(poly[0], lon, lat)) { continue; }
for (k = 1; k < poly.length; k++) {
if (inRing(poly[k], lon, lat)) { break; } // in a hole
}
if (k === poly.length) { return f; }
}
}
return null;
};
/**
* The features of this layer worth drawing in `box`, at or above `minRank`.
*
* Culling by bounding box is what makes a detailed dataset affordable: at
* country scale nothing outside the view is clipped, projected or drawn.
*
* @param {Object} layer A feature layer.
* @param {?Array<Array<number>>} box [[lonW, latS], [lonE, latN]], or
* null for no culling. May straddle the antimeridian.
* @param {number} [minRank] Keep features at or below this rank. One
* without a rank is always kept.
* @returns {Array<Object>} The features that survive.
* @alias JXG.Geography.featuresIn
*/
G.featuresIn = function (layer, box, minRank) {
var w, e, s2, n2, out = [], i, f;
if (box) {
w = box[0][0]; s2 = box[0][1]; e = box[1][0]; n2 = box[1][1];
}
for (i = 0; i < layer.features.length; i++) {
f = layer.features[i];
if (minRank !== undefined && minRank !== null &&
f.rank !== undefined && f.rank > minRank) { continue; }
if (box && f.bbox) {
if (f.bbox[1] > n2 || f.bbox[3] < s2) { continue; }
// longitudes are compared in the box's own frame, so a viewport across
// the antimeridian keeps working
var fw = f.bbox[0], fe = f.bbox[2];
while (fw - w > 180) { fw -= 360; fe -= 360; }
while (w - fw > 180) { fw += 360; fe += 360; }
if (fw > e || fe < w) { continue; }
}
out.push(f);
}
return out;
};
/**
* Every ring of a layer's features, ready for addGeoLayer.
*
* @param {Array<Object>} features Features with a `polys` member.
* @returns {Array<Array<Array<number>>>} Every ring, outer and hole
* alike, in the order they appear.
* @alias JXG.Geography.ringsOf
*/
G.ringsOf = function (features) {
var out = [], i, j, k;
for (i = 0; i < features.length; i++) {
if (!features[i].polys) { continue; }
for (j = 0; j < features[i].polys.length; j++) {
for (k = 0; k < features[i].polys[j].length; k++) {
out.push(features[i].polys[j][k]);
}
}
}
return out;
};
/**
* Bounding box of a group of rings, as [[lonW, latS], [lonE, latN]].
*
* Longitudes are unwrapped against the first point, so a shape straddling
* the antimeridian yields a west edge east of its east edge — Fiji comes out
* 2.9 degrees wide rather than as the whole world, and Russia keeps running
* past 180. A viewport needs that; a label would need the fold back.
*
* @param {Array<Array<Array<number>>>} rings Rings of [lon, lat] in degrees.
* @returns {?Array<Array<number>>} [[lonW, latS], [lonE, latN]] in degrees,
* or null for an empty group. The east edge may exceed 180.
* @alias JXG.Geography.bboxOf
*/
G.bboxOf = function (rings) {
var w = Infinity, e = -Infinity, s = Infinity, n = -Infinity, ref = null;
rings.forEach(function (r) {
r.forEach(function (p) {
var lon = p[0];
if (ref === null) { ref = lon; }
while (lon - ref > 180) { lon -= 360; }
while (ref - lon > 180) { lon += 360; }
if (lon < w) { w = lon; }
if (lon > e) { e = lon; }
if (p[1] < s) { s = p[1]; }
if (p[1] > n) { n = p[1]; }
});
});
if (ref === null) { return null; }
return [[w, s], [e, n]];
};
/**
* The outline of a lon/lat box, sampled along its edges.
*
* Corners alone will not do: every projection bends a straight edge, so a
* box fitted from its four corners comes out too small.
*
* @param {Array<Array<number>>} box [[lonW, latS], [lonE, latN]] in degrees.
* @param {number} [n=60] Samples per edge.
* @returns {Array<Array<number>>} The outline as [lon, lat], closed.
* @alias JXG.Geography.boxLL
*/
G.boxLL = function (box, n) {
var w = box[0][0], s = box[0][1], e = box[1][0], t = box[1][1];
var out = [], i, k = n || 60;
for (i = 0; i <= k; i++) { out.push([w + (e - w) * i / k, s]); }
for (i = 0; i <= k; i++) { out.push([e, s + (t - s) * i / k]); }
for (i = 0; i <= k; i++) { out.push([e - (e - w) * i / k, t]); }
for (i = 0; i <= k; i++) { out.push([w, t - (t - s) * i / k]); }
return out;
};
/**
* A country is either polygons or a point below the resolution.
*
* @param {Object} country A country from the dataset.
* @param {number} smallRadius Radius of the stand-in circle, in degrees.
* @param {boolean} showSmall Whether to draw stand-ins at all. A circle
* is a symbol at the wrong size, not the country's shape, so a map
* meant to be measured from should leave them out.
* @returns {Array<Array<Array<number>>>} Rings ready to draw.
* @alias JXG.Geography.countryRings
*/
G.countryRings = function (country, smallRadius, showSmall) {
var out = [], i, j;
if (country.polys) {
for (i = 0; i < country.polys.length; i++) {
for (j = 0; j < country.polys[i].length; j++) { out.push(country.polys[i][j]); }
}
return out;
}
// A country too small to have an outline at this resolution is stood in
// for by a circle. The circle is a symbol, not a shape: it is the wrong
// size and the wrong outline, so a map that wants only true geometry can
// leave these out.
if (showSmall === false) { return []; }
return [smallCircle(country.point[0], country.point[1], smallRadius || 1.15, 14)];
};
/**
* Countries represented by a symbol rather than an outline.
*
* @param {Object} ds A dataset, as registered with `registerDataset`.
* @returns {Array<string>} Ids of the countries with no outline at this
* resolution.
* @alias JXG.Geography.smallCountries
*/
G.smallCountries = function (ds) {
return (ds.countries || []).filter(function (c) { return !c.polys; });
};
/**
* Is the position inside this ring?
*
* Ray casting in longitude and latitude, so the ring must not straddle the
* antimeridian — the dataset guarantees that, and `validate` checks it.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @returns {boolean} True when the position lies inside.
* @alias JXG.Geography.inRing
*/
function inRing(ring, lon, lat) {
var inside = false, i, j, yi, yj;
for (i = 0, j = ring.length - 1; i < ring.length; j = i++) {
yi = ring[i][1]; yj = ring[j][1];
if ((yi > lat) !== (yj > lat) &&
lon < (ring[j][0] - ring[i][0]) * (lat - yi) / (yj - yi) + ring[i][0]) {
inside = !inside;
}
}
return inside;
}
G.inRing = inRing;
/**
* Point in polygon over the latitude-band index instead of every country.
*
* @param {Object} ds A dataset, as registered with `registerDataset`.
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @returns {?string} The country id, or null for open water.
*
* @alias JXG.Geography.countryAt
*/
G.countryAt = function (ds, lon, lat) {
var band = ds.index && ds.index.bands, ids, i, k, c, byId = {}, poly, j;
for (i = 0; i < ds.countries.length; i++) { byId[ds.countries[i].id] = ds.countries[i]; }
if (band) {
ids = band.buckets[String(Math.floor((lat + 90) / band.size))] || [];
} else {
ids = ds.countries.map(function (c2) { return c2.id; });
}
for (i = 0; i < ids.length; i++) {
c = byId[ids[i]];
if (!c || !c.polys) { continue; }
if (lon < c.bbox[0] || lon > c.bbox[2] || lat < c.bbox[1] || lat > c.bbox[3]) { continue; }
for (j = 0; j < c.polys.length; j++) {
poly = c.polys[j];
if (inRing(poly[0], lon, lat)) {
for (k = 1; k < poly.length; k++) {
if (inRing(poly[k], lon, lat)) { break; }
}
if (k === poly.length) { return c.id; }
}
}
}
return null;
};
/**
* The nearest point of a dataset group to a position.
*
* @param {Object} ds A dataset.
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @param {string} [layer="capitals"] * @param {Object} ds A dataset.
* @param {number} lon Longitude in degrees.
* @param {number} lat Latitude in degrees.
* @param {string} [layer="capitals"] Which group of points to search.
* @returns {?Object} The nearest entry, with `km` added, or null when
* the group is empty.
* @alias JXG.Geography.nearest
*/
G.nearest = function (ds, lon, lat, layer) {
var list = ds[layer || 'capitals'], best = null, bd = Infinity, i, d;
for (i = 0; i < list.length; i++) {
d = angle([lon, lat], [list[i].lon, list[i].lat]);
if (d < bd) { bd = d; best = list[i]; }
}
return best ? { item: best, distance: bd * RAD * R_EARTH } : null;
};
// ------------------------------------------------------------- clip cache
/**
* Rotation and antimeridian clipping depend on the map centre only, never
* on the projection. Results are therefore cached per (dataset, centre) and
* shared between hosts. The same invariant lets the projection morph
* interpolate: with the centre unchanged, point counts match.
*
* @param {Object} ds The dataset being clipped.
* @param {Array<number>} centre Map centre as [lon, lat] in degrees.
* @param {boolean} oblique Whether the rotation tilts the axis.
* @param {number} [maxStep] Densification step, in degrees.
* @param {number} [rhoMax] Disc radius in degrees; 0 for a sheet.
* Two states with different rhoMax cannot share a memo, because a
* disc and a sheet cut different rings in different places.
*/
function ClipState(ds, centre, oblique, maxStep, rhoMax) {
this.dataset = ds;
this.centre = [centre[0], centre[1]];
this.oblique = !!oblique;
this.maxStep = maxStep || ds.maxSegment || 2;
this.lat0 = this.oblique ? centre[1] : 0;
// An azimuthal map cuts against a circle about its centre, so it works in
// pole-rotated coordinates. Everything else cuts against a meridian and
// puts the centre on the equator. The two cannot share a cache entry,
// which is why rhoMax is part of the key.
this.rhoMax = rhoMax || 0;
this.circle = this.rhoMax > 0;
if (this.circle) {
this.rot = poleRotator(centre[0], centre[1]);
this.unrot = poleUnrotator(centre[0], centre[1]);
} else {
this.rot = rotator(centre[0], this.lat0);
this.unrot = unrotator(centre[0], this.lat0);
}
this.refs = 0;
this._memo = {};
}
ClipState.prototype.rings = function (key, source) {
if (this._memo[key]) { return this._memo[key]; }
var parts = [], i, self = this;
for (i = 0; i < source.length; i++) {
parts = parts.concat(self.circle
? clipRingCircle(source[i], self.rot, self.rhoMax, self.maxStep)
: clipRing(source[i], self.rot, self.maxStep));
}
this._memo[key] = join(parts, true);
return this._memo[key];
};
ClipState.prototype.lines = function (key, source) {
if (this._memo[key]) { return this._memo[key]; }
var parts = [], i;
for (i = 0; i < source.length; i++) {
parts = parts.concat(this.circle
? clipLineCircle(source[i], this.rot, this.rhoMax)
: clipLine(source[i], this.rot));
}
this._memo[key] = join(parts, false);
return this._memo[key];
};
ClipState.prototype.forget = function (key) { delete this._memo[key]; };
var cache = {};
G.clipCache = {
key: function (ds, centre, oblique, rhoMax) {
return (ds.id || 'anon') + '|' + centre[0].toFixed(6) + '|' +
(oblique ? centre[1].toFixed(6) : '0') + '|' + (oblique ? 1 : 0) +
'|' + (rhoMax || 0);
},
acquire: function (ds, centre, oblique, rhoMax) {
var k = this.key(ds, centre, oblique, rhoMax);
if (!cache[k]) {
cache[k] = new ClipState(ds, centre, oblique, null, rhoMax);
}
cache[k].refs++;
return cache[k];
},
release: function (state) {
var k = this.key(state.dataset, state.centre, state.oblique);
if (!cache[k]) { return; }
cache[k].refs--;
if (cache[k].refs <= 0) { delete cache[k]; }
},
size: function () { return Object.keys(cache).length; },
clear: function () { cache = {}; }
};
// ------------------------------------------------------------ masked layer
/**
* `occlude` is the distance of the points from the centre in Earth radii.
* On the surface it is 1 and the half-space test is right. Further out a
* point is only hidden when it is behind the plane *and* inside the Earth's
* silhouette — otherwise it is seen past the limb, which is exactly what
* makes a satellite visible above the far side.
*
* @param {Array<Array<number>>} pts Points as [lon, lat] in degrees.
* @param {Array<number>} cam Camera direction as a unit vector.
* @param {number} side +1 for the near hemisphere, -1 for the far one.
* @param {number} radius Sphere radius in world units.
* @param {Array<number>} X Output array for x, NaN between pieces.
* @param {Array<number>} Y Output array for y.
* @param {Array<number>} Z Output array for z.
* @param {number} [occlude] Radius factor a point must clear to be seen
* past the limb; 1 is the surface itself.
* @returns {void} The three arrays are filled in place.
*
* @alias JXG.Geography.maskLine
*/
G.maskLine = function (pts, cam, side, radius, X, Y, Z, occlude) {
var n = pts.length, base = X.length, dPrev = 0, vPrev = false, i, q, d, v, t, a, w, L, k;
var kk = occlude || 1;
// Where the visibility boundary sits. On the surface it is the terminator,
// d = 0. Further out the horizon is the silhouette cylinder, and the
// boundary moves to d = -sqrt(1 - 1/k^2). Inserting the crossing at d = 0
// regardless would place it in the wrong spot and leave a stray segment.
var dLim = kk > 1 ? -Math.sqrt(1 - 1 / (kk * kk)) : 0;
for (i = 0; i < n; i++) {
q = pts[i];
d = (q[0] * cam[0] + q[1] * cam[1] + q[2] * cam[2]) * side;
v = d > dLim;
if (v) {
X.push(q[0] * radius);
Y.push(q[1] * radius);
Z.push(q[2] * radius);
} else {
X.push(NaN); Y.push(NaN); Z.push(NaN);
}
if (i > 0 && v !== vPrev) {
a = pts[i - 1];
t = (dPrev - dLim) / (dPrev - d);
// The linear estimate is exact only for dLim = 0: normalising the
// interpolated vector changes its length, and with it the dot
// product, unless the target is zero. A few bisection steps on the
// normalised point put the crossing back on the silhouette.
if (dLim !== 0) {
var lo = 0, hi = 1, m, u, uL, ud, s2;
for (s2 = 0; s2 < 20; s2++) {
m = (lo + hi) / 2;
u = [a[0] + m * (q[0] - a[0]), a[1] + m * (q[1] - a[1]), a[2] + m * (q[2] - a[2])];
uL = Math.hypot(u[0], u[1], u[2]);
ud = (u[0] * cam[0] + u[1] * cam[1] + u[2] * cam[2]) / uL * side;
if ((ud > dLim) === (dPrev > dLim)) { lo = m; } else { hi = m; }
}
t = (lo + hi) / 2;
}
w = [a[0] + t * (q[0] - a[0]), a[1] + t * (q[1] - a[1]), a[2] + t * (q[2] - a[2])];
L = Math.hypot(w[0], w[1], w[2]);
k = base + (v ? i - 1 : i);
X[k] = w[0] / L * radius; Y[k] = w[1] / L * radius; Z[k] = w[2] / L * radius;
}
dPrev = d; vPrev = v;
}
};
/**
* Visibility mask for a polyline drawn on the globe: the parts behind the
* sphere are cut away and the piece boundaries marked with NaN.
*
* @param {Object} el A JSXGraph curve3d to drive.
* @param {Object} opts How to read the scene:
* `segments()` yields the polylines as [lon, lat];
* `camera()` the view direction as a unit vector;
* `project(lon, lat, r)` places a point in world coordinates;
* `radius()` the sphere radius, default 1;
* `side` +1 for the near hemisphere, -1 for the far one, default +1;
* `occlude` the radius factor a point must clear to be seen past the
* limb, default 1;
* `active` false leaves the element empty.
* @returns {Object} The element, for chaining.
* @alias JXG.Geography.maskedLayer
*/
G.maskedLayer = function (el, opts) {
var segments = opts.segments, camera = opts.camera, project = opts.project,
radius = opts.radius || function () { return 1; },
side = opts.side === undefined ? 1 : opts.side;
var occlude = opts.occlude || 1;
el.geoActive = opts.active !== false;
el.updateDataArray2D = function () {
var X = [], Y = [], Z = [], segs, cam, r, i, out, j, P;
if (!this.geoActive) { return { X: [], Y: [] }; }
r = radius();
if (!(r > 0)) { return { X: [], Y: [] }; } // radius 0 draws nothing
segs = segments(); cam = camera();
for (i = 0; i < segs.length; i++) {
// An empty piece is normal: a trail has none until it is first
// filled, and a layer may be created before its source exists.
if (!segs[i] || segs[i].length === 0) { continue; }
if (segs[i].length > 1) { G.maskLine(segs[i], cam, side, r, X, Y, Z, occlude); } else {
var vis = G.dot(segs[i][0], cam) * side > 0;
X.push(vis ? segs[i][0][0] * r : NaN);
Y.push(vis ? segs[i][0][1] * r : NaN);
Z.push(vis ? segs[i][0][2] * r : NaN);
}
X.push(NaN); Y.push(NaN); Z.push(NaN);
}
out = { X: [], Y: [] };
for (j = 0; j < X.length; j++) {
if (X[j] !== X[j]) { out.X.push(NaN); out.Y.push(NaN); continue; }
P = project([X[j], Y[j], Z[j]]);
out.X.push(P[1]); out.Y.push(P[2]);
}
return out;
};
return el;
};
/**
* Mean direction of a ring, and whether it clears the Earth's silhouette.
*
* @param {Array<Array<number>>} pts The ring, as world vectors.
* @param {Array<number>} cam Camera direction as a unit vector.
* @param {number} side +1 for the near hemisphere, -1 for the far one.
* @param {number} k Radius factor: how far above the surface it sits.
* @returns {boolean} Whether the whole ring is seen. All or nothing,
* because clipping a small marker against the limb would close it
* over the wrong arc and make it jump as it crosses.
*/
function visibleAbove(pts, cam, side, k) {
var m = [0, 0, 0], i, L, d;
for (i = 0; i < pts.length; i++) { m[0] += pts[i][0]; m[1] += pts[i][1]; m[2] += pts[i][2]; }
L = Math.hypot(m[0], m[1], m[2]);
if (L < 1e-12) { return false; }
d = (m[0] * cam[0] + m[1] * cam[1] + m[2] * cam[2]) / L * side;
return d > 0 || k * Math.sqrt(Math.max(0, 1 - d * d)) > 1;
}
/**
* Filled counterpart to `maskedLayer`: rings clipped against the terminator
* and rejoined along the limb, so a coastline running off the edge closes
* against the silhouette instead of across the middle.
*
* @param {Object} el A JSXGraph curve3d to drive.
* @param {Object} opts As `maskedLayer`, plus `orients()` giving each
* ring's winding — stated rather than derived, because a great circle
* spans every longitude and no measure can read its orientation.
* @returns {Object} The element, for chaining.
* @alias JXG.Geography.maskedRings
*/
G.maskedRings = function (el, opts) {
// the camera is read through frame().cam, so opts.camera is not needed here
var rings = opts.rings, frame = opts.frame,
project = opts.project, side = opts.side === undefined ? 1 : opts.side,
radius = opts.radius || function () { return 1; },
orients = opts.orients || null,
occlude = opts.occlude || 1,
arcStep = opts.arcStep || 3 * RAD;
el.geoActive = opts.active !== false;
el.updateDataArray2D = function () {
var out = { X: [], Y: [] }, src, f, r, i, k, parts, q, P;
if (!this.geoActive) { return out; }
r = radius();
if (!(r > 0)) { return out; }
src = rings(); f = frame();
for (i = 0; i < src.length; i++) {
if (!src[i] || src[i].length < 3) { continue; }
if (occlude > 1) {
// A ring above the surface is not cut by the terminator: it is seen
// past the limb wherever it clears the silhouette. These rings are
// small markers, so the honest answer is all or nothing — clipping
// them against the sphere's silhouette would close them over the
// wrong arc and make them jump as they cross it.
parts = visibleAbove(src[i], f.cam, side, occlude) ? [src[i].concat([src[i][0]])] : [];
} else {
parts = clipSphereRing(src[i], f.r1, f.r2, f.cam, side, arcStep,
orients ? orients()[i] : 1);
}
for (k = 0; k < parts.length; k++) {
if (out.X.length) { out.X.push(NaN); out.Y.push(NaN); }
for (q = 0; q < parts[k].length; q++) {
P = project([parts[k][q][0] * r, parts[k][q][1] * r, parts[k][q][2] * r]);
out.X.push(P[1]); out.Y.push(P[2]);
}
}
}
return out;
};
return el;
};
// -------------------------------------------------------------- validation
/**
* The six guarantees the renderer relies on. See SPEC.md §7.
*
* @param {Object} ds A dataset, as registered with `registerDataset`.
* @returns {Array<string>} One message per broken invariant; empty when
* the dataset is sound.
* @alias JXG.Geography.validate
*/
G.validate = function (ds) {
var errs = [], ms = ds.maxSegment || 2;
function checkRing(ring, ccw, where) {
var i, a, b;
if (ring.length < 3) { errs.push(where + ': ring has only ' + ring.length + ' points'); return; }
if (ring[0][0] === ring[ring.length - 1][0] && ring[0][1] === ring[ring.length - 1][1]) {
errs.push(where + ': ring stored closed');
}
if ((shoelace(ring) > 0) !== ccw) { errs.push(where + ': wrong winding order'); }
for (i = 0; i < ring.length; i++) {
a = ring[i]; b = ring[(i + 1) % ring.length];
if (Math.abs(a[0] - b[0]) > 180) { errs.push(where + ': segment spans the antimeridian'); break; }
if (Math.max(Math.abs(a[0] - b[0]), Math.abs(a[1] - b[1])) > ms * 1.5) {
errs.push(where + ': segment longer than maxSegment'); break;
}
if (Math.abs(a[0]) > 180.001 || Math.abs(a[1]) > 90.001) {
errs.push(where + ': coordinate out of range'); break;
}
}
checkPoleClosure(ring, where);
}
/**
* A ring that spans every longitude either circles a pole, in which case
* it must reach it, or it does not, in which case it cannot span every
* longitude. Antarctica in `countries` did the impossible: 360 degrees
* wide and stopping at −85.61, so a map filling that layer showed a
* horizontal cut with whatever lay beneath it shining through. `land` and
* `coast` carry 183 points on the pole line; `countries` carried none.
*
* @param {Array<Array<number>>} ring Open ring of [lon, lat] in degrees.
* @param {string} where How to name it in a complaint.
* @returns {void} Anything wrong is appended to the error list.
*/
function checkPoleClosure(ring, where) {
var lo = 180, hi = -180, s2 = 90, n2 = -90, i;
for (i = 0; i < ring.length; i++) {
if (ring[i][0] < lo) { lo = ring[i][0]; }
if (ring[i][0] > hi) { hi = ring[i][0]; }
if (ring[i][1] < s2) { s2 = ring[i][1]; }
if (ring[i][1] > n2) { n2 = ring[i][1]; }
}
if (hi - lo > 350 && s2 > -89.9 && n2 < 89.9) {
errs.push(where + ': ring spans every longitude but reaches neither ' +
'pole (' + s2.toFixed(2) + ' to ' + n2.toFixed(2) + ')');
}
}
function checkPolys(polys, where) {
var i, j;
for (i = 0; i < polys.length; i++) {
for (j = 0; j < polys[i].length; j++) {
checkRing(polys[i][j], j === 0, where + '[' + i + '][' + j + ']');
}
}
}
checkPolys(ds.land.polys, 'land');
ds.continents.forEach(function (c) { checkPolys(c.polys, c.id); });
ds.countries.forEach(function (c) { if (c.polys) { checkPolys(c.polys, c.id); } });
var ids = ds.countries.map(function (c) { return c.id; });
if (ids.length !== Object.keys(ids.reduce(function (m, i) { m[i] = 1; return m; }, {})).length) {
errs.push('duplicate country ids');
}
var caps = {};
ds.capitals.forEach(function (c) { caps[c.id] = 1; });
ds.countries.forEach(function (c) {
if (c.capital && !caps[c.capital]) { errs.push(c.id + ' refers to an unknown capital'); }
});
var csum = ds.continents.reduce(function (s, c) { return s + c.area; }, 0);
if (Math.abs(csum - ds.land.area) / ds.land.area > 0.02) {
errs.push('continent areas deviate ' +
(100 * (csum - ds.land.area) / ds.land.area).toFixed(2) + ' % from land');
}
return errs;
};
if (typeof module === 'object' && module.exports) { module.exports = G; }
}(typeof globalThis !== 'undefined' ? globalThis : this));