緯度経度を平面直角座標に変換する
位置情報を扱っていると緯度経度の情報を平面座標系に変換したくなる時があるはず
やっぱり慣れ親しんだ平面直角座標系で計算したくなる時あるよね
saji.icon の場合位置情報をもとにコンテンツを表示するARとかやっていて欲しくなった
誤差の問題
なのであんまり広い範囲を考えると地球の丸みが無視できなくなる
具体的には緯度経度を平面直角座標にする際に元々球面だったものを平面に投影するので、基準点から離れるほど平面距離が増大してしまう。 この問題に関しては次に書く国土地理院の平面直角座標系が基準となると思う
国土地理院では以下のような座標系を採用している
誤差 : 投影距離の誤差を相対的に1/10000以内に収める
そのために座標原点を通る子午線上の縮尺係数を0.9999に設定
適用できる範囲 : 座標原点から東西約130km以内
これ以上だと上記の誤差に収まらなくなる
イメージ図
https://gyazo.com/49ec5bc66852955d4e42835f318afb14
上記の通り一つの座標系でカバーできるのは130km以内なので下図のように日本では19の座標系を設定している。
https://gyazo.com/ee4fd8304f42fd2622612099f5d3993d
→ つまるところ(気にする精度によるが)せいぜい原点から100Kmくらいしか平面直角座標は適応できないと考えた方がいいのかも(私見)
計算の前提と計算式
これを地球楕円体という
いくつかあるけど現在日本では測地基準系1980(GRS80)楕円体を採用している
長半径:6,378,137
扁平率の逆数:298.257222101
国土地理院の公開している計算式はこれ
なのだが
その量にちょっとびっくりすると思う
数学苦手な人は逃げ出すかもしれない
そのために本家国土地理院の計算方式を解説してくれている以下記事がが非常にわかりやすい
またPythonでのコードも書いてくれている
一応解説
convertLatLngToXY が本体
第1引数のlatLngが座標を求めたい緯度経度
第2引数のoriginLatLngが座標系の原点
もちろん北上の座標系
code:convertLatLngToXY.ts
interface LatLng {
latitude: number;
longitude: number;
}
// 定数
const m0 = 0.9999;
const a = 6378137;
const F = 298.257222101;
const n = 1 / (2 * F - 1);
const ALPHA_ARRAY = [
(1 / 2) * n -
(2 / 3) * n ** 2 +
(5 / 16) * n ** 3 +
(41 / 180) * n ** 4 -
(127 / 288) * n ** 5,
(13 / 48) * n ** 2 -
(3 / 5) * n ** 3 +
(557 / 1440) * n ** 4 +
(281 / 630) * n ** 5,
(61 / 240) * n ** 3 - (103 / 140) * n ** 4 + (15061 / 26880) * n ** 5,
(49561 / 161280) * n ** 4 - (179 / 168) * n ** 5,
(34729 / 80640) * n ** 5,
];
const A_ARRAY = [
-(3 / 2) * (n - n ** 3 / 8 - n ** 5 / 64),
(15 / 16) * (n ** 2 - n ** 4 / 4),
-(35 / 48) * (n ** 3 - (5 / 16) * n ** 5),
(315 / 512) * n ** 4,
-(693 / 1280) * n ** 5,
];
const A0 = 1 + n ** 2 / 4 + n ** 4 / 64;
const A_ = ((m0 * a) / (1 + n)) * A0;
export const convertLatLngToXY = (latLng: LatLng, originLatLng: LatLng) => {
const latitudeRadian = (Math.PI * latLng.latitude) / 180; //φ
const longitudeRadian = (Math.PI * latLng.longitude) / 180; //λ
const originLatitudeRadian = (Math.PI * originLatLng.latitude) / 180;
const originLongitudeRadian = (Math.PI * originLatLng.longitude) / 180;
const S_ =
((m0 * a) / (1 + n)) *
(A0 * originLatitudeRadian +
A_ARRAY.reduce((acc, cur, index) => {
return acc + cur * Math.sin(2 * (index + 1) * originLatitudeRadian);
}, 0));
const lambda_c = Math.cos(longitudeRadian - originLongitudeRadian);
const lambda_s = Math.sin(longitudeRadian - originLongitudeRadian);
const t = Math.sinh(
Math.atanh(Math.sin(latitudeRadian)) -
((2 * Math.sqrt(n)) / (1 + n)) *
Math.atanh(((2 * Math.sqrt(n)) / (1 + n)) * Math.sin(latitudeRadian)),
);
const t_ = Math.sqrt(1 + t * t);
const xi = Math.atan(t / lambda_c);
const eta = Math.atanh(lambda_s / t_);
const y =
A_ *
(xi +
ALPHA_ARRAY.reduce((acc, cur, index) => {
return (
acc +
cur *
Math.sin(2 * (index + 1) * xi) *
Math.cosh(2 * (index + 1) * eta)
);
}, 0)) -
S_;
const x =
A_ *
(eta +
ALPHA_ARRAY.reduce((acc, cur, index) => {
return (
acc +
cur *
Math.cos(2 * (index + 1) * xi) *
Math.sinh(2 * (index + 1) * eta)
);
}, 0));
return {x, y};
};
参照
国土地理院の正式な解説
sw1227さんによるわかりやすい計算の仕方の解説
同じく sw1227さんがPythonで書いてくださったコード