Aislo/D03_PlanAnalysis/logic/D03_PlanAnalysis_Tm.ts

/* =============================================================================
 * D03_PlanAnalysis_Tm.ts
 * 작업 좌표계(TM m) ↔ 경위도 — 한국 GRS80 TM(EPSG:5185~5188 · 5179) · 화면 한 장 TM 틀(PLAN 586-0).
 *   짝 = pyproj(`common_util_crs.tm_to_lonlat`) · 거울 시험 `test_586_tm_frame.py`.
 *   바탕(위성 · 도엽 겹)은 경위도로 펴고 그은 선 · 급경사 그림은 TM 이라, 두 상자를 같은 직사각형으로 보면
 *   자오선 수렴각만큼 어긋남(D 프로젝트 귀퉁이 63 m) → TM 점을 경위도로 바꿔 바탕과 같은 자리에 놓음.
 *   DOM · fetch · sessionStorage 없음.
 * ========================================================================== */

export interface TmParams {
  lon0: number;
  lat0: number;
  k0: number;
  fe: number;
  fn: number;
}

const KOREA_TM: Record<number, TmParams> = {
  5185: { lon0: 125, lat0: 38, k0: 1, fe: 200000, fn: 600000 },
  5186: { lon0: 127, lat0: 38, k0: 1, fe: 200000, fn: 600000 },
  5187: { lon0: 129, lat0: 38, k0: 1, fe: 200000, fn: 600000 },
  5188: { lon0: 131, lat0: 38, k0: 1, fe: 200000, fn: 600000 },
  5179: { lon0: 127.5, lat0: 38, k0: 0.9996, fe: 1000000, fn: 2000000 },
};

/** 아는 좌표계만(GRS80 TM) — 그 밖(옛 Bessel 등) = null. */
export const tmParams = (epsg: number | null | undefined): TmParams | null =>
  (epsg && KOREA_TM[epsg]) || null;

// GRS80 · Krüger 6차(Karney 2011) — 수 km 안 mm 아래
const F = 1 / 298.257222101;
const N = F / (2 - F);
const E = Math.sqrt(F * (2 - F));
const A = (6378137 / (1 + N)) * (1 + N ** 2 / 4 + N ** 4 / 64 + N ** 6 / 256);
const ALPHA = [
  N / 2 - (2 * N ** 2) / 3 + (5 * N ** 3) / 16 + (41 * N ** 4) / 180 - (127 * N ** 5) / 288,
  (13 * N ** 2) / 48 - (3 * N ** 3) / 5 + (557 * N ** 4) / 1440 + (281 * N ** 5) / 630,
  (61 * N ** 3) / 240 - (103 * N ** 4) / 140 + (15061 * N ** 5) / 26880,
  (49561 * N ** 4) / 161280 - (179 * N ** 5) / 168,
  (34729 * N ** 5) / 80640,
];
const BETA = [
  N / 2 - (2 * N ** 2) / 3 + (37 * N ** 3) / 96 - N ** 4 / 360 - (81 * N ** 5) / 512,
  N ** 2 / 48 + N ** 3 / 15 - (437 * N ** 4) / 1440 + (46 * N ** 5) / 105,
  (17 * N ** 3) / 480 - (37 * N ** 4) / 840 - (209 * N ** 5) / 4480,
  (4397 * N ** 4) / 161280 - (11 * N ** 5) / 504,
  (4583 * N ** 5) / 161280,
];
const DELTA = [
  2 * N - (2 * N ** 2) / 3 - 2 * N ** 3 + (116 * N ** 4) / 45,
  (7 * N ** 2) / 3 - (8 * N ** 3) / 5 - (227 * N ** 4) / 45,
  (56 * N ** 3) / 15 - (136 * N ** 4) / 35,
  (4279 * N ** 4) / 630,
];
const RAD = Math.PI / 180;

/** 경위도 → 등각 (ξ, η). */
function conformal(lonRad: number, latRad: number): [number, number] {
  const s = Math.sin(latRad);
  const t = Math.sinh(Math.atanh(s) - E * Math.atanh(E * s));
  const xi1 = Math.atan2(t, Math.cos(lonRad));
  const eta1 = Math.atanh(Math.sin(lonRad) / Math.sqrt(1 + t * t));
  let xi = xi1;
  let eta = eta1;
  ALPHA.forEach((a, i) => {
    const j = 2 * (i + 1);
    xi += a * Math.sin(j * xi1) * Math.cosh(j * eta1);
    eta += a * Math.cos(j * xi1) * Math.sinh(j * eta1);
  });
  return [xi, eta];
}

const xi0Of = (p: TmParams): number => conformal(0, p.lat0 * RAD)[0];

export function lonLatToTm(p: TmParams, lon: number, lat: number): [number, number] {
  const [xi, eta] = conformal((lon - p.lon0) * RAD, lat * RAD);
  return [p.fe + p.k0 * A * eta, p.fn + p.k0 * A * (xi - xi0Of(p))];
}

export function tmToLonLat(p: TmParams, x: number, y: number): [number, number] {
  const xi = (y - p.fn) / (p.k0 * A) + xi0Of(p);
  const eta = (x - p.fe) / (p.k0 * A);
  let xi1 = xi;
  let eta1 = eta;
  BETA.forEach((b, i) => {
    const j = 2 * (i + 1);
    xi1 -= b * Math.sin(j * xi) * Math.cosh(j * eta);
    eta1 -= b * Math.cos(j * xi) * Math.sinh(j * eta);
  });
  const chi = Math.asin(Math.sin(xi1) / Math.cosh(eta1));
  let lat = chi;
  DELTA.forEach((d, i) => (lat += d * Math.sin(2 * (i + 1) * chi)));
  return [p.lon0 + Math.atan2(Math.sinh(eta1), Math.cos(xi1)) / RAD, lat / RAD];
}

/** 화면 한 장 TM 틀 — 화면 가운데 둘레 1차 맞춤(그 화면 안 1 px 아래 · 그림 겹은 세 귀퉁이로 놓음). */
export interface TmFrame {
  toScreen: (x: number, y: number) => [number, number];
  toMetric: (px: number, py: number) => [number, number];
  pxPerMeter: number;
}

/** 미분 걸음(m) — 가운데 차분. */
const STEP_M = 25;

/**
 * 바탕이 경위도를 화면에 펴는 길(`project` · `unproject`)을 받아, 화면 가운데(`cx`, `cy`) 둘레에서
 * TM → 화면을 1차식으로 맞춤. 바탕이 어떤 꼴(캔버스 직사각형 · MapLibre 메르카토르)이든 같은 식.
 */
export function localTmFrame(
  p: TmParams,
  project: (lon: number, lat: number) => [number, number],
  unproject: (px: number, py: number) => [number, number],
  cx: number,
  cy: number,
): TmFrame {
  const [mx, my] = lonLatToTm(p, ...unproject(cx, cy));
  const at = (x: number, y: number): [number, number] => project(...tmToLonLat(p, x, y));
  const [ox, oy] = at(mx, my);
  const [x1, y1] = at(mx + STEP_M, my);
  const [x0, y0] = at(mx - STEP_M, my);
  const [x3, y3] = at(mx, my + STEP_M);
  const [x2, y2] = at(mx, my - STEP_M);
  const a = (x1 - x0) / (2 * STEP_M);
  const b = (y1 - y0) / (2 * STEP_M);
  const c = (x3 - x2) / (2 * STEP_M);
  const d = (y3 - y2) / (2 * STEP_M);
  const det = a * d - b * c || 1;
  return {
    toScreen: (x, y) => [ox + a * (x - mx) + c * (y - my), oy + b * (x - mx) + d * (y - my)],
    toMetric: (px, py) => {
      const dx = px - ox;
      const dy = py - oy;
      return [mx + (d * dx - c * dy) / det, my + (a * dy - b * dx) / det];
    },
    pxPerMeter: Math.sqrt(Math.abs(det)),
  };
}