Polygon Area Acreage (c + + code) on the Earth's ellipsoidal surface

Yesterday when the test was suddenly found in the previous product written in the Earth's ellipsoid area of the calculation of the code a bit of a problem, so today is a complete correction, from the Qgis to key out the code to rewrite with C + +, the new code can be more accurate calculation of the surface of the ellipsoid polygon area, This basic function is very important to the measurement of the area in the function of space acreage, which is shared here for your reference and even used directly.

The header files are as follows:

/*** @file distancearea.h* interface file for calculating polygon area on @brief ellipsoid * @details * @author zxg* @date May 15, 2015 * @version 1.0.0.1* @par Copyright (c): Zhou Xuquang * @par histor y:*/#ifndef __distancearea_h_c0c6c4c8_b1c4_49c3_87e4_eef1abc862f0__#define __distancearea_h_c0c6c4c8_b1c4_49c3_    87e4_eef1abc862f0__class Distancearea{public:distancearea (); ~distancearea ();/*** @brief setellipsepara set parameters for calculated area * @param [in] double a long half axis * @param [in] double b short half axis * @return void* @aut Hor zxg* @date May 15, 2015 * @note */void Setellipsepara (double a,double b)/*** @brief Computepolygonarea Calculate the area of the polygon on the Earth's ellipsoid * @ Param[in] Const Double *PADX x-coordinate array * @param [in] const double* Pady y-coordinate array * @param [in] int ncount point number * @return Double return value is polygon Product * @author zxg* @date May 15, 2015 * @note */double computepolygonarea (const double *padx,const double* pady,int ncount);p R    Ivate:double msemimajor, Msemiminor, minvflattening;    Double GetQ (double x); Double Getqbar (double x); void Computeareainit ();    Area calculation temporary variable 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__

The implementation code is as follows:

#include "DistanceArea.h" #define DEG2RAD (x) ((x) *m_pi/180) Distancearea::D Istancearea () {//}distancearea::~ Distancearea () {}void Distancearea::setellipsepara (double a,double b) {msemimajor = A;msemiminor = B;//mInvFlattening = Msemimajorcomputeareainit ();} 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 + C OSX2 * (M_qbarb + cosx2 * (M_qbarc + cosx2 * m_qbard)));} void Distancearea::computeareainit () {Double a2 = (msemimajor * msemimajor);d ouble e2 = 1-(A2/(Msemiminor * Msemimi Nor));d ouble 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_q  B = (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 Q Bar1, qbar2;if (NULL = = PADX | | NULL = = pady) {return 0;} if (Ncount < 3) {return 0;} Double x2 = Deg2rad (padx[ncount-1]);d ouble y2 = Deg2rad (pady[ncount-1]); QBAR2 = Getqbar (y2);d ouble 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;} 

The main function is Computepolygonarea, the calculated area unit is square meters.

The test examples are as follows:

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);D Istancearea Area;area. Setellipsepara (6378140.0,6356755.0);d ouble AAA = Area.computepolygonarea (&vecx[0],&vecy[0],vecy.size ());

Tested to meet the requirements.

Polygon Area Acreage (c + + code) on the Earth's ellipsoidal surface