计算两个经纬度间的距离(c++)

double D_jw(double wd1,double jd1,double wd2,double jd2)
 {
     double x,y,out;
     double PI=3.14159265;
     double R=6.371229*1e6;

    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;
 }

==

  一个经纬度相关计算的C 类
写了一个经纬度距离计算的类

 

--------------CJWD.h--------------

#ifndef __JWD_AND_HELPER_20051005
 #define __JWD_AND_HELPER_20051005

#include "stdafx.h"
 #include <math.h>
 #include <iostream>
 using namespace std;

#ifndef PI

#define PI 3.14159265;

#endif
 static double Rc = 6378137;  // 赤道半径

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;
 }