經緯度距離計算方法

因為一些專案的需求,需要計算兩個 GPS 的點之間的距離。
因此拜 G神 的偉大能力,讓我找到了 經緯度計算距離公式 這一篇文章。

大致上,由兩點的經緯度去算換其距離,就是球面座標的距離計算。
其實這東西以前大學時代就有學過了,不過公式我是已經忘的一乾二淨了。因此還是要在找尋一下資料還解決問題啦。
不過,簡單的正圓球面座標系統的計算法,對於地球來說是有一點點的不正確,因為地球其實是一個橢圓形的物體。不過網路就這麼神奇與偉大的,已經有人 Thaddeus Vincenty 研究出來橢圓球面的計算法 Vincenty’s formulae

當然,皇天不負苦心人,讓我找到了 線上運算 的網站。
而線上計算的網站,使用的是 Javascript 的方式進行計算。但我是專案內需要使用到,因此需要轉換成 php 的程式碼。
不過還好,Javascript 跟 php 的程式寫法沒有差多少,最主要就是花時間替變數前面加上 $ 字號可以運算。

分別有兩總計算
已知兩點的經緯度,計算兩點之間的距離 計算網頁

function distVincenty($lat1, $long1, $lat2, $long2) {
	$a = 6378137;
	$b = 6356752.314245;
	$f = 1/298.257223563;
	$L = deg2rad($long2 - $long1);
	$U1 = atan((1 - $f) * tan(deg2rad($lat1)));
	$U2 = atan((1 - $f) * tan(deg2rad($lat2)));
	$sinU1 = sin($U1);
	$cosU1 = cos($U1);
	$sinU2 = sin($U2);
	$cosU2 = cos($U2);
	$lambda = $L;
	$iterLimit = 100;
	do {
		$sinLambda = sin($lambda);
		$cosLambda = cos($lambda);
		$sinSigma = sqrt(($cosU2 * $sinLambda) * ($cosU2 * $sinLambda) + ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));
		if( $sinSigma == 0 ) {
			return 0;
		}
		$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
		$sigma = atan2($sinSigma, $cosSigma);
		$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
		$cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
		if( $cosSqAlpha == 0 ) {
			$cos2SigmaM = 0;
		} else {
			$cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
		}
		$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
		$lambdaP = $lambda;
		$lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
	} while(abs($lambda - $lambdaP) > 0.000000000001 && --$iterLimit > 0);
	if( $iterLimit == 0 ) {
		return 0;
	}
	$uSq = $cosSqAlpha * ($a * $a - $b * $b) / ($b * $b);
	$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
	$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
	$deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));
	return $b * $A * ($sigma - $deltaSigma);
}

已知 點1 的經緯度和到 點2 的角度與距離,求 點2 的經緯度 計算網頁

function place_len_to($lat, $long, $brng, $dist) {
	$a = 6378137;
	$b = 6356752.3142;
	$f = 1/298.257223563;
	$alpha1 = deg2rad($brng);
	$sinAlpha1 = sin($alpha1);
	$cosAlpha1 = cos($alpha1);
	$tanU1 = (1 - $f) * tan(deg2rad($lat));
	$cosU1 = 1 / sqrt((1 + $tanU1*$tanU1));
	$sinU1 = $tanU1 * $cosU1;
	$sigma1 = atan2($tanU1, $cosAlpha1);
	$sinAlpha = $cosU1 * $sinAlpha1;
	$cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
	$uSq = $cosSqAlpha * ($a * $a - $b * $b) / ($b * $b);
	$A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
	$B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
	$sigma = $dist / ($b * $A);
	$pi = pi();
	$sigmaP = 2 * $pi;
	while( abs($sigma - $sigmaP) > 0.000000000001 ) {
		$cos2SigmaM = cos(2 * $sigma1 + $sigma);
		$sinSigma = sin($sigma);
		$cosSigma = cos($sigma);
		$deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));
		$sigmaP = $sigma;
		$sigma = $dist / ($b * $A) + $deltaSigma;
	}
	$tmp = $sinU1 * $sinSigma - $cosU1 * $cosSigma * $cosAlpha1;
	$lat_end = atan2($sinU1 * $cosSigma + $cosU1 * $sinSigma * $cosAlpha1, (1 - $f) * sqrt($sinAlpha * $sinAlpha + $tmp * $tmp));
	$lambda = atan2($sinSigma * $sinAlpha1, $cosU1 * $cosSigma - $sinU1 * $sinSigma * $cosAlpha1);
	$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
	$L = $lambda - (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
	$long_end = fmod((deg2rad($long) + $L + 3 * $pi), 2 * $pi) - $pi;
	return array(rad2deg($lat_end), rad2deg($long_end));
}

需要注意的是,一般經緯度是以角度來表現的,但是 php 的三角函式是需要傳入弳度的,所以要記得轉換才可以正確的顯示與運算。
在第二個運算過程中,最後一個步驟求 點2 的經度,php 的 % 求餘算運算符號,計算出來的結果為 整數,但本運算中需要使用到符點數才會有正確的結果,所以需使用 fmod 才可以計算出正確的答案。

我也正在找經緯度距離計算方法的PHP Script, 謝謝您的代碼!! 很好用