worldmap-viewer

git clone https://git.godosa.eu/worldmap-viewer

master

raw · 4180 bytes

// Sphere maths for World Maps. Scene convention (y up): x = cos φ cos λ, y = sin φ, z = −cos φ sin λ.
import { PLANET } from './planet.js';
export const R_KM = PLANET.radius_km;
const D = Math.PI / 180;

export function toVec(lat, lon) {
  const la = lat * D, lo = lon * D, c = Math.cos(la);
  return [c * Math.cos(lo), Math.sin(la), -c * Math.sin(lo)];
}

export function toLatLon(v) {
  const n = Math.hypot(v[0], v[1], v[2]) || 1;
  return { lat: Math.asin(Math.max(-1, Math.min(1, v[1] / n))) / D, lon: Math.atan2(-v[2], v[0]) / D };
}

const dot = (a, b) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
const cross = (a, b) => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];

export function angleRad(p, q) {
  const a = toVec(p.lat, p.lon), b = toVec(q.lat, q.lon);
  return Math.atan2(Math.hypot(...cross(a, b)), dot(a, b));
}

export const distanceKm = (p, q) => R_KM * angleRad(p, q);

export function interpolate(p, q, stepDeg = 1) {
  const a = toVec(p.lat, p.lon), b = toVec(q.lat, q.lon), w = angleRad(p, q);
  const n = Math.max(1, Math.ceil(w / (stepDeg * D) - 1e-9));   // tolerate rounding: 3° at 1° steps = 3 pieces
  const out = [];
  for (let i = 0; i <= n; i++) {
    if (w < 1e-12) { out.push({ lat: p.lat, lon: p.lon }); continue; }
    const t = i / n, s1 = Math.sin((1 - t) * w) / Math.sin(w), s2 = Math.sin(t * w) / Math.sin(w);
    out.push(toLatLon([a[0] * s1 + b[0] * s2, a[1] * s1 + b[1] * s2, a[2] * s1 + b[2] * s2]));
  }
  return out;
}

export function polygonAreaKm2(pts) {
  if (pts.length < 3) return 0;
  const v = pts.map(p => toVec(p.lat, p.lon));
  let e = 0;
  for (let i = 1; i + 1 < v.length; i++) {       // signed spherical excess of a triangle fan (Van Oosterom–Strackee)
    const a = v[0], b = v[i], c = v[i + 1];
    e += 2 * Math.atan2(dot(a, cross(b, c)), 1 + dot(a, b) + dot(b, c) + dot(c, a));
  }
  return Math.abs(e) * R_KM * R_KM;
}

// Equal Earth (Šavrič et al. 2018): same constants as map/mapgen/projections.py
const A1 = 1.340264, A2 = -0.081106, A3 = 0.000893, A4 = 0.003796, M = Math.sqrt(3) / 2, P3 = Math.PI / 3;
export const EE_XMAX = (2 * Math.sqrt(3) * Math.PI) / (3 * A1);
export const EE_YMAX = A1 * P3 + A2 * P3 ** 3 + A3 * P3 ** 7 + A4 * P3 ** 9;
export const PC_XMAX = Math.PI;
export const PC_YMAX = Math.PI / 2;
const eeD = t => A1 + 3 * A2 * t ** 2 + 7 * A3 * t ** 6 + 9 * A4 * t ** 8;

export function project(view, lat, lon) {
  if (view === 'plate_carree') return { x: lon * D, y: lat * D };
  const t = Math.asin(M * Math.sin(lat * D));
  return { x: (2 * Math.sqrt(3) * lon * D * Math.cos(t)) / (3 * eeD(t)), y: A1 * t + A2 * t ** 3 + A3 * t ** 7 + A4 * t ** 9 };
}

export function unproject(view, x, y) {
  if (view === 'plate_carree') {
    if (Math.abs(x) > PC_XMAX || Math.abs(y) > PC_YMAX) return null;
    return { lat: y / D, lon: x / D };
  }
  if (Math.abs(y) > EE_YMAX) return null;
  let t = y / A1;
  for (let i = 0; i < 12; i++) t -= (A1 * t + A2 * t ** 3 + A3 * t ** 7 + A4 * t ** 9 - y) / eeD(t);
  const s = Math.sin(t) / M;
  if (Math.abs(s) > 1 + 1e-12) return null;
  const lon = (3 * x * eeD(t)) / (2 * Math.sqrt(3) * Math.cos(t)) / D;
  if (!Number.isFinite(lon) || Math.abs(lon) > 180 + 1e-9) return null;
  return { lat: Math.asin(Math.max(-1, Math.min(1, s))) / D, lon: Math.max(-180, Math.min(180, lon)) };
}

export function sphereMesh(latSegs, lonSegs) {
  const nv = (latSegs + 1) * (lonSegs + 1);
  const positions = new Float32Array(nv * 3), uvs = new Float32Array(nv * 2);
  for (let i = 0; i <= latSegs; i++) {
    const lat = 90 - (180 * i) / latSegs;
    for (let j = 0; j <= lonSegs; j++) {
      const k = i * (lonSegs + 1) + j;
      positions.set(toVec(lat, -180 + (360 * j) / lonSegs), 3 * k);
      uvs[2 * k] = j / lonSegs;
      uvs[2 * k + 1] = 1 - i / latSegs;
    }
  }
  const indices = new Uint32Array(latSegs * lonSegs * 6);
  let n = 0;
  for (let i = 0; i < latSegs; i++) {
    for (let j = 0; j < lonSegs; j++) {
      const a = i * (lonSegs + 1) + j, b = a + lonSegs + 1, c = b + 1, d = a + 1;
      indices.set([a, b, c, a, c, d], n);
      n += 6;
    }
  }
  return { positions, uvs, indices };
}