地球椭球体基本要素的计算

地球椭球体基本要素的计算，主要包括纬线弧长、子午线弧长、椭球面上梯形面积，以及同一个椭球体下大地坐标和空间直角坐标之间的转换等。为了方便，写了一个类如下，方便调用，在此也分享给大家：

头文件如下：

/*******************************************************************************
* 版权所有(C) 福建省空间信息工程研究中心 2013
* 文件名称  : GeoEllipse.h
* 当前版本  : 1.0.0.1
* 作    者    : 周旭光 (zxg)
* 设计日期  : 2013年9月22日
* 内容摘要  : 地球椭球体计算，主要有一些参数的计算和转换关系
* 修改记录  :
* 日    期        版    本      修改人     修改摘要********************************************************************************/
#ifndef __GEOELLIPSE_H_05958735_A7BF_409C_A444_CF9F202FCBE9__
#define __GEOELLIPSE_H_05958735_A7BF_409C_A444_CF9F202FCBE9__/**********************************  头文件 ************************************/#include <math.h>
#include <assert.h>
#include "CoordCommon.h"class COORDTRAN_API CGeoEllipse
{
public:explicit CGeoEllipse(double dbRadiusA,double dbRadiusB);~CGeoEllipse(void);double GetRadiusA() const{return m_dbRadiusA;}double GetRadiusB() const{return m_dbRadiusB;}double GetE1() const{return m_dbE1;}double GetE2() const{return m_dbE2;}//子午圈曲率半径,Mdouble GetMeridianRadius(double lat) const;//卯酉圈曲率半径,Ndouble GetPrimeRadius(double lat) const;//纬线圈的半径double GetLatCycleRadius(double lat) const;//平均曲率半径double GetAveRadius() const;//纬线弧长，即在指定纬度dbLat上经度差为dbDLon的纬线弧长double GetParallelArcLen(double dbLat,double dbDLon) const;//子午线弧长double GetMeridianArcLen(double dbLat) const;//椭球体梯形面积,即两条纬线和两条经线之间所夹的梯形面积double GetTrapeziumArea(double dbLat1,double dbLat2,double dbLon1,double dbLon2) const;//计算U值，dbLat以弧度为单位，等量纬度double GetValueU(double dbLat);//大地坐标转换为空间直角坐标bool BLH_XYZ(double dbLon,double dbLat,double dbHei,double& x,double &y,double &z);//空间直角坐标转换为大地坐标bool XYZ_BLH(double X,double Y,double Z, double &dbLon,double &dbLat,double &dbHei);private:double m_dbRadiusA;          //长半轴double m_dbRadiusB;            //短半轴double m_dbE1;             //第一偏心率double m_dbE2;               //第二偏心率
};#endif // end of __GEOELLIPSE_H_05958735_A7BF_409C_A444_CF9F202FCBE9__

实现文件如下：

#include "GeoEllipse.h"CGeoEllipse::CGeoEllipse(double dbRadiusA,double dbRadiusB)
:m_dbRadiusA(dbRadiusA),m_dbRadiusB(dbRadiusB)
{//m_dbRadiusA = dbRadiusA;//m_dbRadiusB = dbRadiusB;//计算第一和第二偏心率m_dbE1 = sqrt(m_dbRadiusA*m_dbRadiusA - m_dbRadiusB*m_dbRadiusB) / m_dbRadiusA;m_dbE2 = sqrt(m_dbRadiusA*m_dbRadiusA - m_dbRadiusB*m_dbRadiusB) / m_dbRadiusB;
}CGeoEllipse::~CGeoEllipse(void)
{
}double CGeoEllipse::GetMeridianRadius(double lat) const
{//先将纬度化作弧度制lat = lat *  acos(-1.0)/180.0;double dbNumerator = m_dbRadiusA*(1.0-m_dbE1*m_dbE1); double dbDenominator = pow(1.0 - m_dbE1*m_dbE1 * sin(lat)*sin(lat),1.5);return dbNumerator/dbDenominator;
}double CGeoEllipse::GetPrimeRadius(double lat) const
{assert(lat >= -90 && lat <= 90);//先将纬度化作弧度制lat = lat * acos(-1.0)/180.0;double dbDenominator = sqrt(1.0-m_dbE1*m_dbE1 * sin(lat)*sin(lat));return m_dbRadiusA/dbDenominator;
}double CGeoEllipse::GetLatCycleRadius(double lat) const
{assert(lat >= -90 && lat <= 90);double dbN = GetPrimeRadius(lat);lat = lat * acos(-1.0)/180.0;return dbN * cos(lat);
}double CGeoEllipse::GetParallelArcLen(double dbLat,double dbDLon) const
{//检查参数assert(dbLat >= -90 && dbLat <= 90&& dbDLon <= 180 && dbDLon >= 0);double dbN = GetPrimeRadius(dbLat);dbLat = dbLat * acos(-1.0)/180.0;dbDLon = dbDLon * acos(-1.0)/180.0;return dbN*cos(dbLat)*dbDLon;
}double CGeoEllipse::GetMeridianArcLen(double dbLat) const
{//检测参数是否合法assert(dbLat >= -90 && dbLat <= 90);//首先计算各个系数的值double dbA1 = 1.0 + 3.0/4.0*pow(m_dbE1,2.0) + 45.0/64.0*pow(m_dbE1,4.0)+175.0/256.0*pow(m_dbE1,6.0) + 11025.0/16384.0*pow(m_dbE1,8.0)+ 43659.0/65536.0*pow(m_dbE1,10.0) + 693693/1048576.0*pow(m_dbE1,12.0);double dbB1 = 3.0/4.0*pow(m_dbE1,2.0) + 45.0/64.0*pow(m_dbE1,4.0)+175.0/256.0*pow(m_dbE1,6.0) + 11025.0/16384.0*pow(m_dbE1,8.0)+ 43659.0/65536.0*pow(m_dbE1,10.0) + 693693/1048576.0*pow(m_dbE1,12.0);double dbC1 = 15.0/32.0*pow(m_dbE1,4.0) + 175.0/384.0*pow(m_dbE1,6.0) + 3675.0/8192.0*pow(m_dbE1,8.0) + 14553.0/32768.0*pow(m_dbE1,10.0) + 231231.0/524288.0*pow(m_dbE1,12.0);double dbD1 = 35.0/96.0*pow(m_dbE1,6.0) + 735.0/2048.0*pow(m_dbE1,8.0) + 14553.0/40960.0*pow(m_dbE1,10.0) + 231231.0/655360.0*pow(m_dbE1,12.0);double dbE1 = 315.0/1024.0*pow(m_dbE1,8.0) + 6237.0/20480.0*pow(m_dbE1,10.0) + 99099.0/327680.0*pow(m_dbE1,12.0);double dbF1 = 693.0/2560.0*pow(m_dbE1,10.0) + 11011.0/40960.0*pow(m_dbE1,12.0);double dbG1 = 1001.0/4096.0*pow(m_dbE1,12.0);//将角度转换为弧度dbLat = fabs(dbLat * acos(-1.0)/180.0);return m_dbRadiusA*(1.0-pow(m_dbE1,2.0)) * (dbA1*dbLat - cos(dbLat) * (dbB1*sin(dbLat) + dbC1*pow(sin(dbLat),3.0) + dbD1*pow(sin(dbLat),5.0) + dbE1*pow(sin(dbLat),7.0) + dbF1*pow(sin(dbLat),9.0) + dbG1*pow(sin(dbLat),11.0)));
}double CGeoEllipse::GetTrapeziumArea(double dbLat1,double dbLat2,double dbLon1,double dbLon2) const
{//检查参数是否合法assert(dbLat1 >= -90 && dbLat1 <= 90&& dbLon1 <= 180 && dbLon1 >= -180);assert(dbLat2 >= -90 && dbLat2 <= 90&& dbLon2 <= 180 && dbLon2 >= -180);//计算相关的系数dbLon1 = dbLon1 * acos(-1.0)/180.0;dbLon2 = dbLon2 * acos(-1.0)/180.0;double dbK = 2 * m_dbRadiusA*m_dbRadiusA * (1.0-m_dbE1*m_dbE1) * (dbLon2-dbLon1);double dbA = 1.0 + 0.5*m_dbE1*m_dbE1 + 3.0/8.0*pow(m_dbE1,4.0) + 5.0/16.0*pow(m_dbE1,6.0) + 630.0/2304.0*pow(m_dbE1,8.0);double dbB = 1.0/6.0*m_dbE1*m_dbE1 + 0.3*pow(m_dbE1,4.0) + 3.0/16.0*pow(m_dbE1,6.0) + 420.0/2304.0*pow(m_dbE1,8.0);double dbC = 3.0/80.0*pow(m_dbE1,4.0) + 1.0/16.0*pow(m_dbE1,6.0) + 180.0/2304.0*pow(m_dbE1,8.0);double dbD = 1.0/112.0*pow(m_dbE1,6.0) + 45.0/2304.0*pow(m_dbE1,8.0);double dbE = 5.0/2304.0*pow(m_dbE1,8.0);//计算纬度差，中间纬度等dbLat1 = dbLat1 * acos(-1.0)/180.0;dbLat2 = dbLat2 * acos(-1.0)/180.0;double dbDetLat = fabs(dbLat2-dbLat1);double dbMidLat = (dbLat1+dbLat2)/2.0;//计算结果return dbK * (dbA*sin(dbDetLat/2.0)*cos(dbMidLat) - dbB*sin(3*dbDetLat/2.0)*cos(3*dbMidLat)+ dbC*sin(5*dbDetLat/2.0)*cos(5*dbMidLat) - dbD*sin(7*dbDetLat/2.0)*cos(7*dbMidLat) + dbE*sin(9*dbDetLat/2.0)*cos(9*dbMidLat)) ;
}double CGeoEllipse::GetValueU(double dbLat)
{double dbTemp = (1.0-m_dbE1*sin(dbLat))/(1.0+m_dbE1*sin(dbLat));double dbExp = m_dbE1/2.0;double dbEquard = pow(dbTemp , dbExp);return tan(M_PI_4 + dbLat/2.0) * dbEquard;
}bool CGeoEllipse::BLH_XYZ(double dbLon,double dbLat,double dbHei,double& x,double &y,double &z)
{//参数合法性检查assert(dbLat >= -90 && dbLat <= 90);assert(dbLon >= -180 && dbLon <= 180);//计算卯酉圈曲率半径double dbN = GetPrimeRadius(dbLat);double B = dbLat * DEG_TO_RAD; //转换为弧度制double L = dbLon * DEG_TO_RAD; double H = dbHei;//计算空间直角坐标x = (dbN + H)*cos(B)*cos(L);y = (dbN + H)*cos(B)*sin(L);z = (dbN*(1.0-m_dbE1*m_dbE1)+ H)*sin(B);//转换成功return true;
}bool CGeoEllipse::XYZ_BLH(double X,double Y,double Z, double &dbLon,double &dbLat,double &dbHei)
{//double dDet = 0.0000000001;     //阈值
//#ifdef 0
//  //求解经度
//  dbLon = atan(Y/X);
//  if (X < 0)
//  {
//      dbLon += M_PI;
//  }
//
//  //求解纬度
//  double r = /*pow(X,2.0)+pow(Y,2.0)*/X*X + Y*Y;
//  double R = Z*Z + r;
//
//  double dbU = atan(Z/sqrt(r)*sqrt(1.0-m_dbE1*m_dbE1));
//
//  double dbTemp = (Z + m_dbE2*m_dbE2*m_dbRadiusB*pow(sin(dbU),3.0) ) /
//      (sqrt(r) - m_dbE1*m_dbE1*m_dbRadiusA*pow(cos(dbU),3.0));
//
//  dbLat = atan(dbTemp);
//
//  //求解高程
//  dbHei = sqrt(r)*cos(dbLat) + Z*sin(dbLat) -
//      m_dbRadiusA*(sqrt(1.0-m_dbE1*m_dbE1*pow(sin(dbLat),2.0)));
//
//  //转换为角度值
//  dbLon = (dbLon/acos(-1.0))*180;
//  dbLat = (dbLat/acos(-1.0))*180;
//
//  //转换成功
//  return true;
//
//#else//求解经度dbLon = atan2(Y,X);/*if (X < 0){dbLon += M_PI;}*///求解纬度double r = X*X + Y*Y;double R = Z*Z + r;//求解一些中间变量double dbPhi = atan(Z/sqrt(r));double dbU = ( m_dbRadiusB *Z * (1 + m_dbE2*m_dbE2*m_dbRadiusB/sqrt(R)) ) / (m_dbRadiusA*sqrt(r));dbU = atan(dbU);double dbTemp = (Z + m_dbE2*m_dbE2*m_dbRadiusB* sin(dbU) * sin(dbU) * sin(dbU) ) / (sqrt(r) - m_dbE1*m_dbE1*m_dbRadiusA*cos(dbU) * cos(dbU) * cos(dbU));dbLat = atan(dbTemp);//求解高程dbHei = sqrt(r)*cos(dbLat) + Z*sin(dbLat) - m_dbRadiusA*(sqrt(1.0-m_dbE1*m_dbE1*sin(dbLat) * sin(dbLat)));//转换为角度值dbLon = dbLon*RAD_TO_DEG;dbLat = dbLat*RAD_TO_DEG;return true;//#endif
}double CGeoEllipse::GetAveRadius() const
{return (2*m_dbRadiusA + m_dbRadiusB)/3;
}

希望能对大家有用，写这些主要是个人兴趣，其实有更好的Proj4这个开源库可以实现这些，或者GDAL里面也有。

地球椭球体基本要素的计算相关推荐

地球椭球体(Ellipsoid)、大地基准面(Datum)及地图投影(Projection)三者的基本概念
地球椭球体(Ellipsoid) 众所周知我们的地球表面是一个凸凹不平的表面,而对于地球测量而言,地表是一个无法用数学公式表达的曲面,这样的曲面不能作为测量和制图的基准面.假想一个扁率极小的椭圆,绕大 ...
大地水准面地球椭球体大地基准面地图投影理解
大地水准面地球椭球体大地基准面地图投影理解地球的自然表面是崎岖不平的,在地理课本上我们会看到对地球形状的描述:地球是一个两极稍扁,赤道略鼓的不规则球体. 大地水准面: 对不规则地球的第一层抽象 ...
（二）地理信息中对地球的描述-地球的大地水准面、地球椭球体、大地基准面
目录 1 概述 1.1 关于本文 1.2 进一步的认识地球 1.3 地球的三级逼近 2 大地水准面.地球椭球体.大地基准面 2.1 大地水准面 2.2 地球椭球体 2.3 大地基准面 1 概述 1.1 ...
电磁场与电磁波面电流和体电流磁感应强度的计算
1.面电流磁感应强度的计算 2.体电流磁感应强度的计算 2.6 矢量磁位的引入与计算磁通量的连续性矢量磁位的引入矢量磁位的计算 4 结论:穿过任意闭合曲面的磁通量恒为零.这就是磁通连续性原理.它 ...
地球椭球体与大地基准面的关系
地球椭球体与大地基准面的关系前言本文讲解什么是地球椭球体?什么是大地基准面?为什么需要地球椭球体?又为什么需要大地基准面? 一.地球表面几何模型有哪些? 为了研究地理空间,就必须对研究的对象--- ...
建筑幕墙单元体.组装件.零部件计算.查询器（Excel VBA版）
前言: 本篇是建筑幕墙零部件计算器(Excel版)的使用说明书. 零部件计算与查询器为建筑幕墙设计专用程序,适用于有查询及算料需求的幕墙设计或工厂算料人员使用.(程序亦适用于迭代递归的材料计算) 程 ...
结构体占用内存的计算
目录前言一.结构体二.计算方法总结前言大家都知道我们的数据类型如 int,char都有他们各自占据的内存大小,int 是4个字节,char 是一个字节,当我们接触到结构体这个东西的时候,或 ...
第10章结构体01——结构体字节大小的计算
1.该篇笔记详见C提高笔记(传智播客) 文章目录博文01:(常考面试题) 三步解决C语言中struct字节对齐问题,结构体的字节大小问题第1步.先确定结构体实际对齐单位, 第2步.除结构体的第一个 ...
白话地图投影之初识地球
本文是Koala带你进入GIS世界的开篇,Koala打算用简单通俗的语言为大家介绍地图投影,帮助GISer理解地图投影的概念.作为进入GIS世界多年的老鸟,Koala也是在不断的实战中才真正理解和掌握 ...

地球椭球体基本要素的计算

地球椭球体基本要素的计算相关推荐

最新文章

热门文章