使用qt计算地理平面下椭圆坐标的不规则不封闭多边形面积

实例代码

#include <QApplication>
#include <QPainter>
#include <QPixmap>
#include <QPolygon>
#include <QPen>
#include "DistanceArea.h"

int main(int argc, char *argv[])
{
	//    QApplication a(argc, argv);
	//    QPixmap pixmap(1024,800);

	double dSize = 0;//最终计算的多边形面积，单位平方米
	std::vector<double> vecX;
	std::vector<double> vecY;
	DistanceArea distanceAreaTool;//地理面积计算工具

								  //    QPen pen(QColor(Qt::blue));
								  //    pen.setWidth(20);

								  //    QPainter* painter = new QPainter(&pixmap);
								  //    painter->setPen(pen);
	QPolygonF polygonF;
	polygonF.append(QPointF(112857055, 34329828));
	polygonF.append(QPointF(115098266, 34293530));
	polygonF.append(QPointF(113378906, 32911873));
	polygonF.append(QPointF(114049072, 35456195));
	polygonF.append(QPointF(114884033, 32588477));
	polygonF.append(QPointF(112857055, 34329828));
	//需要计算的面积的相交不封闭不规则多边形
	QPainterPath pathNeedCalcAreaPolygon;
	pathNeedCalcAreaPolygon.addPolygon(polygonF);
	//    painter->fillPath(pathNeedCalcAreaPolygon, QBrush(painter->pen().color(), Qt::DiagCrossPattern));

	//多边形的最大外界矩形
	QRectF rectfNeedCalcAreaPolygonMaxBoundingBox = pathNeedCalcAreaPolygon.boundingRect();
	QPainterPath pathPolygonMaxBoundingBox;
	pathPolygonMaxBoundingBox.addRect(rectfNeedCalcAreaPolygonMaxBoundingBox);

	//多边形的外界矩形与多边形的补集
	QList<QPolygonF> listSubtracted = pathPolygonMaxBoundingBox.subtracted(pathNeedCalcAreaPolygon).toSubpathPolygons();

	//绘制多边形，需要经过坐标转换转成像素点
	//    for (int i = 0; i < listSubtracted.size(); ++i)
	//    {
	//        pen.setColor(QColor(qrand() % 255, qrand() % 255, qrand() % 255));
	//        painter->setPen(pen);
	//        painter->drawPolygon(listSubtracted.at(i));
	//    }
	//    pixmap.save("polygon.png");

	//外接矩形的面积
	double dRectfArea = 0;
	vecX.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topLeft().x() / 1000000.0);
	vecX.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topRight().x() / 1000000.0);
	vecX.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.bottomRight().x() / 1000000.0);
	vecX.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.bottomLeft().x() / 1000000.0);
	vecX.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topLeft().x() / 1000000.0);

	vecY.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topLeft().y() / 1000000.0);
	vecY.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topRight().y() / 1000000.0);
	vecY.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.bottomRight().y() / 1000000.0);
	vecY.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.bottomLeft().y() / 1000000.0);
	vecY.push_back(rectfNeedCalcAreaPolygonMaxBoundingBox.topLeft().y() / 1000000.0);

	dRectfArea = distanceAreaTool.ComputePolygonArea(&vecX[0], &vecY[0], vecX.size());

	//补集多边形的面积
	double dOtherPolygonArea = 0;
	for (int i = 0; i < listSubtracted.size(); ++i)
	{
		vecX.clear();
		vecY.clear();

		for (int j = 0; j < listSubtracted.at(i).size(); ++j)
		{
			vecX.push_back(listSubtracted.at(i).at(j).x() / 1000000.0);
			vecY.push_back(listSubtracted.at(i).at(j).y() / 1000000.0);
		}
		dOtherPolygonArea += distanceAreaTool.ComputePolygonArea(&vecX[0], &vecY[0], vecX.size());
	}


	dSize = dRectfArea - dOtherPolygonArea;

	return 0;
}

椭圆封闭多边形面积计算

头文件

#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"
#include <qmath.h>
#define DEG2RAD(x)    ((x)*M_PI/180)


DistanceArea::DistanceArea()
{
    SetEllipsePara(6378140.0,6356755.0);
}


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 = qSin( 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 = qCos( 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;
}