地図について

測地系の種類

　測地系とは、地球上の位置を経度・緯度によって表すための方式で、地球を楕円体に近似して処理されています。そのため、同じ場所でも測地系の違いにより緯度経度が異なるため、正しく測地系を処理しないと位置がずれてしまうという問題が発生します。
　日本では、主に4種類の測地系が使われていて、近似する楕円体が異なります。また、測量は日本水準原点を基準にして行われてきましたが、天文観測による測量の誤差のため離島毎に水準原点が存在しています。

世界測地系
GRS80、JGD2000とも呼ばれ、国土地理院の新測量図で使われている
ITRF座標系GRS80楕円体による近似
世界的に最も広く使われている測地系
日本測地系
Tokyo Datumとも呼ばれ、国土地理院の旧測量図で使われていた
ベッセル楕円体による近似
測量法で離島毎に水準原点が決められている
GPS用世界測地系
WGS84 Datumとも呼ばれ、GPSなどで使われている
WGS84楕円体による近似
正確にはGRS80、JGD2000と異なるが、同じと考えて問題はない
海域の測地系として標準的に使われている
GPS用日本測地系
楕円体の変換のみを行った測地系で、GPSなどで使われている
Tokyo Datum、日本測地系とも呼ばれ紛らわしいため、GPS用日本測地系と呼んでいる
ベッセル楕円体による近似
離島毎の水準原点に対応していないため、日本測地系の地図で数100mずれることがある

測地系の変換

　測地系の変換には、2つの変換処理が必要となります。 1つ目は近似する楕円体の違いによる変換で、2つ目は水準原点の違いによる変換です。

近似する楕円体の違いによる変換

　近似する楕円体の違いによる変換をすることで、世界測地系（GPS用世界測地系）とGPS用日本測地系の変換を行うことができます。世界測地系（GPS用世界測地系）とGPS用日本測地系では、位置は450mほどずれています。

各測地系が近似に利用している楕円体の値

楕円体の種類	赤道半径	扁平率
GRS80楕円体	6,378,137.000m	1/298.257222
WGS84楕円体	6,378,137.000m	1/298.257224
ベッセル楕円体	6,377,397.155m	1/299.152813

世界測地系（GPS用世界測地系）からGPS用日本測地系への変換式

　世界測地系（GPS用世界測地系）からGPS用日本測地系への変換は、下記の式で行うことができます。三角関数を使わない簡易式ですので、誤差は5m程度あります。

$lat_gtk = $lat_jgd * 1.000106960 - $lng_jgd * 0.000017467 - 0.0046020;
$lng_gtk = $lng_jgd * 1.000083049 + $lat_jgd * 0.000046047 - 0.0100411;

GPS用日本測地系から世界測地系（GPS用世界測地系）への変換式

$lat_jgd = $lat_gtk * 0.999893050 + $lng_gtk * 0.000017464 + 0.0046017;
$lng_jgd = $lng_gtk * 0.999916957 - $lat_gtk * 0.000046038 + 0.0100405;

水準原点の違いによる変換

　水準原点の違いは、測量法によって決まっているだけですので、数学的に変換することはできません。日本測地系とGPS用日本測地系では、下記の表のように780mもずれることがあります。石垣島や与那国島では誤差が大きいのに対し、間にある西表島では誤差が少ないのは、測量を開始した時期の天文測量の精度の違いではないかと思われます。
　日本測地系とGPS用日本測地系の両方に対応していない地図システムでは、これだけのずれが発生する可能性があります。

水準原点の差（日本測地系とGPS用日本測地系の差）

区域	南北誤差(m)	東西誤差(m)
三宅島	-5.8	+4.6
御蔵島	-7.1	5.9
八丈島	-14.5	5.7
青ヶ島	-77.9	1.3
小笠原	-1.5	13.2
硫黄島	783.1	-143.0
隠岐	4.9	-4.3
男女	0.7	4.8
大東	372.9	-505.6
トカラ	5.0	4.9
奄美	7.7	6.5
沖縄	12.6	9.1
尖閣	0.0	-55.7
宮古	56.0	-47.2
多良間	280.3	-462.9
石垣	-147.4	-185.1
西表	7.4	13.3
与那国	-149.9	-178.3

国土地理院の変換パラメータファイル

　国土地理院では、世界測地系と日本測地系の変換を正確に行うために、TKY2JGD.parというパラメータファイルを提供しています。

TKY2JGD.parの形式

　TKY2JGD.parは、MeshCode、dB(sec)、dL(sec)の3レコードからなるファイルです。国土地理院のメッシュは、緯度30秒、経度45秒の大きさのメッシュとなっています。 MeshCodeは、緯度方向1次メッシュ2桁、経度方向1次メッシュ2桁、緯度方向2次メッシュ1桁、経度方向2次メッシュ1桁、緯度方向3次メッシュ1桁、経度方向3次メッシュ1桁となっていますので、それぞれを$mc11、$mc12、$mc21、$mc22、$mc31、$mc32とすると、下記の計算式で緯度経度が求まります。 dB、dLには、緯度、経度方向の補正値（度）が格納されていますので、この値を使って測地系の変換をすることができます。

$lat = ($mc11 + $mc21 / 8 + $mc31 / 64) / 1.5;
$lng = ($mc12 + $mc22 / 8 + $mc32 / 64) + 100;

　水準原点の差（日本測地系とGPS用日本測地系の差）の表は、TKY2JGD.par内のすべてのメッシュコードに対して、近似する楕円体の違いによる測地系変換を行った結果と、TKY2JGD.parの値によって測地系変換を行った結果で、0.2秒以上のずれがある地点を調べて作成してありますので、それなりの精度があります。
　実際に、この変換パラメータを利用して測地系変換を行っていますが、変換誤差は10m以内に収まっています。

緯度経度と距離

各緯度における距離

　経線は北に行くほど狭くなりますので、経度1秒あたりの距離は緯度によって異なります。また、地球は球ではなく楕円体に近いので、緯度1秒あたりの距離も緯度によって異なります。
　そこで、各緯度における1秒あたりの距離と1mあたりの秒数を表にしました。
　誤差を10m以下にする場合は、緯度経度の誤差は0.3秒以下であれば十分です。

緯度	南北距離	東西距離	緯度の長さ	経度の長さ
度	m/秒	m/秒	秒/m	秒/m
0	30.7	30.9	0.033	0.032
5	30.7	30.8	0.033	0.032
10	30.7	30.5	0.033	0.033
15	30.7	29.9	0.033	0.033
20	30.8	29.1	0.032	0.034
25	30.8	28.0	0.032	0.036
30	30.8	26.8	0.032	0.037
35	30.8	25.4	0.032	0.039
40	30.8	23.7	0.032	0.042
45	30.9	21.9	0.032	0.046
50	30.9	19.9	0.032	0.050
55	30.9	17.8	0.032	0.056
60	30.9	15.5	0.032	0.065
65	31.0	13.1	0.032	0.076
70	31.0	10.6	0.032	0.094
75	31.0	8.0	0.032	0.125
80	31.0	5.4	0.032	0.185
85	31.0	2.7	0.032	0.370

ヒュベニの公式（2地点間の距離）

　2地点間の距離は、ヒュベニの公式を使って求めることができます。
　公式は、MathMLで記述しているため、MathJaxを使って表示しています。

ヒュベニの公式（2地点間の距離）	$dist = \sqrt{{(d_{y} M)}^{2} + {(d_{x} N cos a_{y})}^{2}}$
経度の差（ラジアン）	$d_{x} = x_{2} - x_{1}$
緯度の差（ラジアン）	$d_{y} = y_{2} - y_{1}$
緯度の平均値（ラジアン）	$a_{y} = \frac{y_{1} + y_{2}}{2}$
子午線曲率半径	$M = \frac{a (1 - e^{2})}{W^{3}} = \frac{6335439.327}{W^{3}}$
卯酉線曲率半径	$N = \frac{a}{W} = \frac{6378137.000}{W}$
子午線・卯酉線曲率半径の分母	$W = \sqrt{1 - e^{2} {sin}^{2} a_{y}}$
地球の長半径（m）	$a = 6378137.000$
地球の短半径（m）	$b = 6356752.314$
地球の扁平率	$f = \frac{a - b}{a} = \frac{1}{298.2572221}$
地球の離心率の2乗	$e^{2} = f (2 - f) = 0.006694380023$

度とラジアンの変換

度からラジアンへの変換	rad = deg * PI / 180 = deg * 0.01745329252
ラジアンから度への変換	deg = rad / PI * 180 = rad * 57.29577951
PI（円周率）	PI = 3.141592654