根據經緯度計算多邊形面積

計算方法比較簡單，主要是求出多邊形外接矩形已米為單位面積和已經緯度為單位面積比值，然后用這個比值乘以多邊形經緯度為單位面積，即可得出這個多邊形以米為單位面積。

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，基本上滿足一些要求精度不是很高的應用。