通过坐标系求覆盖物面积

首先先供上原文链接：https://blog.csdn.net/neil89/article/details/49331641

这是我在项目上遇到的问题，在地图上标注多边形覆盖物，同时求出覆盖物的面积。在查找方法的过程中，首先了解了百度api中的通过坐标获取面积的方式，解读代码发现求解方式是依照地球是曲面的这个标准计算的面积，当我将api中的就是方法拷贝出来自己封装调用后发现，面积大于2000平方米的区域面积大小计算误差较小，但是当覆盖较小区域时计算的误差特别大，甚至会出现负数的情况。针对这一概念，我进行了深度思考，结合百度参考发现，有可能是在曲面计算过程中向量的方向错误导致数据有误，还有可能曲率问题导致较小面积的数据不准确问题。通过进一步的思考，我想较小的面积可以看做是近似于平面多边形的方式计算（这就可以做两个方法，分别是大面积和小面积的不同求解方法），沿着这样的思路，回忆起可以将多边形分成多个三角形然后分别通过海伦公式（坐标间距离可以求出）求出三角形面积并求和得到多边形的面积。（这一多边形依然不包括有交叉边界的多边形）（这只是我个人的想法，实践起来可能要花一些时间，下面是百度到的一篇一样分为两种面积求法的代码）

【网站上的算法比较复杂，我还没有钻研过深，如果有理解的可以反馈大致思路】

下面是大神代码原文供大家欣赏：

var earthRadiusMeters = 6371000.0;
var metersPerDegree = 2.0 * Math.PI * earthRadiusMeters / 360.0;
var radiansPerDegree = Math.PI / 180.0;
var degreesPerRadian = 180.0 / Math.PI;
var pointArr;

$(document).ready(function() {
pointArr = new Array();
b();
});

function calculateArea(points) {
if (points.length > 2) {
var areaMeters2 = PlanarPolygonAreaMeters2(points);
if (areaMeters2 > 1000000.0) {
areaMeters2 = SphericalPolygonAreaMeters2(points);
alert("面积为" + areaMeters2 + "平方米");
}
}
}

/球面多边形面积计算/
function SphericalPolygonAreaMeters2(points) {
var totalAngle = 0;
for (var i = 0; i < points.length; i++) {
var j = (i + 1) % points.length;
var k = (i + 2) % points.length;
totalAngle += Angle(points[i], points[j], points[k]);
}
var planarTotalAngle = (points.length - 2) * 180.0;
var sphericalExcess = totalAngle - planarTotalAngle;
if (sphericalExcess > 420.0) {
totalAngle = points.length * 360.0 - totalAngle;
sphericalExcess = totalAngle - planarTotalAngle;
} else if (sphericalExcess > 300.0 && sphericalExcess < 420.0) {
sphericalExcess = Math.abs(360.0 - sphericalExcess);
}
return sphericalExcess * radiansPerDegree * earthRadiusMeters * earthRadiusMeters;
}

/角度/
function Angle(p1, p2, p3) {
var bearing21 = Bearing(p2, p1);
var bearing23 = Bearing(p2, p3);
var angle = bearing21 - bearing23;
if (angle < 0) {
angle += 360;
}
return angle;
}

/方向/
function Bearing(from, to) {
var lat1 = from[1] * radiansPerDegree;
var lon1 = from[0] * radiansPerDegree;
var lat2 = to[1] * radiansPerDegree;
var lon2 = to[0] * radiansPerDegree;
var angle = -Math.atan2(Math.sin(lon1 - lon2) * Math.cos(lat2), Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(lon1 - lon2));
if (angle < 0) {
angle += Math.PI * 2.0;
}
angle = angle * degreesPerRadian;
return angle;
}

/平面多边形面积/
function PlanarPolygonAreaMeters2(points) {
var a = 0;
for (var i = 0; i < points.length; ++i) {
var j = (i + 1) % points.length;
var xi = points[i][0] * metersPerDegree * Math.cos(points[i][1] * radiansPerDegree);
var yi = points[i][1] * metersPerDegree;
var xj = points[j][0] * metersPerDegree * Math.cos(points[j][1] * radiansPerDegree);
var yj = points[j][1] * metersPerDegree;
a += xi * yj - xj * yi;
}
return Math.abs(a / 2);
}

function b() {
var s = "112.523197631836,37.868892669677734;112.5170669555664,37.8605842590332;112.52099609375,37.849857330322266;112.54137420654297,37.8512732521875;112.5351180302734,37.858699798583984";
var s1 = new Array()
s1 = s.split(";");
for (var i = 0; i < s1.length; i++) {
var ss = s1[i];
var temp = ss.split(",");
var point = new Array();
point.push(Number(temp[0]), Number(temp[1]));
pointArr.push(point);
}
calculateArea(pointArr);
}

后记：后面会放百度api的相关求面积代码，供大家参考
//百度api提供方法
// 检查类型：既不是百度类型的范围又不是数组类型的数据，直接返回0

    //console.log(pts.length);

    
     

    //如果是百度类型的，得到点集合，不是的话，另说。
    function getArea(polygon) {
    var pts = new Array();
    if (polygon instanceof BMap.Polygon) {
        pts = polygon.getPath();
    }
    if (!(polygon instanceof BMap.Polygon) && !(polygon instanceof Array)) {
       return 0;
    }

    ////判断数组的长度，如果是小于3的话，不构成面，无法计算面积
    if (pts.length < 3) {
        return 0;
    }

    var totalArea = 0;// 初始化总面积
    var LowX = 0.0;
    var LowY = 0.0;
    var MiddleX = 0.0;
    var MiddleY = 0.0;
    var HighX = 0.0;
    var HighY = 0.0;
    var AM = 0.0;
    var BM = 0.0;
    var CM = 0.0;
    var AL = 0.0;
    var BL = 0.0;
    var CL = 0.0;
    var AH = 0.0;
    var BH = 0.0;
    var CH = 0.0;
    var CoefficientL = 0.0;
    var CoefficientH = 0.0;
    var ALtangent = 0.0;
    var BLtangent = 0.0;
    var CLtangent = 0.0;
    var AHtangent = 0.0;
    var BHtangent = 0.0;
    var CHtangent = 0.0;
    var ANormalLine = 0.0;
    var BNormalLine = 0.0;
    var CNormalLine = 0.0;
    var OrientationValue = 0.0;
    var AngleCos = 0.0;
    var Sum1 = 0.0;
    var Sum2 = 0.0;
    var Count2 = 0;
    var Count1 = 0;
    var Sum = 0.0;
    var Radius = 6378137.0;// ,WGS84椭球半径
    var Count = pts.length;
    for (var i = 0; i < Count; i++) {
        if (i == 0) {
            LowX = pts[Count - 1].lng * Math.PI / 180;
            LowY = pts[Count - 1].lat * Math.PI / 180;
            MiddleX = pts[0].lng * Math.PI / 180;
            MiddleY = pts[0].lat * Math.PI / 180;
            HighX = pts[1].lng * Math.PI / 180;
            HighY = pts[1].lat * Math.PI / 180;
        } else if (i == Count - 1) {
            LowX = pts[Count - 2].lng * Math.PI / 180;
            LowY = pts[Count - 2].lat * Math.PI / 180;
            MiddleX = pts[Count - 1].lng * Math.PI / 180;
            MiddleY = pts[Count - 1].lat * Math.PI / 180;
            HighX = pts[0].lng * Math.PI / 180;
            HighY = pts[0].lat * Math.PI / 180;
        } else {
            LowX = pts[i - 1].lng * Math.PI / 180;
            LowY = pts[i - 1].lat * Math.PI / 180;
            MiddleX = pts[i].lng * Math.PI / 180;
            MiddleY = pts[i].lat * Math.PI / 180;
            HighX = pts[i + 1].lng * Math.PI / 180;
            HighY = pts[i + 1].lat * Math.PI / 180;
        }
        AM = Math.cos(MiddleY) * Math.cos(MiddleX);
        BM = Math.cos(MiddleY) * Math.sin(MiddleX);
        CM = Math.sin(MiddleY);
        AL = Math.cos(LowY) * Math.cos(LowX);
        BL = Math.cos(LowY) * Math.sin(LowX);
        CL = Math.sin(LowY);
        AH = Math.cos(HighY) * Math.cos(HighX);
        BH = Math.cos(HighY) * Math.sin(HighX);
        CH = Math.sin(HighY);
        CoefficientL = (AM * AM + BM * BM + CM * CM) / (AM * AL + BM * BL + CM * CL);
        CoefficientH = (AM * AM + BM * BM + CM * CM) / (AM * AH + BM * BH + CM * CH);
        ALtangent = CoefficientL * AL - AM;
        BLtangent = CoefficientL * BL - BM;
        CLtangent = CoefficientL * CL - CM;
        AHtangent = CoefficientH * AH - AM;
        BHtangent = CoefficientH * BH - BM;
        CHtangent = CoefficientH * CH - CM;
        AngleCos = (AHtangent * ALtangent + BHtangent * BLtangent + CHtangent * CLtangent) / (Math.sqrt(AHtangent * AHtangent + BHtangent * BHtangent + CHtangent * CHtangent) * Math.sqrt(ALtangent * ALtangent + BLtangent * BLtangent + CLtangent * CLtangent));
        AngleCos = Math.acos(AngleCos);
        ANormalLine = BHtangent * CLtangent - CHtangent * BLtangent;
        BNormalLine = 0 - (AHtangent * CLtangent - CHtangent * ALtangent);
        CNormalLine = AHtangent * BLtangent - BHtangent * ALtangent;
        if (AM != 0)
            OrientationValue = ANormalLine / AM;
        else if (BM != 0)
            OrientationValue = BNormalLine / BM;
        else
            OrientationValue = CNormalLine / CM;
        if (OrientationValue > 0) {
            Sum1 += AngleCos;
            Count1++;
        } else {
            Sum2 += AngleCos;
            Count2++;
        }
    }

    var tempSum1, tempSum2;
    tempSum1 = Sum1 + (2 * Math.PI * Count2 - Sum2);
    tempSum2 = (2 * Math.PI * Count1 - Sum1) + Sum2;
    if (Sum1 > Sum2) {
        if ((tempSum1 - (Count - 2) * Math.PI) < 1)
            Sum = tempSum1;
        else
            Sum = tempSum2;
    } else {
        if ((tempSum2 - (Count - 2) * Math.PI) < 1)
            Sum = tempSum2;
        else
            Sum = tempSum1;
    }
    totalArea = (Sum - (Count - 2) * Math.PI) * Radius * Radius;

    console.log((totalArea.toFixed(4))/10000); // 返回总面积

    return totalArea.toFixed(4);
}