Include computation of scale factor (k) and convergence angle (γ)

[yongjun21]

I used this library not for lat/lon to UTM x/y conversion but for deriving the transformation from UTM x/y offset to Web Mercator x/y offset from a ref point.

Current process involves

1. converting 3 nearby points in UTM to lat/lon then
2. calculating the Web Mercator meter offsets of the 3 transformed points then
3. solving the linear equations to get the transformation that maps UTM offsets (from a ref point) to Web Mercator offsets
4. the result is only an approximation since it relies on solving the transformation from the empirical values of 3 transformed points

After digging into the theory behind UTM, I realize this can be done more simply and accurately but computing the convergence angle and scale factor for the ref point: https://en.wikipedia.org/wiki/Transverse_Mercator_projection#convergence.

Like most other UTM libs, this does not include functions to compute k and γ but is probably not hard to extend since most of the logic for computing preliminary and intermediate values are already inside.

If maintainer is not able to commit time, I might be able to create a PR for this

[bartvanandel]

For completeness, Wikipedia links for both formulas:

• Point scale ($k$)
• Convergence ($γ$)

If you're still interested, a PR would be nice!

[bartvanandel]

I've looked into this, and have implemented the algorithm found on that Wikipedia page, which indeed includes convergence and scale factor, and compared it to the currently implemented algorithm, for which I still haven't found the actual source (although I did find some articles mentioning the same constants being used). Unfortunately, the Wikipedia algorithm is about twice as slow, from the top of my head. Accuracy might be slightly better. I haven't decided yet whether to offer this algorithm or what to do with it.

Unfortunately, your assumption that it's "probably not hard to extend" the existing algorithm isn't true as the two algorithms take quite a different approach. Unless I'm missing something, that is.

[yongjun-datature]

Is it this one: https://en.wikipedia.org/wiki/Transverse_Mercator:_Redfearn_series you referring to?

[yongjun-datature]

Unfortunately, I have moved on from the project where I've previously needed this feature. So I do not think I will be looking at this anymore.

But back then I have gotten an intern to help me implement this in JS:

import { TO_DEGREE } from './math';

/*
UTM bases itself on WGS 84. Hence, we reference the Semi-major axis and Semi-minor axis of WGS84.
Args:
    x: Easting from UTM
    y: Northing from UTM
Returns:
    longitude, latitude, scaling and convergence (will create object here upon confirmation where you can get what you want via object.longitude, etc.)
Variables:
    Semi Major Axis = a
    Semi Minor Axis = b
    Flattening = f
    Third Flattening = n
    Rectifying Latitude for Footpoint Latitude = mu
    Meridian distance from equator to pole = mp
    Central Meridian Scale Factor = k0
    Footpoint Latitude = phi_1
    First Eccentricity Squared = e_squared
*/

// Semi Major Axis
const a = 6378137;

// Semi Minor Axis
const b = 6356752.314245;

// Flattening
const f = (a - b) / a;

// First eccentricity squared
const e_squared = 2 * f - f ** 2;

// Third Flattening
const n = (a - b) / (a + b);

// Caching reused n values
const n_2 = n * n;
const n_3 = n * n_2;
const n_4 = n * n_3;

// Meridian distance from equator to pole
const b0 = b * (1 + n + (5 / 4) * n_2 + (5 / 4) * n_3);
const mp = (Math.PI * b0) / 2;

// Central Meridian Scale Factor
const k0 = 0.9996;

// Footpoint Latitude
const D2 = (3 / 2) * n - (27 / 32) * n_3;
const D4 = (21 / 16) * n_2 - (55 / 32) * n_4;
const D6 = (151 / 96) * n_3;
const D8 = (1097 / 512) * n_4;

export function utmToLatLon(
  x: number,
  y: number,
  zoneNumber: number,
  northern: boolean
) {
  x -= 500000;

  if (!northern) {
    y -= 10000000;
  }

  // Rectifying Latitude using Footpoint Latitude
  const mu = (Math.PI * y) / (2 * mp * k0);

  const phi_1 =
    mu +
    D2 * Math.sin(2 * mu) +
    D4 * Math.sin(4 * mu) +
    D6 * Math.sin(6 * mu) +
    D8 * Math.sin(8 * mu);

  const t1 = Math.tan(phi_1);
  const v1 = a / Math.sqrt(1 - e_squared * Math.sin(phi_1) ** 2);
  const p1 = (v1 ** 3 * (1 - e_squared)) / a ** 2;
  const b1 = v1 / p1;
  const c1 = Math.cos(phi_1);

  // Caching reused values
  const t1_2 = t1 * t1;
  const t1_4 = t1_2 * t1_2;
  const t1_6 = t1_4 * t1_2;

  const b1_2 = b1 * b1;
  const b1_3 = b1_2 * b1;
  const b1_4 = b1_2 * b1_2;

  const x_2 = x * x;
  const x_3 = x * x_2;
  const x_4 = x * x_3;
  const x_5 = x * x_4;
  const x_6 = x * x_5;
  const x_7 = x * x_6;
  const x_8 = x * x_7;

  const k0_v1 = k0 * v1;
  const k0_v1_2 = k0_v1 * k0_v1;
  const k0_v1_3 = k0_v1 * k0_v1_2;
  const k0_v1_4 = k0_v1 * k0_v1_3;
  const k0_v1_5 = k0_v1 * k0_v1_4;
  const k0_v1_6 = k0_v1 * k0_v1_5;
  const k0_v1_7 = k0_v1 * k0_v1_6;
  const k0_v1_8 = k0_v1 * k0_v1_7;

  // Calculate Longitude
  const V3 = b1 + 2 * t1_2;
  const V5 =
    4 * b1_3 * (1 - 6 * t1_2) -
    b1_2 * (9 - 68 * t1_2) -
    72 * b1 * t1_2 -
    24 * t1_4;
  const V7 = 61 + 662 * t1_2 + 1320 * t1_4 + 720 * t1_6;
  let longitude =
    x / (c1 * k0_v1) -
    (x_3 * V3) / (6 * c1 * k0_v1_3) -
    (x_5 * V5) / (120 * c1 * k0_v1_5) -
    (x_7 * V7) / (5040 * c1 * k0_v1_7);
  longitude = longitude * TO_DEGREE - 183 + zoneNumber * 6;

  // Calculate Latitude
  const U4 = 4 * b1_2 - 9 * b1 * (1 - t1_2) - 12 * t1_2;
  const U6 =
    8 * b1_4 * (11 - 24 * t1_2) -
    12 * b1_3 * (21 - 71 * t1_2) +
    15 * b1_2 * (15 - 98 * t1_2 + 15 * t1_4) +
    180 * b1 * (5 * t1_2 - 3 * t1_4) +
    360 * t1_4;
  const U8 = -1385 - 3633 * t1_2 - 4095 * t1_4 - 1575 * t1_6;
  let latitude =
    phi_1 -
    (x_2 * b1 * t1) / (2 * k0_v1_2) -
    (x_4 * b1 * t1 * U4) / (24 * k0_v1_4) -
    (x_6 * b1 * t1 * U6) / (720 * k0_v1_6) -
    (x_8 * b1 * t1 * U8) / (40320 * k0_v1_8);
  latitude = latitude * TO_DEGREE;

  // Calculate Scaling
  const K2 = b1;
  const K4 =
    4 * b1_3 * (1 - 6 * t1_2) - 3 * b1_2 * (1 - 16 * t1_2) - 24 * b1 * t1_2;
  const K6 = 1;
  const scaleFactor =
    k0 *
    (1 +
      (x_2 * K2) / (2 * k0_v1_2) +
      (x_4 * K4) / (24 * k0_v1_4) +
      (x_6 * K6) / (720 * k0_v1_6));

  // Calculate Convergence
  const K3 = 2 * b1_2 - 3 * b1 - t1_2;
  const K5 =
    b1_4 * (11 - 24 * t1_2) -
    3 * b1_3 * (8 - 23 * t1_2) +
    5 * b1_2 * (3 - 14 * t1_2) +
    30 * b1 * t1_2 +
    3 * t1_4;
  const K7 = -17 - 77 * t1_2 - 105 * t1_4 - 45 * t1_6;
  const convergence =
    (x * t1) / k0_v1 +
    (x_3 * t1 * K3) / (3 * k0_v1_3) +
    (x_5 * t1 * K5) / (15 * k0_v1_5) +
    (x_7 * t1 * K7) / (315 * k0_v1_7);

  return { latitude, longitude, scaleFactor, convergence };
}

I'm pretty sure this is the same approach taken by the library source so extending shouldn't be too hard.

[bartvanandel]

Thanks! I'll have a look later.