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

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

double GetArea(const vector<Coordinate>& ls)
{
    if (ls.size() < 4)
        return 0;

    double sum = 0;

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

    return sum/2;
}

double GetPrjArea(const vector<Coordinate> &ls)
{
    if (ls.size() < 4)
        return 0;

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

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

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

    Coordinate p1, 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;
    return dArea;
}

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