根據兩站點的經緯度求兩站點間的距離
/**** 根據兩站點的經緯度求兩站點間的距離 ****/
double D_jw(double wd1,double jd1,double wd2,double jd2)
{
double x,y,out;
double PI=3.14159265;
double R=6.371229*1e6;ios
x=(jd2-jd1)*PI*R*cos( ((wd1 wd2)/2) *PI/180)/180;
y=(wd2-wd1)*PI*R/180;
out=hypot(x,y);
return out/1000;
}git
==算法
一個經緯度相關計算的C 類
寫了一個經緯度距離計算的類函數
--------------CJWD.h--------------測試
#ifndef __JWD_AND_HELPER_20051005
#define __JWD_AND_HELPER_20051005google
#include "stdafx.h"
#include <math.h>
#include <iostream>
using namespace std;spa
#ifndef PI.net
#define PI 3.14159265;blog
#endif
static double Rc = 6378137; // 赤道半徑get
static double Rj = 6356725; // 極半徑
namespace CDYW{
class JWD
{
public:
double m_LoDeg, m_LoMin, m_LoSec; // longtitude 經度
double m_LaDeg, m_LaMin, m_LaSec;
double m_Longitude, m_Latitude;
double m_RadLo, m_RadLa;
double Ec;
double Ed;
public:
// 構造函數, 經度: loDeg 度, loMin 分, loSec 秒; 緯度: laDeg 度, laMin 分, laSec秒
JWD(double loDeg, double loMin, double loSec, double laDeg, double laMin, double laSec)
{
m_LoDeg=loDeg; m_LoMin=loMin; m_LoSec=loSec; m_LaDeg=laDeg; m_LaMin=laMin; m_LaSec=laSec;
m_Longitude = m_LoDeg m_LoMin / 60 m_LoSec / 3600;
m_Latitude = m_LaDeg m_LaMin / 60 m_LaSec / 3600;
m_RadLo = m_Longitude * PI / 180.;
m_RadLa = m_Latitude * PI / 180.;
Ec = Rj (Rc - Rj) * (90.- m_Latitude) / 90.;
Ed = Ec * cos(m_RadLa);
}
//!
JWD(double longitude, double latitude)
{
m_LoDeg = int(longitude);
m_LoMin = int((longitude - m_LoDeg)*60);
m_LoSec = (longitude - m_LoDeg - m_LoMin/60.)*3600;
m_LaDeg = int(latitude);
m_LaMin = int((latitude - m_LaDeg)*60);
m_LaSec = (latitude - m_LaDeg - m_LaMin/60.)*3600;
m_Longitude = longitude;
m_Latitude = latitude;
m_RadLo = longitude * PI/180.;
m_RadLa = latitude * PI/180.;
Ec = Rj (Rc - Rj) * (90.-m_Latitude) / 90.;
Ed = Ec * cos(m_RadLa);
}
};
class CJWDHelper
{
public:
CJWDHelper() {};
~CJWDHelper() {};
//! 計算點A 和 點B的經緯度,求他們的距離和點B相對於點A的方位
/*!
* \param A A點經緯度
* \param B B點經緯度
* \param angle B相對於A的方位, 不須要返回該值,則將其設爲空
* \return A點B點的距離
*/
static double distance(JWD A, JWD B, double *angle)
{
double dx = (B.m_RadLo - A.m_RadLo) * A.Ed;
double dy = (B.m_RadLa - A.m_RadLa) * A.Ec;
double out = sqrt(dx * dx dy * dy);
if( angle != NULL)
{
*angle = atan(fabs(dx/dy))*180./PI;
// 判斷象限
double dLo = B.m_Longitude - A.m_Longitude;
double dLa = B.m_Latitude - A.m_Latitude;
if(dLo > 0 && dLa <= 0) {
*angle = (90. - *angle) 90.;
}
else if(dLo <= 0 && dLa < 0) {
*angle = *angle 180.;
}
else if(dLo < 0 && dLa >= 0) {
*angle = (90. - *angle) 270;
}
}
return out/1000;
}
//! 計算點A 和 點B的經緯度,求他們的距離和點B相對於點A的方位
/*!
* \param longitude1 A點經度
* \param latitude1 A點緯度
* \param longitude2 B點經度
* \param latitude2 B點緯度
* \param angle B相對於A的方位, 不須要返回該值,則將其設爲空
* \return A點B點的距離
*/
static double distance(
double longitude1, double latitude1,
double longitude2, double latitude2,
double *angle)
{
JWD A(longitude1,latitude1);
JWD B(longitude2,latitude2);
return distance(A, B, angle);
}
//! 已知點A經緯度,根據B點據A點的距離,和方位,求B點的經緯度
/*!
* \param A 已知點A
* \param distance B點到A點的距離
* \param angle B點相對於A點的方位
* \return B點的經緯度座標
*/
static JWD GetJWDB(JWD A, double distance, double angle)
{
double dx = distance*1000 * sin(angle * PI /180.);
double dy = distance*1000 * cos(angle * PI /180.);
//double dx = (B.m_RadLo - A.m_RadLo) * A.Ed;
//double dy = (B.m_RadLa - A.m_RadLa) * A.Ec;
double BJD = (dx/A.Ed A.m_RadLo) * 180./PI;
double BWD = (dy/A.Ec A.m_RadLa) * 180./PI;
JWD B(BJD, BWD);
return B;
}
//! 已知點A經緯度,根據B點據A點的距離,和方位,求B點的經緯度
/*!
* \param longitude 已知點A經度
* \param latitude 已知點A緯度
* \param distance B點到A點的距離
* \param angle B點相對於A點的方位
* \return B點的經緯度座標
*/
static JWD GetJWDB(double longitude, double latitude, double distance, double angle)
{
JWD A(longitude,latitude);
return GetJWDB(A, distance, angle);
}
};
}
#endif
=========== 測試程序==========
#include "stdafx.h"
#include <math.h>
#include <iostream>#include "CJWD.h"
using namespace std;using namespace CDYW;
double Rc = 6378137; // 赤道半徑
double Rj = 6356725; // 極半徑// 綿陽
double jd1 = 104.740999999;
double wd1 = 31.4337;// 成都
double jd2 = 104.01;
double wd2 = 30.40; int main(int argc, char* argv[])
{
double angle = 0;
cout << "A(綿陽): JD = " << jd1 << " WD = " << wd1 << endl;
cout << "B(成都): JD = " << jd2 << " WD = " << wd2 << endl;
cout << "--------------------" << endl;
cout << D_jw(wd1,jd1,wd2,jd2, angle) << endl;
cout << "angle: " << angle <<endl;
cout << "==============" <<endl;
JWD A(jd1,wd1),B(jd2,wd2);
double distance = CJWDHelper::distance(jd1,wd1,jd2,wd2, &angle);
//cout << CJWDHelper::distance(A,B, &angle) << endl;
cout << distance << endl;
cout << "angle: " << angle <<endl;
cout << "==============" <<endl;
JWD C = CJWDHelper::GetJWDB(A, distance, angle);
cout << "JD = " << C.m_Longitude << " WD = " << C.m_Latitude << endl;
cout << "==============" <<endl;
cout << A.m_LoDeg << " " << A.m_LoMin << " " << A.m_LoSec << endl; return 0;
}
=====
經過兩個點的經緯度計算距離
關鍵詞: gis
從google maps的腳本里扒了段代碼,沒準啥時會用上。你們一塊看看是怎麼算的。
private const double EARTH_RADIUS = 6378.137;
private static double rad(double d)
{
return d * Math.PI / 180.0;
}
public static double GetDistance(double lat1, double lng1, double lat2, double lng2)
{
double radLat1 = rad(lat1);
double radLat2 = rad(lat2);
double a = radLat1 - radLat2;
double b = rad(lng1) - rad(lng2);
double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a/2),2)
Math.Cos(radLat1)*Math.Cos(radLat2)*Math.Pow(Math.Sin(b/2),2)));
s = s * EARTH_RADIUS;
s = Math.Round(s * 10000) / 10000;
return s;
}
很是感謝,幫了我大忙了:)
雖然我也沒看明白到底原理是什麼,但驗算了A(60,30),B(60,90)兩點之間,此段代碼和我用餘弦定理算出來的結果很一致。
餘弦定理的步驟是:一、算A、B弦長:地球半徑R*cos(經度差60)=R/2;
二、算角AOB,O爲地球圓心,利用餘弦定理,
cosAOB=(2R*R-(R/2)^2) /2*R*R=7/8;
三、弧AB的長爲:R*arc cos(7/8);求畢 。