gendesign/frontend/src/lib/geo-local.ts

/**
 * geo-local.ts — dependency-free geo helpers for the §7 «Концепция» 3D massing
 * viewer.
 *
 * Two jobs:
 *  1. Project WGS84 [lon, lat] → local metres about a parcel centroid
 *     (equirectangular), so building footprints and the OSM ground tile can be
 *     placed in a shared Three.js metre-space and stay geo-aligned.
 *  2. Slippy-tile (OSM XYZ) math to pick + fetch + geo-pin a raster ground tile.
 *
 * No Three.js / DOM imports — pure functions, trivially unit-testable. `R` matches
 * `polygonAreaSqm` in concept-api.ts (WGS84 equatorial radius) for consistency.
 */

import type { Position } from "geojson";

/** WGS84 equatorial radius (m) — same constant as `polygonAreaSqm`. */
export const R = 6_378_137;
const DEG = Math.PI / 180;

// ── bounds / centroid ─────────────────────────────────────────────────────────

export interface Bounds {
  minLon: number;
  minLat: number;
  maxLon: number;
  maxLat: number;
}

/** Axis-aligned WGS84 bounds of a ring (or rings). Null on an empty/degenerate ring. */
export function ringBounds(ring: Position[]): Bounds | null {
  let minLon = Infinity;
  let minLat = Infinity;
  let maxLon = -Infinity;
  let maxLat = -Infinity;
  for (const [lon, lat] of ring) {
    if (lon < minLon) minLon = lon;
    if (lon > maxLon) maxLon = lon;
    if (lat < minLat) minLat = lat;
    if (lat > maxLat) maxLat = lat;
  }
  if (!Number.isFinite(minLon) || !Number.isFinite(minLat)) return null;
  return { minLon, minLat, maxLon, maxLat };
}

/**
 * Ring-average of unique vertices (skips the closing vertex), mirroring
 * `polygonCentroidWkt`. Returns [lon, lat]. Falls back to [0, 0] on an empty ring.
 */
export function ringCentroid(ring: Position[]): [number, number] {
  if (ring.length === 0) return [0, 0];
  // Drop the closing vertex when the ring is closed (first == last).
  const closed =
    ring.length > 1 &&
    ring[0][0] === ring[ring.length - 1][0] &&
    ring[0][1] === ring[ring.length - 1][1];
  const pts = closed ? ring.slice(0, ring.length - 1) : ring;
  if (pts.length === 0) return [0, 0];
  let sLon = 0;
  let sLat = 0;
  for (const [lon, lat] of pts) {
    sLon += lon;
    sLat += lat;
  }
  return [sLon / pts.length, sLat / pts.length];
}

// ── projection (WGS84 → local metres) ─────────────────────────────────────────

export interface LocalPoint {
  /** East offset (metres), east → +x. */
  x: number;
  /** South offset (metres), north → −z (Three.js right-handed). */
  z: number;
}

export type Projector = (pos: Position) => LocalPoint;

/**
 * Equirectangular projector about (lon0, lat0). Sub-metre accurate at lot scale
 * (<1 km). Returns [lon, lat] → { x (east +), z (south +) } in metres.
 *
 * Sign discipline: north → −z, east → +x, height → +y — matches the azimuth
 * convention in `sunFromHour` so the sun arc lands east→south→west correctly.
 */
export function makeProjector(lon0: number, lat0: number): Projector {
  const kx = R * DEG * Math.cos(lat0 * DEG); // metres per degree lon at this lat
  const kz = R * DEG; // metres per degree lat
  return ([lon, lat]: Position) => ({
    x: (lon - lon0) * kx,
    z: -(lat - lat0) * kz,
  });
}

// ── slippy-tile (OSM XYZ) math ────────────────────────────────────────────────

export function lon2tile(lon: number, z: number): number {
  return Math.floor(((lon + 180) / 360) * 2 ** z);
}

export function lat2tile(lat: number, z: number): number {
  const r = lat * DEG;
  return Math.floor(
    ((1 - Math.log(Math.tan(r) + 1 / Math.cos(r)) / Math.PI) / 2) * 2 ** z,
  );
}

export function tile2lon(x: number, z: number): number {
  return (x / 2 ** z) * 360 - 180;
}

export function tile2lat(y: number, z: number): number {
  const n = Math.PI - (2 * Math.PI * y) / 2 ** z;
  return (180 / Math.PI) * Math.atan(0.5 * (Math.exp(n) - Math.exp(-n)));
}

export interface TileRange {
  z: number;
  x0: number;
  x1: number;
  y0: number;
  y1: number;
  cols: number;
  rows: number;
  /** Exact WGS84 span of the stitched tile block (NW + SE corners). */
  westLon: number;
  eastLon: number;
  northLat: number;
  southLat: number;
}

/**
 * Pick the OSM zoom for a bbox such that the stitched tile block stays ≤ maxTiles
 * (default 4 → a 2×2 cap), clamped to [minZ, maxZ]. Starts at maxZ and steps the
 * zoom down until the tile count fits.
 *
 * Returns the integer tile range plus the EXACT WGS84 corners of the stitched
 * block (tile2lon/lat of the inclusive range), which the scene projects with the
 * SAME `makeProjector` as the buildings to geo-pin the ground plane.
 */
export function osmTileRangeForBbox(
  bbox: Bounds,
  opts: { minZ?: number; maxZ?: number; maxTiles?: number } = {},
): TileRange {
  const minZ = opts.minZ ?? 16;
  const maxZ = opts.maxZ ?? 19;
  const maxTiles = opts.maxTiles ?? 4;

  let z = maxZ;
  let x0 = 0;
  let x1 = 0;
  let y0 = 0;
  let y1 = 0;
  for (; z >= minZ; z--) {
    x0 = lon2tile(bbox.minLon, z);
    x1 = lon2tile(bbox.maxLon, z);
    // lat → tile is inverted (north = smaller y), so y0 (north) uses maxLat.
    y0 = lat2tile(bbox.maxLat, z);
    y1 = lat2tile(bbox.minLat, z);
    const cols = x1 - x0 + 1;
    const rows = y1 - y0 + 1;
    if (cols * rows <= maxTiles || z === minZ) break;
  }
  const cols = x1 - x0 + 1;
  const rows = y1 - y0 + 1;
  return {
    z,
    x0,
    x1,
    y0,
    y1,
    cols,
    rows,
    westLon: tile2lon(x0, z),
    eastLon: tile2lon(x1 + 1, z),
    northLat: tile2lat(y0, z),
    southLat: tile2lat(y1 + 1, z),
  };
}

/** Pad a bbox by a fraction of its span on every side (e.g. 0.3 → +30% context). */
export function padBounds(bbox: Bounds, frac: number): Bounds {
  const dLon = (bbox.maxLon - bbox.minLon) * frac;
  const dLat = (bbox.maxLat - bbox.minLat) * frac;
  return {
    minLon: bbox.minLon - dLon,
    minLat: bbox.minLat - dLat,
    maxLon: bbox.maxLon + dLon,
    maxLat: bbox.maxLat + dLat,
  };
}