xybl.inc.php
下記ファイルの内容を表示しています。 ダウンロードを行うにはファイル名をクリックしてください。
<?php
//============================================================================
// 緯度経度⇔平面座標 変換
//============================================================================
// 原点座標位置を戻す
function gentenset($prj, &$b0,&$L0) {
switch ($prj) {
case 1:$b0=33.0; $L0=129.5; break;
case 2:$b0=33.0; $L0=131.0; break;
case 3:$b0=36.0; $L0=132.1666666666666; break;
case 4:$b0=33.0; $L0=133.5; break;
case 5:$b0=36.0; $L0=134.3333333333333; break;
case 6:$b0=36.0; $L0=136.0; break;
case 7:$b0=36.0; $L0=137.1666666666666; break;
case 8:$b0=36.0; $L0=138.5; break;
case 9:$b0=36.0; $L0=139.833333333333; break;
case 10:$b0=40.0; $L0=140.833333333333; break;
case 11:$b0=44.0; $L0=140.25; break;
case 12:$b0=44.0; $L0=142.25; break;
case 13:$b0=44.0; $L0=144.25; break;
case 14:$b0=26.0; $L0=142.0; break;
case 15:$b0=26.0; $L0=127.5; break;
case 16:$b0=26.0; $L0=124.0; break;
case 17:$b0=26.0; $L0=131.0; break;
case 18:$b0=20.0; $L0=136.0; break;
case 19:$b0=26.0; $L0=154.0; break;
}
}
// BLtoXY 緯度経度を平面座標（測地座標）に変換
//
// K :19座標系の系番号
// Lp:経度の度点
// bp:緯度の度点
// var rX:19座標系のX座標
// var rY:19座標系のY座標
// 戻り値=子午線収差(度点表記)
//
function BLtoXY($K, $Lp,$bp, &$rX,&$rY) {
gentenset($K,$b0p,$L0p);
$ra=6377397.155;
$rb=6356078.965;
$k1=1.00503730604855;
$k2=5.0478492403e-3;
$k3=1.0563786831e-5;
$k4=2.0633322e-8;
$m0=0.9999;
$b=$bp; $L=$Lp; $b0=$b0p; $L0=$L0p;
$r0=180/pi();
$b0=$b0/$r0;
$L0=$L0/$r0;
$b=$b/$r0;
$L=$L/$r0;
$e =sqrt(1-pow(($rb/$ra),2));
$e1=sqrt(pow(($ra/$rb),2)-1);
$w=sqrt(1-pow(($e*sin($b)),2));
$n=$ra/$w;
$a1=$k1*($b-$b0);
$a2=$k2*(sin(2*$b)-sin(2*$b0))/2;
$a3=$k3*(sin(4*$b)-sin(4*$b0))/4;
$a4=$k4*(sin(6*$b)-sin(6*$b0))/6;
$db=$ra*(1-pow($e,2))*($a1-$a2+$a3-$a4);
$t=tan($b);
$h=$e1*cos($b);
$dl=$L-$L0;
$x=$db+pow($dl,2)*$n/2*sin($b)*cos($b);
$x=$x+pow($dl,4)*$n/24*sin($b)*pow(cos($b),3)*(5-pow($t,2)+9*pow($h,2)+4*pow($h,4));
$x=$x+pow($dl,6)*$n/720*sin($b)*pow(cos($b),5)*(61-58*pow($t,2)+pow($t,4)+270*pow($h,2)-330*pow(($h*$t),2));
$y=$dl*$n*cos($b)+pow($dl*cos($b),3)*$n/6*(1-pow($t,2)+pow($h,2));
$y=$y+pow(($dl*cos($b)),5)*$n/120*(5-18*pow($t,2)+pow($t,4)+14*pow($h,2)-58*pow(($h*$t),2));
$g=$dl*sin($b)+pow($dl,3)/3*sin($b)*pow(cos($b),2)*(1+3*pow($h,2)+2*pow($h,4))+pow($dl,5)/15*sin($b)*pow(cos($b),4)*(2-pow($t,2));
$x=$x*$m0;
$y=$y*$m0;
$g=$g*$r0;
$rX=$y;
$rY=$x;
return($g);
}
// XYtoBL 測地座標を緯度経度に変換する
//
// K :19座標系の系番号
// yp:測地Ｘ座標
// xp:測地Ｙ座標
// var oB:経度
// var oL:緯度
// 戻り値=子午線収差
//
function XYtoBL($K, $yp,$xp, &$oB,&$oL) {
gentenset($K,$b0p,$L0p);
$ra=6377397.155;
$rb=6356078.963;
$k1=1.005037306048555;
$k2=5.0478492403e-3;
$k3=1.0563786831e-5;
$k4=2.0633322e-8;
$m0=0.9999;
$x=$xp; $y=$yp; $b0=$b0p; $L0=$L0p;
$r0=180/pi();
$b0=$b0/$r0;
$L0=$L0/$r0;
$e =sqrt(1-pow(($rb/$ra),2));
$e1=sqrt(pow(($ra/$rb),2)-1);
$a=$ra*(1-pow($e,2));
$d=$a*($k1*$b0 - $k2/2*sin(2*$b0) + $k3/4*sin(4*$b0) - $k4/6*sin(6*$b0)) + $x/$m0;
for ($i=1; $i<=6; $i++) {
if ($i==1) {
$p[$i]=$b0;
} else {
$p[$i]=$p[$i-1]-($s[$i-1]-$d)/$m[$i-1];
}
$s[$i]=$a*($k1*$p[$i] - $k2/2*sin(2*$p[$i]) + $k3/4*sin(4*$p[$i]) - $k4/6*sin(6*$p[$i]));
$m[$i]=$a/(pow((1-pow(($e*sin($p[$i])),2)),1.5));
}
$pp=$p[5];
$y1=$y/$m0;
$t =tan($pp);
$h =$e1*cos($pp);
$n =$ra/(pow((1-pow(($e*sin($pp)),2)),0.5));
$mm=$m[5];
$LB=$pp-(pow($y1,2)*$t/(2*$mm*$n)-pow($y1,4)*$t*(5+3*pow($t,2)
+pow($h,2)-9*pow(($t*$h),2))/(24*$mm*pow($n,3)));
$LB=$LB-(pow($y1,6)*$t*(61+90*pow($t,2)+45*pow($t,4))/(720*$mm*pow($n,5)));
$L=($y1/($n*cos($pp))-pow($y1,3)*(1+2*pow($t,2)+pow($h,2))/(6*pow($n,3)*cos($pp)));
$L=$L+ (pow($y1,5)*(5+28*pow($t,2)+24*pow($t,4))/(120*pow($n,5)*cos($pp)));
$L=$L+$L0;
$G=($y1*$t/$n-pow($y1,3)*$t*(1+pow($t,2)-pow($h,2))/(3*pow($n,3))
+pow(($y1/$n),5)*$t*(2+5*pow($t,2)+3*pow($t,4))/15);
$L=$L*$r0;
$LB=$LB*$r0;
$G=$G*$r0;
$oL=$LB;
$oB=$L;
return($G);
}
?>