File "Pos2Loc.php"
Full path: /www/ansys/qgis/common/class/Pos2Loc.php
File size: 8.97 KiB (9186 bytes)
MIME-type: text/x-php
Charset: utf-8
<?php
//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
// 緯度経度→平面直角座標
// クラス名： cPos2Loc
// 備考 ：
// 更新 ： 2012/12/06 KMK-Teraoka
//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
class Pos2Loc
{
private static $mEp_a;
private static $mEp_f1;
private static $mEp_e;
private static $mAEE;
private static $mCEE;
private static $mEp2;
private static $mPlaneIndex = -1;
// 1～19系の原点北緯
private static $LAT_ZERO_LIST = array(
"0" => 33.000000,
"1" => 33.000000,
"2" => 36.000000,
"3" => 33.000000,
"4" => 36.000000,
"5" => 36.000000,
"6" => 36.000000,
"7" => 36.000000,
"8" => 35.553333,
"9" => 40.000000,
"10" => 44.000000,
"11" => 44.000000,
"12" => 44.000000,
"13" => 26.000000,
"14" => 26.000000,
"15" => 26.000000,
"16" => 26.000000,
"17" => 20.000000,
"18" => 26.000000
);
// 1～19系の原点東経
private static $LON_ZERO_LIST = array(
"0" => 129.500000,
"1" => 131.000000,
"2" => 132.166667,
"3" => 133.500000,
"4" => 134.333333,
"5" => 136.000000,
"6" => 137.166667,
"7" => 138.500000,
"8" => 139.781111,
"9" => 140.833333,
"10" => 140.250000,
"11" => 142.250000,
"12" => 144.250000,
"13" => 142.000000,
"14" => 127.500000,
"15" => 124.000000,
"16" => 131.000000,
"17" => 136.000000,
"18" => 154.000000
);
private static $mJ = array(
"0" => 0.0,
"1" => 0.0,
"2" => 0.0,
"3" => 0.0,
"4" => 0.0,
"5" => 0.0,
"6" => 0.0,
"7" => 0.0,
"8" => 0.0
);
//
// 座標変換を使用するための初期化
// 平面直角座標を求める際に使用する変数を計算する
//
public static function init($planeIndex_)
{
if (0 != strcmp($planeIndex_, ""))
{
self::$mPlaneIndex = $planeIndex_;
$e2 = 0.0;
$e4 = 0.0;
$e6 = 0.0;
$e8 = 0.0;
$e10 = 0.0;
$e12 = 0.0;
$e14 = 0.0;
$e16 = 0.0;
self::$mEp_a = 6378137.0;
self::$mEp_f1 = 298.257222101;
self::$mEp_e = (2.0 * self::$mEp_f1 - 1.0) / self::$mEp_f1 / self::$mEp_f1;
$e2 = self::$mEp_e;
$e4 = $e2 * $e2;
$e6 = $e4 * $e2;
$e8 = $e4 * $e4;
$e10 = $e8 * $e2;
$e12 = $e8 * $e4;
$e14 = $e8 * $e6;
$e16 = $e8 * $e8;
// 定数項
self::$mAEE = self::$mEp_a * (1 - self::$mEp_e); // a(1 - e2)
self::$mCEE = self::$mEp_a / sqrt(1 - self::$mEp_e); // C = a * sqrt(1 + e'2) = a / sqrt(1 - e2)
self::$mEp2 = self::$mEp_e / (1 - self::$mEp_e); // e'2 (e prime 2) Eta2phiを計算するため
// 「緯度を与えて赤道からの子午線弧長を求める計算」のための９つの係数を求める。
// 「精密測地網一次基準点測量計算式」P55,P56より。係数チェック済み1999/10/19。
self::$mJ["0"] = 4927697775.0 / 7516192768.0 * $e16;
self::$mJ["0"] = self::$mJ["0"] + 19324305.0 / 29360128.0 * $e14;
self::$mJ["0"] = self::$mJ["0"] + 693693.0 / 1048576.0 * $e12;
self::$mJ["0"] = self::$mJ["0"] + 43659.0 / 65536.0 * $e10;
self::$mJ["0"] = self::$mJ["0"] + 11025.0 / 16384.0 * $e8;
self::$mJ["0"] = self::$mJ["0"] + 175.0 / 256.0 * $e6;
self::$mJ["0"] = self::$mJ["0"] + 45.0 / 64.0 * $e4;
self::$mJ["0"] = self::$mJ["0"] + 3.0 / 4.0 * $e2;
self::$mJ["0"] = self::$mJ["0"] + 1.0;
self::$mJ["1"] = 547521975.0 / 469762048.0 * $e16;
self::$mJ["1"] = self::$mJ["1"] + 135270135.0 / 117440512.0 * $e14;
self::$mJ["1"] = self::$mJ["1"] + 297297.0 / 262144.0 * $e12;
self::$mJ["1"] = self::$mJ["1"] + 72765.0 / 65536.0 * $e10;
self::$mJ["1"] = self::$mJ["1"] + 2205.0 / 2048.0 * $e8;
self::$mJ["1"] = self::$mJ["1"] + 525.0 / 512.0 * $e6;
self::$mJ["1"] = self::$mJ["1"] + 15.0 / 16.0 * $e4;
self::$mJ["1"] = self::$mJ["1"] + 3.0 / 4.0 * $e2;
self::$mJ["2"] = 766530765.0 / 939524096.0 * $e16;
self::$mJ["2"] = self::$mJ["2"] + 45090045.0 / 58720256.0 * $e14;
self::$mJ["2"] = self::$mJ["2"] + 1486485.0 / 2097152.0 * $e12;
self::$mJ["2"] = self::$mJ["2"] + 10395.0 / 16384.0 * $e10;
self::$mJ["2"] = self::$mJ["2"] + 2205.0 / 4096.0 * $e8;
self::$mJ["2"] = self::$mJ["2"] + 105.0 / 256.0 * $e6;
self::$mJ["2"] = self::$mJ["2"] + 15.0 / 64.0 * $e4;
self::$mJ["3"] = 209053845.0 / 469762048.0 * $e16;
self::$mJ["3"] = self::$mJ["3"] + 45090045.0 / 117440512.0 * $e14;
self::$mJ["3"] = self::$mJ["3"] + 165165.0 / 524288.0 * $e12;
self::$mJ["3"] = self::$mJ["3"] + 31185.0 / 131072.0 * $e10;
self::$mJ["3"] = self::$mJ["3"] + 315.0 / 2048.0 * $e8;
self::$mJ["3"] = self::$mJ["3"] + 35.0 / 512.0 * $e6;
self::$mJ["4"] = 348423075.0 / 1879048192.0 * $e16;
self::$mJ["4"] = self::$mJ["4"] + 4099095.0 / 29360128.0 * $e14;
self::$mJ["4"] = self::$mJ["4"] + 99099.0 / 1048576.0 * $e12;
self::$mJ["4"] = self::$mJ["4"] + 3465.0 / 65536.0 * $e10;
self::$mJ["4"] = self::$mJ["4"] + 315.0 / 16384.0 * $e8;
self::$mJ["5"] = 26801775.0 / 469762048.0 * $e16;
self::$mJ["5"] = self::$mJ["5"] + 4099095.0 / 117440512.0 * $e14;
self::$mJ["5"] = self::$mJ["5"] + 9009.0 / 524288.0 * $e12;
self::$mJ["5"] = self::$mJ["5"] + 693.0 / 131072.0 * $e10;
self::$mJ["6"] = 11486475.0 / 939524096.0 * $e16;
self::$mJ["6"] = self::$mJ["6"] + 315315.0 / 58720256.0 * $e14;
self::$mJ["6"] = self::$mJ["6"] + 3003.0 / 2097152.0 * $e12;
self::$mJ["7"] = 765765.0 / 469762048.0 * $e16;
self::$mJ["7"] = self::$mJ["7"] + 45045.0 / 117440512.0 * $e14;
self::$mJ["8"] = 765765.0 / 7516192768.0 * $e16;
}
else
{
print "系番号が指定されていません";
}
}
//
// 変換処理
//
public static function convert($lat_, $lon_)
{
$ret = array(
"0" => 0.0,
"1" => 0.0
);
if (0 != $lat_ && 0 != $lon_ && -1 != self::$mPlaneIndex)
{
$DL = 0.0;
$DL2 = 0.0;
$DL4 = 0.0;
$DL6 = 0.0;
$s0 = 0.0; // 赤道から座標系原点までの子午線長の計算
$s = 0.0; // 赤道から求点までの子午線長の計算
$cos2 = 0.0;
$eta2phi = 0.0;
$nphi = 0.0;
$t = 0.0;
$t2 = 0.0;
$t4 = 0.0;
$t6 = 0.0;
$m0 = 1.0000;
$degTorad = 1.74532925199433E-02;
// 上記計算は定数に格納可能
$latRad = $lat_ * $degTorad;
$lonRad = $lon_ * $degTorad;
$lat0Rad = self::$LAT_ZERO_LIST[self::$mPlaneIndex - 1] * $degTorad;
$lon0Rad = self::$LON_ZERO_LIST[self::$mPlaneIndex - 1] * $degTorad;
$DL = $lonRad - $lon0Rad; // Δλ
//print $DL;
// 赤道から座標系原点までの子午線長の計算
$s0 = self::meridS($lat0Rad, self::$mAEE, self::$mJ["0"], self::$mJ["1"], self::$mJ["2"], self::$mJ["3"], self::$mJ["4"], self::$mJ["5"], self::$mJ["6"], self::$mJ["7"], self::$mJ["8"]);
// 何度も使う式を変数に代入
$t = tan($latRad);
$t2 = $t * $t;
$t4 = $t2 * $t2;
$t6 = $t4 * $t2;
$cos2 = cos($latRad) * cos($latRad);
$eta2phi = self::$mEp2 * $cos2; // =η1*η1
$nphi = self::$mCEE / sqrt(1 + $eta2phi); //「精密測地網一次基準点測量計算式」P52のN(phi)
$DL2 = $DL * $DL;
$DL4 = $DL2 * $DL2;
$DL6 = $DL4 * $DL2;
// 赤道から求点までの子午線長の計算
$s = self::meridS($latRad, self::$mAEE, self::$mJ["0"], self::$mJ["1"], self::$mJ["2"], self::$mJ["3"], self::$mJ["4"], self::$mJ["5"], self::$mJ["6"], self::$mJ["7"], self::$mJ["8"]);
// X，Yの計算 in meter
// 「精密測地網一次基準点測量計算式」P52,53のx,yを求める式より
$convertX = 0.0;
$convertX = -(-1385.0 + 3111.0 * $t2 - 543.0 * $t4 + $t6) * $DL6 * pow($cos2, 3.0) / 40320.0;
$convertX = $convertX - (-61.0 + 58.0 * $t2 - $t4 - 270.0 * $eta2phi + 330.0 * $t2 * $eta2phi) * $DL4 * pow($cos2, 2.0) / 720.0;
$convertX = $convertX + (5 - $t2 + 9 * $eta2phi + 4 * $eta2phi * $eta2phi) * $DL2 * $cos2 / 24;
$convertX = $convertX + 1.0 / 2.0;
$convertX = $convertX * $nphi * $cos2 * $t * $DL2;
$convertX = $convertX + $s - $s0;
$convertX = $convertX * $m0;
$convertY = 0.0;
$convertY = -(-61.0 + 479.0 * $t2 - 179.0 * $t4 + $t6) * $DL6 * pow($cos2, 3.0) / 5040.0;
$convertY = $convertY - (-5.0 + 18.0 * $t2 - $t4 - 14.0 * $eta2phi + 58.0 * $t2 * $eta2phi) * $DL4 * pow($cos2, 2.0) / 120.0;
$convertY = $convertY - (-1.0 + $t2 - $eta2phi) * $DL2 * $cos2 / 6.0;
$convertY = $convertY + 1.0;
$convertY = $convertY * $nphi * cos($latRad) * $DL;
$convertY = $convertY * $m0;
$ret["0"] = $convertX;
$ret["1"] = $convertY;
}
else
{
print "緯度経度が指定されていません";
}
return $ret;
}
//
// 子午線弧長を計算
//
private static function meridS($phi_, $aee_, $aJ_, $bJ_, $cJ_, $dJ_, $eJ_, $fJ_, $gJ_, $hJ_, $iJ_)
{
$ss = 0.0;
$ss = $iJ_ / 16.0 * sin(16.0 * $phi_);
$ss = $ss - $hJ_ / 14.0 * sin(14.0 * $phi_);
$ss = $ss + $gJ_ / 12.0 * sin(12.0 * $phi_);
$ss = $ss - $fJ_ / 10.0 * sin(10.0 * $phi_);
$ss = $ss + $eJ_ / 8.0 * sin(8.0 * $phi_);
$ss = $ss - $dJ_ / 6.0 * sin(6.0 * $phi_);
$ss = $ss + $cJ_ / 4.0 * sin(4.0 * $phi_);
$ss = $ss - $bJ_ / 2.0 * sin(2.0 * $phi_);
$ss = $ss + $aJ_ * $phi_;
$ss = $aee_ * $ss;
return $ss;
}
}
?>