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

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座標への変換が必要になります。
長半径や離心率はメートル座標への変換と同じです。

変換用定数

数式

長半径と逆扁平率

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

第一離心率

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

平坦化定数

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

子午線半径

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

クルーガー式

\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スケール

k=0.9996

コード

/// <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,
];

入力数

緯度$(degree)$

latitude（南緯-,北緯+）

経度$(degree)$

longitude（西経-,東経+）

変換式

数式

緯度バンド

\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}

経度ゾーン

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)$

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

差分経度$(radian)$

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

緯度$(radian)$

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

途中計算

\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座標

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

Y座標

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

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

コード

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);

/// <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));
}

UTMからMGRSへの変換

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

数式

100km格子識別子

緯度方向

\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}

緯度方向

\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}

変換式

緯度方向格子識別子

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

経度方向格子識別子

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

表示桁制御

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

東方向座標

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

北方向座標

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

MGRS文字列

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

コード

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}";
}