緯度経度からMGRSへの変換

以前の記事で、メートル座標への変換を書きましたが、今回はMGRSへの変換になります。

https://qiita.com/PonSuke9/items/22d038a5d213f9391b6e

## MGRSとは（ChatGPTより）
MGRS（Military Grid Reference System）は、軍用グリッド参照システムの略で、地図上の位置を正確に特定するための座標系の一つです。主に軍事や測量で使用されることが多いですが、その他の地理的な用途にも使われることがあります。

MGRSは、次の要素で構成されます：

- グリッドゾーン識別子（緯度バンドと経度ゾーン）
世界を60の経度帯（各帯は6度の経度幅）に分け、各帯をAからZまでの文字で識別します。また、北緯と南緯を表すためにさらに細かく分割され、例えば「33T」などの形式になります。

- 100km識別子
各グリッドゾーン内をさらに100km四方のグリッドに分け、その各グリッドを2つの英字で識別します。

- 東方向座標と北方向座標
100kmグリッド内の具体的な位置を示すために、通常5桁の数値で表されます。この座標は具体的な地点のメートル単位の位置を示します。

例えば、MGRS座標「33TWN882721」では、次のようになります：

- 33T：グリッドゾーン識別子
- WN：100,000メートル識別子
- 882：東方向座標
- 721：北方向座標

このようにして、MGRSは地球上の位置を正確に示すことができます。

## 緯度経度からUTMへの変換

MGRSへ変換するためには、まずUTM座標への変換が必要になります。
長半径や離心率はメートル座標への変換と同じです。

### 変換用定数

<details open><summary>数式</summary>

- 長半径と逆扁平率

||長半径($a$)|逆扁平率($F$)|
|:-:|:-:|:-:|
|GRS80|6378137|298.257222101|
|WGS-84|6378137|298.257223563|

- 第一離心率

```math
e=\frac{\sqrt{2F-1}}{F}\\
```

- 平坦化定数

```math
n=\frac{1}{2F-1}\\
```

- 子午線半径

```math
A = \frac{a}{1+n} \left( 1 + \frac{n^2}{2^2} + \frac{n^4}{2^6} + \frac{n^6}{2^8} \right)
```

- クルーガー式

```math
\begin{aligned}
\alpha _1 &=\frac{1}{2}n
          - \frac{2}{3}n^2
          + \frac{5}{16}n^3
          + \frac{41}{180}n^4
          - \frac{127}{288}n^5
          + \frac{7891}{37800}n^6
\\\\
\alpha _2 &=\frac{13}{48}n^2
          - \frac{3}{5}n^3
          + \frac{557}{1440}n^4
          + \frac{281}{630}n^5
          - \frac{1983433}{1935360}n^6
\\\\
\alpha _3 &=\frac{61}{240}n^3
          - \frac{103}{140}n^4
          + \frac{15061}{26880}n^5
          + \frac{167603}{181440}n^6
\\\\
\alpha _4 &=\frac{49561}{161280}n^4
          - \frac{179}{168}n^5
          + \frac{6601661}{7257600}n^6
\\\\
\alpha _5 &=\frac{34729}{80640}n^5
          - \frac{3418889}{1995840}n^5
\\\\
\alpha _6 &=\frac{212378941}{319334400}n^6
\end{aligned}
```

- UTMスケール

```math
k=0.9996
```

</details>
<details open><summary>コード</summary>

```cs
/// <summary>長半径[GRS80,WGS-84]</summary>
const double a = 6378137;

/// <summary>逆扁平率[GRS80]</summary>
const double F = 298.257222101;

/// <summary>第一離心率</summary>
static readonly double e = Math.Sqrt(2 * F - 1) / F;

/// <summary>平坦化定数</summary>
const double n = 1 / (2 * F - 1), n2 = n * n, n3 = n2 * n, n4 = n3 * n, n5 = n4 * n, n6 = n5 * n;

/// <summary>子午線半径</summary>
const double A = a / (1 + n) * (1 + n2 / 4 + n4 / 64 + n6 / 256);

/// <summary>6次のクルーガー式</summary>
static readonly double[] α = [
    0,
    1 / 2d * n - 2 / 3d * n2 + 5 / 16d * n3 + 41 / 180d * n4 - 127 / 288d * n5 + 7891 / 37800d * n6,
    13 / 48d * n2 - 3 / 5d * n3 + 557 / 1440d * n4 + 281 / 630d * n5 - 1983433 / 1935360d * n6,
    61 / 240d * n3 - 103 / 140d * n4 + 15061 / 26880d * n5 + 167603 / 181440d * n6,
    49561 / 161280d * n4 - 179 / 168d * n5 + 6601661 / 7257600d * n6,
    34729 / 80640d * n5 - 3418889 / 1995840d * n6,
    212378941 / 319334400d * n6,
];
```

</details>

### 入力数

- 緯度$(degree)$

```math
latitude（南緯-,北緯+）
```

- 経度$(degree)$

```math
longitude（西経-,東経+）
```

### 変換式

<details open><summary>数式</summary>

- 緯度バンド

```math
\begin{aligned}
band = &B \left[ \frac{latitude+180}{8} \right]
\left(-80 \leq latitude \leq 84 \right)\\
&B = \left[C,D,E,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,W,X \right]
\end{aligned}
```

- 経度ゾーン

```math
zone = \frac{longitude+180}{6}+1
```

- ゾーンの調整

|$band$|$zone$|$longitude$|$zone調整$|
|:-:|:-:|:-:|:-:|
|$V$|$31$|$\geq3$ |$+1$|
|$X$|$32$|$<9$    |$-1$|
|^  |^   |$\geq9$ |$+1$|
|^  |$34$|$<21$   |$-1$|
|^  |^   |$\geq21$|$+1$|
|^  |$36$|$<33$   |$-1$|
|^  |^   |$\geq33$|$+1$|

- 基準子午線$(degree)$

```math
lon0 = 6(zone - 1) + 3 - 180
```

- 差分経度$(radian)$

```math
\lambda = \frac{longitude - lon0}{180} \pi
```

- 緯度$(radian)$

```math
\phi = \frac{latitude}{180} \pi
```

- 途中計算

```math
\begin{aligned}
\tau &= \tan \phi\\
\sigma &= \sinh \left(e \tanh^{-1} \left(e \frac{\tau}{\sqrt{1+\tau^2}} \right) \right)\\
\\
\tau' &= \tau \sqrt{1 + \sigma^2} - \sigma \sqrt{1 + \tau^2}\\
\xi' &= \tan^{-1} \frac{\cos \lambda}{\tau} \left(-\pi < \xi' < \pi \right)\\
\eta' &= \sinh^{-1} \left(\frac{\sin \lambda}{\sqrt{\tau'^2 + \cos^2 \lambda}} \right)\\
\\
\xi &=\xi' + \sum_{j=1}^6 \left(\alpha_j \cdot \sin{2 j \xi'} \cdot \cosh{2 j \eta'} \right)\\
\eta &=\eta' + \sum_{j=1}^6 \left(\alpha_j \cdot \cos{2 j \xi'} \cdot \sinh{2 j \eta'} \right)\\
\end{aligned}
```

- X座標

```math
x = k0 \cdot A \cdot \eta + 500,000
```

- Y座標

```math
\begin{aligned}
y &= k0 \cdot A \cdot \xi\\
(y < 0:y &= y + 10,000,000)
\end{aligned}
```

UTMでは、$zone,x,y$の3つの数値で位置を表します。

</details>
<details open><summary>コード</summary>

```cs
public struct DegRad<T> where T : IBinaryFloatingPointIeee754<T>
{
    static readonly T DEG2RAD = T.Pi / T.CreateChecked(180);
    public T Deg { get; set; }
    public T Rad { readonly get => Deg * DEG2RAD; set => Deg = value / DEG2RAD; }

    /// <summary>角度のTanを返します</summary>
    public readonly T Tan() => T.Tan(Rad);

    /// <summary>角度のSinCosを返します</summary>
    public readonly (T Sin, T Cos) SinCos() => T.SinCos(Rad);
}

public record Geocode(DegRad<double> Latitude, DegRad<double> Longitude);
```

```cs
/// <summary>グリッド</summary>
const int LO_GRID = 6, LA_GRID = 8;
const string LA_BAND = "CDEFGHJKLMNPQRSTUVWX";

/// <param name="lon">経度(Degree)</param>
static int UTMZone(DegRad<double> lon) => (int)Math.Floor((lon.Deg + 180) / LO_GRID) + 1;

/// <param name="lat">緯度(Degree)</param>
static char UTMBand(DegRad<double> lat)
{
    if (lat.Deg is < -80 or > 84) return UTM.NaN.Band;
    return LA_BAND[Math.Max(0, Math.Min((int)Math.Floor(lat.Deg / LA_GRID) + 10, LA_BAND.Length - 1))];
}

public record UTM(int Zone, char Band, double X, double Y)
{
    public static readonly UTM NaN = new(0, (char)0, double.NaN, double.NaN);

    public bool IsNaN => double.IsNaN(X) || double.IsNaN(Y);
}

/// <summary>指定したGPS座標をUTM座標に変換して返します。</summary>
public static UTM ConvertUTM(this Geocode geo)
{
    var band = UTMBand(geo.Latitude);
    if (band == UTM.NaN.Band) return UTM.NaN;

    var zone = UTMZone(geo.Longitude);

    // ノルウェーのゾーン調整
    if (band == 'V' && zone == 31 && geo.Longitude.Deg >= 3) zone++;
    // スバールバル諸島のゾーン調整
    if (band == 'X')
        zone += zone switch
        {
            32 => geo.Longitude.Deg < 9 ? -1 : 1,
            34 => geo.Longitude.Deg < 21 ? -1 : 1,
            36 => geo.Longitude.Deg < 33 ? -1 : 1,
            _ => 0,
        };

    var lon0 = (zone - 1) * LO_GRID - 180 + LO_GRID / 2;

    var λ = new DegRad<double> { Deg = geo.Longitude.Deg - lon0 };
    var (sinλ, cosλ) = λ.SinCos();

    var τ = geo.Latitude.Tan();
    var σ = Math.Sinh(e * Math.Atanh(e * τ / Math.Sqrt(1 + τ * τ)));

    // (ʹ)は、共形球体上の角度を示します
    var τʹ = τ * Math.Sqrt(1 + σ * σ) - σ * Math.Sqrt(1 + τ * τ);

    var ξʹ = Math.Atan2(τʹ, cosλ);
    var ηʹ = Math.Asinh(sinλ / Math.Sqrt(τʹ * τʹ + cosλ * cosλ));

    double ξ = ξʹ, η = ηʹ;
    for (var j = 1; j < α.Length; j++)
    {
        ξ += α[j] * Math.Sin(2 * j * ξʹ) * Math.Cosh(2 * j * ηʹ);
        η += α[j] * Math.Cos(2 * j * ξʹ) * Math.Sinh(2 * j * ηʹ);
    }

    var x = k0 * A * η + 500e3;
    var y = k0 * A * ξ + (ξ < 0 ? 10000e3 : 0);

    return new UTM(zone, band, Math.Round(x, 3), Math.Round(y, 3));
}
```

</details>

## UTMからMGRSへの変換

UTM変換で算出した$band$と$zone$はそのまま使います。
$x,y$座標を更に変換して、MGRSを作成します。

<details open><summary>数式</summary>

### 100km格子識別子

- 緯度方向

```math
\begin{aligned}
LaID_0 &= [A,B,C,D,E,F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V]\\
LaID_1 &= [F,G,H,J,K,L,M,N,P,Q,R,S,T,U,V,A,B,C,D,E]
\end{aligned}
```

- 緯度方向

```math
\begin{aligned}
LoID_0 &= [A,B,C,D,E,F,G,H]\\
LoID_1 &= [J,K,L,M,N,P,Q,R]\\
LoID_2 &= [S,T,U,V,W,X,Y,Z]
\end{aligned}
```

### 変換式

- 緯度方向格子識別子

```math
LaID = LaID_{(zone-1)\%2} \left[\frac{y}{100,000}\%20 \right]
```

- 経度方向格子識別子

```math
LoID = LoID_{(zone-1)\%3} \left[\frac{x}{100,000}-1 \right]
```

- 表示桁制御

```math
m = 10^{5-i}:(i=表示桁数)
```

- 東方向座標

```math
easting = \frac{x\%100,000}{m}
```

- 北方向座標

```math
northing = \frac{y\%100,000}{m}
```

- MGRS文字列

```math
\{zone\}\{band\}\{LaID\}\{LoID\}\{easting\}\{northing\}
```

</details>
<details open><summary>コード</summary>

```cs
static readonly string[] LaID = ["ABCDEFGHJKLMNPQRSTUV", "FGHJKLMNPQRSTUVABCDE"];
static readonly string[] LoID = ["ABCDEFGH", "JKLMNPQR", "STUVWXYZ"];

/// <summary>指定したGPS座標をMGRS文字列に変換して返します。</summary>
/// <param name="digits">距離文字数</param>
public static string ConvertMGRS(this Geocode geo, int digits = 5) =>
    geo.ConvertUTM().ConvertMGRS(digits);

/// <summary>指定したGPS座標をMGRS文字列に変換して返します。</summary>
/// <param name="digits">距離文字数</param>
public static string ConvertMGRS(this UTM utm, int digits = 5)
{
    if (utm == null || utm.IsNaN) return string.Empty;

    var xID = LoID[(utm.Zone - 1) % 3][(int)Math.Floor(utm.X / 100e3) - 1];
    var yID = LaID[(utm.Zone - 1) % 2][(int)Math.Floor(utm.Y / 100e3) % 20];

    var mul = Math.Pow(10, 5 - digits);
    var padding = new string('0', digits);
    var easting = Math.Floor(utm.X % 100e3 / mul).ToString(padding);
    var northing = Math.Floor(utm.Y % 100e3 / mul).ToString(padding);

    return $"{utm.Zone:00}{utm.Band}{xID}{yID}{easting}{northing}";
}
```

</details>