Map layers: web Mercator, tile grids and world files

// ------------------------------------------------------------------ web Mercator
export function lonLatToMerc(lon, lat) {
  const l = Math.max(-85.0511287798066, Math.min(85.0511287798066, lat));
  return [R_MERC * lon * Math.PI / 180, R_MERC * Math.log(Math.tan(Math.PI / 4 + l * Math.PI / 360))];
}

export function mercToLonLat(x, y) {
  return [x / R_MERC * 180 / Math.PI, (2 * Math.atan(Math.exp(y / R_MERC)) - Math.PI / 2) * 180 / Math.PI];
}

/** The standard XYZ ("Google") tile grid: 256 px tiles, zoom 0 = one tile for the world. */
export const XYZ_GRID = Object.freeze({
  code: 3857, origin: [-HALF_MERC, HALF_MERC],
  matrices: Array.from({ length: 23 }, (_, z) => Object.freeze({ id: String(z), res: (2 * HALF_MERC) / 256 / 2 ** z, tw: 256, th: 256, mw: 2 ** z, mh: 2 ** z })),
});

/** Tile column/row of a lon/lat at zoom z (whole numbers). */
export function tileAt(lon, lat, z) {
  const [x, y] = lonLatToMerc(lon, lat), m = XYZ_GRID.matrices[z];
  return [Math.floor((x + HALF_MERC) / (m.res * 256)), Math.floor((HALF_MERC - y) / (m.res * 256))];
}

/** Tiles of `grid` covering bbox [minX, minY, maxX, maxY] (grid units) at screen resolution `res`; a coarser matrix when
 *  more than `cap` tiles would be needed. */
export function gridTiles(grid, bbox, res, { cap = 64 } = {}) {
  const ms = [...grid.matrices].sort((a, b) => b.res - a.res);
  let i = ms.indexOf(pickMatrix(grid, res));
  for (; i >= 0; i--) {
    const m = ms[i], tw = m.res * m.tw, th = m.res * m.th;
    let c0 = Math.floor((bbox[0] - grid.origin[0]) / tw), c1 = Math.floor((bbox[2] - grid.origin[0]) / tw);
    let r0 = Math.floor((grid.origin[1] - bbox[3]) / th), r1 = Math.floor((grid.origin[1] - bbox[1]) / th);
    c0 = Math.max(0, c0); r0 = Math.max(0, r0);
    if (Number.isFinite(m.mw)) c1 = Math.min(m.mw - 1, c1);
    if (Number.isFinite(m.mh)) r1 = Math.min(m.mh - 1, r1);
    if (!(c1 >= c0 && r1 >= r0)) return { matrix: m, tiles: [] };
    if ((c1 - c0 + 1) * (r1 - r0 + 1) > cap && i > 0) continue;
    const tiles = [];
    for (let r = r0; r <= r1; r++) for (let c = c0; c <= c1; c++) {
      if (tiles.length >= cap) break;
      tiles.push({ col: c, row: r, bounds: [grid.origin[0] + c * tw, grid.origin[1] - (r + 1) * th, grid.origin[0] + (c + 1) * tw, grid.origin[1] - r * th] });
    }
    // the middle first: what the eye looks at arrives first
    const cx = (bbox[0] + bbox[2]) / 2, cy = (bbox[1] + bbox[3]) / 2;
    tiles.sort((a, b) => Math.hypot((a.bounds[0] + a.bounds[2]) / 2 - cx, (a.bounds[1] + a.bounds[3]) / 2 - cy) - Math.hypot((b.bounds[0] + b.bounds[2]) / 2 - cx, (b.bounds[1] + b.bounds[3]) / 2 - cy));
    return { matrix: m, tiles };
  }
  return { matrix: ms[0], tiles: [] };
}

/** Metres per unit of a system for WMTS scale denominators (0.28 mm pixels). Geographic: metres per degree. */
export const metresPerUnit = (unit) => (unit === 'deg' ? 2 * Math.PI * 6378137 / 360 : unit === 'ft' ? 0.3048 : unit === 'ftUS' ? 1200 / 3937 : 1);

/** A WMTS tile matrix set as a grid with x/y corners. `toSet(lon, lat)` maps WGS 84 to the set's system and `bbox84`
 *  (the layer's) decides the axis order of TopLeftCorner when the set's own CRS does not say it clearly. */
export function wmtsGrid(set, unit, { toSet = null, bbox84 = null } = {}) {
  const k = metresPerUnit(unit);
  let swap = set.latFirst;
  if (!set.crs84 && !set.latFirst && toSet && bbox84) {
    // projected systems with northing-first axes (NZTM…) write TopLeftCorner northing first: try both, keep the one whose
    // matrix holds the layer's middle
    try {
      const [mx, my] = toSet((bbox84[0] + bbox84[2]) / 2, (bbox84[1] + bbox84[3]) / 2);
      const m = set.matrices[0], w = m.scale * 0.00028 / k * m.tw * m.mw, h = m.scale * 0.00028 / k * m.th * m.mh;
      const holds = ([x, y]) => mx >= x && mx <= x + w && my <= y && my >= y - h;
      if (!holds(m.tl) && holds([m.tl[1], m.tl[0]])) swap = true;
    } catch { /* keep the CRS's own order */ }
  }
  const matrices = set.matrices.map((m) => ({ id: m.id, res: m.scale * 0.00028 / k, tw: m.tw, th: m.th, mw: m.mw, mh: m.mh, origin: swap ? [m.tl[1], m.tl[0]] : [m.tl[0], m.tl[1]] }));
  // one origin per grid; WMTS allows one per matrix — keep each matrix's own
  return { code: set.code, origin: matrices[0].origin, matrices };
}

/** A world file (.pgw .jgw .tfw .wld …): six numbers A D B E C F — x = A·col + B·row + C, y = D·col + E·row + F, for the
 *  CENTRE of pixel (col, row). → GDAL geotransform [x0, dxCol, dxRow, y0, dyCol, dyRow] for pixel CORNERS. */
export function parseWorldFile(text) {
  const s = String(text || '');
  if (s.length > 10000) throw new Error(t('basemap.world.long'));
  const v = s.trim().split(/\s+/).map(Number);
  if (v.length !== 6 || !v.every(Number.isFinite)) throw new Error(t('basemap.world.six'));
  const [A, D, B, E, C, F] = v;
  if (Math.abs(A * E - B * D) < 1e-30) throw new Error(t('basemap.world.zero'));
  return [C - A / 2 - B / 2, A, B, F - D / 2 - E / 2, D, E];
}

export const pixelToSystem = (gt, u, v) => [gt[0] + u * gt[1] + v * gt[2], gt[3] + u * gt[4] + v * gt[5]];

const fmtNum = (v) => (Math.abs(v) >= 1e4 ? v.toFixed(3) : String(+v.toPrecision(12)));