日本測地系から世界測地系に変換する

超お久しぶりのポスト。

じゃらんWebサービス - 空室検索API

じゃらんWebサービスは、じゃらんnetの提供する宿泊施設情報を活用して、開発者やサイト制作者、アフィリエイトサイトオーナーの皆様のために情報提供をしています。

www.jalan.net

最近じゃらんWebサービスを使って遊んだのだけれど、APIから戻ってくる緯度・経度が日本測地系でそのまんま地図にプロットするとズレッズレやないか！　で困っちゃった。

1 困った時の他人頼み
2 困った時の諦めない人頼み
3 けっこうガチでPHPへ再移植したった
4 ねえ君、これが誤差の範囲だと思うんだ

困った時の他人頼み

どうするかなんてググれば誰かが教えてくれるさ！
さあこい！　変換式こい！

じゃらんWebサービス - チュートリアル

www.jalan.net

というか、ここに書いてあるじゃん！
ググるといっぱい出てくるの 0.00010695 とかをかけたり足したりする簡易計算なるモノで変換してみたんです。

で、プロットしてみた。
いやいや、これ・・・　ズレッズレやないかあああああああああああ！

数百ｍ単位でずれてるずれてる。
まあ誤差の範囲ですね！で締めくくっちゃってる人。
範囲外だと思うよ！　諦めないで！！　最後まで走り抜けて！！！

困った時の諦めない人頼み

いやいや、きっといるはずだよ。
ズレッズレじゃない方法を諦めずに確立して提供してくれる人が！

で、たどり着きました。

日本測地系から世界測地系に変換する（けっこうガチで計算するVer） – 東京伊勢海老通信

マジ神や！　ありがとうきたむらさん！
何々？ PHPで動くのがあって、それをRubyに移植したんですね。
そのPHPのやつください・・・。
ありませんか、そうですか・・・。

けっこうガチでPHPへ再移植したった

もっと探せばその元になったブツにたどり着くかもだけど、0.010041046 とかを引いたりかけたりするやつばっかりの検索結果に嫌気が指したので、素直にPHPへ再移植。
たぶんまともに動くはずだよ。
何をどう計算してるかなんて正直全然わかんないよ。

class jgd2wgs
{
	private $RD;
	private $E2_JP;
	private $E2_WS;
	private $R_JP = 6377397.155;
	private $F_JP = 1 / 299.1528128;
	private $R_WS = 6378137.0;
	private $F_WS = 1 / 298.257223563;
	private $DX = - 148.0;
	private $DY = 507.0;
	private $DZ = 681.0;
	private $HEIGHT = 0;

	public function __construct()
	{
		$this->RD = pi() / 180;
		$this->E2_JP = 2 * $this->F_JP - $this->F_JP * $this->F_JP;
		$this->E2_WS = 2 * $this->F_WS - $this->F_WS * $this->F_WS;
	}

	private function llh2xyz($lat, $lng, $h, $a, $e2, $rd)
	{
		$lat *= $rd;
		$lng *= $rd;
		$sb = sin($lat);
		$cb = cos($lat);
		$rn = $a / sqrt(1 - $e2 * $sb * $sb);
		$x = ($rn + $h) * $cb * cos($lng);
		$y = ($rn + $h) * $cb * sin($lng);
		$z = ($rn * (1 - $e2) + $h) * $sb;
		return array(
			$x,
			$y,
			$z
		);
	}

	private function xyz2llh($x, $y, $z, $a, $e2, $rd)
	{
		$bda = sqrt(1 - $e2);
		$p = sqrt($x * $x + $y * $y);
		$t = atan2($z, $p * $bda);
		$st = sin($t);
		$ct = cos($t);
		$b = atan2($z + $e2 * $a / $bda * $st * $st * $st, $p - $e2 * $a * $ct * $ct * $ct);
		$l = atan2($y, $x);
		$sb = sin($b);
		$rn = $a / sqrt(1 - $e2 * $sb * $sb);
		$h = $p / cos($b) - $rn;
		$l1 = $b / $rd;
		$l2 = $l / $rd;
		return array(
			$l1,
			$l2,
			$h
		);
	}

	public function convert(&$lat, &$lng)
	{
		list ($x, $y, $z) = $this->llh2xyz($lat, $lng, $this->HEIGHT, $this->R_JP, $this->E2_JP, $this->RD);
		$x += $this->DX;
		$y += $this->DY;
		$z += $this->DZ;
		list ($lat, $lng, $h) = $this->xyz2llh($x, $y, $z, $this->R_WS, $this->E2_WS, $this->RD);
	}
}

// じゃらんのミリ秒緯度経度を
$lat = 127830116;
$lng = 496212164;
// 度の緯度経度に変換して
$lat /= 3600000;
$lng /= 3600000;
// そっから世界測地系へ変換しちゃうぜ！
$jgd2wgs = new jgd2wgs();
$jgd2wgs->convert($lat, $lng);

echo $lat . ',' . $lng;

class jgd2wgs

{

private $RD;

private $E2_JP;

private $E2_WS;

private $R_JP = 6377397.155;

private $F_JP = 1 / 299.1528128;

private $R_WS = 6378137.0;

private $F_WS = 1 / 298.257223563;

private $DX = - 148.0;

private $DY = 507.0;

private $DZ = 681.0;

private $HEIGHT = 0;

public function __construct()

{

$this->RD = pi() / 180;

$this->E2_JP = 2 * $this->F_JP - $this->F_JP * $this->F_JP;

$this->E2_WS = 2 * $this->F_WS - $this->F_WS * $this->F_WS;

}

private function llh2xyz($lat, $lng, $h, $a, $e2, $rd)

{

$lat *= $rd;

$lng *= $rd;

$sb = sin($lat);

$cb = cos($lat);

$rn = $a / sqrt(1 - $e2 * $sb * $sb);

$x = ($rn + $h) * $cb * cos($lng);

$y = ($rn + $h) * $cb * sin($lng);

$z = ($rn * (1 - $e2) + $h) * $sb;

return array(

$x,

$y,

);

}

private function xyz2llh($x, $y, $z, $a, $e2, $rd)

{

$bda = sqrt(1 - $e2);

$p = sqrt($x * $x + $y * $y);

$t = atan2($z, $p * $bda);

$st = sin($t);

$ct = cos($t);

$b = atan2($z + $e2 * $a / $bda * $st * $st * $st, $p - $e2 * $a * $ct * $ct * $ct);

$l = atan2($y, $x);

$sb = sin($b);

$rn = $a / sqrt(1 - $e2 * $sb * $sb);

$h = $p / cos($b) - $rn;

$l1 = $b / $rd;

$l2 = $l / $rd;

return array(

$l1,

$l2,

);

}

public function convert(&$lat, &$lng)

{

list ($x, $y, $z) = $this->llh2xyz($lat, $lng, $this->HEIGHT, $this->R_JP, $this->E2_JP, $this->RD);

$x += $this->DX;

$y += $this->DY;

$z += $this->DZ;

list ($lat, $lng, $h) = $this->xyz2llh($x, $y, $z, $this->R_WS, $this->E2_WS, $this->RD);

}

// じゃらんのミリ秒緯度経度を

$lat = 127830116;

$lng = 496212164;

// 度の緯度経度に変換して

$lat /= 3600000;

$lng /= 3600000;

// そっから世界測地系へ変換しちゃうぜ！

$jgd2wgs = new jgd2wgs();

$jgd2wgs->convert($lat, $lng);

echo $lat . ',' . $lng;

ねえ君、これが誤差の範囲だと思うんだ

計算された、35.511572171709,137.83366470446 で見ると

Google Maps

Find local businesses, view maps and get driving directions in Google Maps.

www.google.com

なかなかいいトコついてるついてる。