測地系について

参考

GISで使う経緯度の基礎知識とかTIPSとか
日本測地系と世界測地系の変換式
日本測地系のミリ秒を世界測地系の度に変換する計算

概要

位置データの表現には以下の要素がある

• 世界測地系／日本測地系
• 10進法表記(度のみ)／60進法表記(度分秒)／ミリ秒表記

ケーススタディ

Q1. 提供された位置データが日本測地系か世界測地系かわからない

Yahoo地図でそれぞれの測地系で経緯度を入力した位置が地図上で見れるので、ズレていない方の測地系データとわかる。
http://maps.loco.yahoo.co.jp/maps?lat=41.7973055555556&lon=140.773833333333&datum=tky
上記URLのリクエストパラメータを変更して確認できる。(詳細はここ)

lat   : 緯度
lon   : 経度
datum : 測地系(wgs=世界測地系／tky=日本測地系)

Google Mapは世界測地系で統一している。

Q2. 二地点の緯度・経度から距離を算出したい

結構複雑な計算が必要。

値	備考
$x_1, y_1$	地点１の経度と緯度
$x_2, y_2$	地点２の経度と緯度
$a = 6,378,137.000$	赤道半径
$b = 6,356,752.314 245$	極半径
$d_x = x_1 - x_2$	経度の差
$d_y = y_1 - y_2$	緯度の差
$\mu_y = \frac{ y_1 + y_2 }{ 2 }$	緯度の平均値
$e = \sqrt{ \frac{ a^2 - b^2 }{ a^2 } }$	第一離心率
$W = \sqrt{ 1 - e^2 \sin^2 \mu_y }$
$N = \frac{ a }{ W }$	卯酉線曲率半径
$M = \frac{ a (1 - e^2) }{ W^3 }$	子午線曲率半径
$d = \sqrt{ (d_y M)^2 + (d_x N \cos \mu_y)^2 }$	距離($m$)

詳しくはこちら。

WGS位置を扱うクラス実装。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
2点間のWGS緯度経度から距離(メートル)に変換
http://yamadarake.jp/trdi/report000001.html
"""

from math import pi, sin, cos, sqrt

# WGS経度緯度クラス
class WgsPoint:

    def __init__(self, kdo, ido):
        self._x = float(kdo)
        self._y = float(ido)

    def getx(self):
        return self._x

    def setx(self, v):
        self._x = float(v)

    def gety(self):
        return self._y

    def sety(self, v):
        self._y = float(v)

    x = property(getx, setx)
    y = property(gety, sety)

# WGS測量クラス
class Measurement:

    r1 = 6378137.000000                              # 赤道半径
    r2 = 6356752.314245                              # 極半径
    e  = sqrt(((r1 ** 2) - (r2 ** 2)) / (r1 ** 2))   # 第一離心率

    def __init__(self, p1, p2):
        if isinstance(p1, WgsPoint) and isinstance(p2, WgsPoint):
            self.points = [p1, p2]
        else:
            raise TypeError('argument type invalid.')

    # ラジアン変換
    def _deg2rad(self, d):
        return d * pi / 180

    # 経度の差
    def dx(self):
        return self._deg2rad(self.points[0].x - self.points[1].x)

    # 緯度の差
    def dy(self):
        return self._deg2rad(self.points[0].y - self.points[1].y)

    # 緯度の平均値
    def my(self):
        return self._deg2rad((self.points[0].y + self.points[1].y) / 2.0)

    def e2(self):
        return self.e ** 2
        #return ((self.r1 ** 2) - (self.r2 ** 2)) / (self.r1 ** 2)

    def W(self):
        return sqrt(1 - self.e2 * (sin(self.my) ** 2))

    # 卯酉線曲率半径
    def M(self):
        return self.r1 * (1 - self.e2) / (self.W ** 3)

    # 子午線曲率半径
    def N(self):
        return self.r1 / self.W

    # 距離
    def distance(self):
        return sqrt((self.dy * self.M) ** 2 + (self.dx * self.N * cos(self.my)) ** 2)

    dx = property(dx)
    dy = property(dy)
    my = property(my)
    e2 = property(e2)
    W  = property(W)
    M  = property(M)
    N  = property(N)
    distance = property(distance)

if __name__ == '__main__':
    tokyo     = WgsPoint(139.74472, 35.65500)
    tsukuba   = WgsPoint(140.09111, 36.10056)
    fukuoka   = WgsPoint(130.36208, 33.59532)
    north     = WgsPoint(140.380034, 35.802739)
    south     = WgsPoint(140.392265, 35.785796)
    hakodate  = WgsPoint(140.726413, 41.773709)
    goryokaku = WgsPoint(140.733539, 41.803557)

    m1 = Measurement(tokyo, tsukuba)
    m2 = Measurement(tokyo, fukuoka)
    m3 = Measurement(north, south)
    m4 = Measurement(hakodate, goryokaku)

    print('東京-筑波                  : ' + str(round(m1.distance / 1000, 3)) + ' km.')
    print('東京-福岡ドーム            : ' + str(round(m2.distance / 1000, 3)) + ' km.')
    print('成田空港北側滑走路北端-南端: ' + str(round(m3.distance / 1000, 3)) + ' km.')
    print('函館-五稜郭                : ' + str(round(m4.distance / 1000, 3)) + ' km.')