通过坐标系求覆盖物面积
首先先供上原文链接：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的相关求面积代码，供大家参考
function getArea(polygon) {
var pts = new Array();
if (polygon instanceof BMap.Polygon) {
pts = polygon.getPath();
}
console.log(pts.length);
// 检查类型：既不是百度类型的范围又不是数组类型的数据，直接返回0
//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);
}