根据经纬度计算多边形面积

高冷的哥哥 · 发表于 2017-7-27 16:54

计算方法比较简单，主要是求出多边形外接矩形已米为单位面积和已经纬度为单位面积比值，然后用这个比值乘以多边形经纬度为单位面积，即可得出这个多边形以米为单位面积。

double GetArea(constvector<Coordinate>&ls)
{
if(ls.size() < 4)
   return0;

double sum =0;

for (size_ti=0; i<ls.size()-1; ++i)
{
   const Coordinate& p = ls[i];
   constCoordinate& q = ls[i+1];
   sum += (p.x + q.x) * (q.y - p.y);
}

returnsum/2;
}

double GetPrjArea(constvector<Coordinate>&ls)
{
if(ls.size() < 4)
   return0;

double dArea= GetArea(ls);
dArea =abs(dArea);
if (dArea ==0)
   return0;

double xmin,ymin, xmax, ymax;
xmin = xmax= ls[0].x;
ymax = ymin= ls[0].y;

for (size_ti=1; i<ls.size(); ++i)
{
   constCoordinate& p = ls[i];
   xmin =min(xmin, p.x);
   ymin =min(ymin, p.y);
   xmax =max(xmax, p.x);
   ymax =max(ymax, p.y);
}

Coordinatep1, p2;
p1.x =xmin;
p1.y =(ymin+ymax)/2;
p2.x =xmax;
p2.y =(ymin+ymax)/2;
double dx =GetPrjDistance(p1, p2);

p1.x = p2.x= xmin;
p1.y =ymin;
p2.y =ymax;
double dy =GetPrjDistance(p1, p2);
dy *=dx;

dx =(xmax-xmin)*(ymax-ymin);
dy /=dx;
dArea *=dy;
returndArea;
}

做了简单的测试，用此方法计算出来的面积和投影变换后计算的面积误差大约为1/1000，基本上满足一些要求精度不是很高的应用。