Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

お題は不問！Qiita Engineer Festa 2024で記事投稿！
Qiita Engineer Festa20242024年7月17日まで開催中！
@PonSuke9

緯度経度からMGRSへの変換

Last updated at Posted at 2024-07-11

以前の記事で、メートル座標への変換を書きましたが、今回は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}";
}
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?