地球椭球面上多边形面积量算（C++代码）

昨天突然测试的时候发现以前产品中写的地球椭球面上面积计算的代码有点问题，于是今天就彻底修正，从QGIS中抠出代码来用C++重写了一下，新代码可以比较准确计算椭球面上多边形的面积，这个基础函数对空间量算功能中的面积量测非常重要，在这里共享出来供大家参考甚至直接拿过去用。

头文件如下：

/**
* @file              DistanceArea.h
* @brief             椭球面上计算多边形面积的接口文件
* @details
* @author            zxg
* @date              2015年5月15日
* @version           1.0.0.1
* @par               Copyright (c):周旭光
* @par               History:
*/
#ifndef __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__
#define __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__

class DistanceArea
{
public:
	DistanceArea();

    ~DistanceArea();

	/**
	* @brief SetEllipsePara 设置计算面积的参数
	* @param[in] double a 长半轴
	* @param[in] double b 短半轴
	* @return void
	* @author zxg
	* @date 2015年5月15日
	* @note
	*/
	void SetEllipsePara(double a,double b);

	/**
	* @brief ComputePolygonArea 计算地球椭球面上多边形的面积
	* @param[in]  const double *padX X坐标数组
	* @param[in] const double* padY Y坐标数组
	* @param[in] int nCount 点的个数
	* @return double 返回值是面积
	* @author zxg
	* @date 2015年5月15日
	* @note
	*/
	double ComputePolygonArea( const double *padX,const double* padY,int nCount );

private:

    double mSemiMajor, mSemiMinor, mInvFlattening;

    double GetQ( double x );
    double GetQbar( double x );

	void ComputeAreaInit();

    // 面积计算临时变量

    double m_QA, m_QB, m_QC;
    double m_QbarA, m_QbarB, m_QbarC, m_QbarD;
    double m_AE;  /* a^2(1-e^2) */
    double m_Qp;
    double m_E;
    double m_TwoPI;

};

#endif // end of __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__

实现代码如下：

#include "DistanceArea.h"

#define DEG2RAD(x)    ((x)*M_PI/180)

DistanceArea::DistanceArea()
{
	//
}

DistanceArea::~DistanceArea()
{
}

void DistanceArea::SetEllipsePara(double a,double b)
{
	mSemiMajor = a;
	mSemiMinor = b;
	//mInvFlattening = mSemiMajor

	ComputeAreaInit();
}

double DistanceArea::GetQ( double x )
{
	double sinx, sinx2;

	sinx = sin( x );
	sinx2 = sinx * sinx;

	return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
}

double DistanceArea::GetQbar( double x )
{
	double cosx, cosx2;

	cosx = cos( x );
	cosx2 = cosx * cosx;

	return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
}

void DistanceArea::ComputeAreaInit()
{
	double a2 = ( mSemiMajor * mSemiMajor );
	double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
	double e4, e6;

	m_TwoPI = M_PI + M_PI;

	e4 = e2 * e2;
	e6 = e4 * e2;

	m_AE = a2 * ( 1 - e2 );

	m_QA = ( 2.0 / 3.0 ) * e2;
	m_QB = ( 3.0 / 5.0 ) * e4;
	m_QC = ( 4.0 / 7.0 ) * e6;

	m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4  - ( 4.0 / 7.0 ) * e6;
	m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4  + ( 4.0 / 7.0 ) * e6;
	m_QbarC =                     - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
	m_QbarD = ( 4.0 / 49.0 ) * e6;

	m_Qp = GetQ( M_PI / 2 );
	m_E  = 4 * M_PI * m_Qp * m_AE;
	if ( m_E < 0.0 )
	m_E = -m_E;
}

double DistanceArea::ComputePolygonArea( const double *padX,const double* padY,int nCount )
{
	double x1, y1, dx, dy;
	double Qbar1, Qbar2;

	if (NULL == padX || NULL == padY)
	{
		return 0;
	}

	if (nCount < 3)
	{
		return 0;
	}

	double x2 = DEG2RAD( padX[nCount-1] );
	double y2 = DEG2RAD( padY[nCount-1] );
	Qbar2 = GetQbar( y2 );

	double area = 0.0;

	for ( int i = 0; i < nCount; i++ )
	{
		x1 = x2;
		y1 = y2;
		Qbar1 = Qbar2;

		x2 = DEG2RAD( padX[i] );
		y2 = DEG2RAD( padY[i] );
		Qbar2 = GetQbar( y2 );

		if ( x1 > x2 )
		  while ( x1 - x2 > M_PI )
			x2 += m_TwoPI;
		else if ( x2 > x1 )
		  while ( x2 - x1 > M_PI )
			x1 += m_TwoPI;

		dx = x2 - x1;
		area += dx * ( m_Qp - GetQ( y2 ) );

		if (( dy = y2 - y1 ) != 0.0 )
		  area += dx * GetQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
	}
	if (( area *= m_AE ) < 0.0 )
		area = -area;

	if ( area > m_E )
		area = m_E;
	if ( area > m_E / 2 )
		area = m_E - area;

	return area;
}

主要的函数是ComputePolygonArea，计算出来的面积单位是平方米。

测试示例如下：

std::vector<double> vecX;
	std::vector<double> vecY;
	vecX.push_back(116.0120);
	vecX.push_back(116.0121);
	vecX.push_back(116.01205);

	vecY.push_back(40.004);
	vecY.push_back(40.004);
	vecY.push_back(40.0041);

	DistanceArea area;
	area.SetEllipsePara(6378140.0,6356755.0);
	double aaa = area.ComputePolygonArea(&vecX[0],&vecY[0],vecY.size());

经过测试，可以满足要求。