这两天看到有人询问WGS84下的多边形面积计算。一般说来，经纬度坐标多边形面积指的是球面多边形面积。我曾经在作ArcIMS项目时写了一个Javascript函数，特贴出来，大家需要时可以参考。为方便大家直接调用，我做了简单修改，如果有问题，请批评指正。还需要注意的是，该函数不适用于自交叉多边形。需要注意的是其中PointX和PointY虽是多边形的节点坐标，但本公式已将首尾相同的那个节点剔除掉了一个，作为一个点纳入计算。
[Copy to clipboard] [ - ]CODE:
// calculate Area
function calcArea(PointX,PointY,MapUnits) {
var Count = PointX.length
if (Count>2) {
var mtotalArea = 0;
if (MapUnits=="DEGREES")//经纬度坐标下的球面多边形
{
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 = 6378000;
for(i=0;i
if(i==0)
{
LowX = PointX[Count-1] * Math.PI / 180;
LowY = PointY[Count-1] * Math.PI / 180;
MiddleX = PointX[0] * Math.PI / 180;
MiddleY = PointY[0] * Math.PI / 180;
HighX = PointX[1] * Math.PI / 180;
HighY = PointY[1] * Math.PI / 180;
}
else if(i==Count-1)
{
LowX = PointX[Count-2] * Math.PI / 180;
LowY = PointY[Count-2] * Math.PI / 180;
MiddleX = PointX[Count-1] * Math.PI / 180;
MiddleY = PointY[Count-1] * Math.PI / 180;
HighX = PointX[0] * Math.PI / 180;
HighY = PointY[0] * Math.PI / 180;
}
else
{
LowX = PointX[i-1] * Math.PI / 180;
LowY = PointY[i-1] * Math.PI / 180;
MiddleX = PointX[i] * Math.PI / 180;
MiddleY = PointY[i] * Math.PI / 180;
HighX = PointX[i+1] * Math.PI / 180;
HighY = PointY[i+1] * 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 ++;
//Sum +=2*Math.PI-AngleCos;
}
}
if(Sum1>Sum2){
Sum = Sum1+(2*Math.PI*Count2-Sum2);
}
else{
Sum = (2*Math.PI*Count1-Sum1)+Sum2;
}
//平方米
mtotalArea = (Sum-(Count-2)*Math.PI)*Radius*Radius;
}
else { //非经纬度坐标下的平面多边形
var i,j;
var j;
var p1x,p1y;
var p2x,p2y;
for(i=Count-1, j=0; j
p1x = PointX[i];
p1y = PointY[i];
p2x = PointX[j];
p2y = PointY[j];
mtotalArea +=p1x*p2y-p2x*p1y;
}
mtotalArea /= 2.0;
}
return mtotalArea;
}
return;
}
不太好注释，具体原理请参考古人的定理：
球面多边形计算面积的关键在于计算多边形所有角的度数.
对于球面n边形,所有角的和为S,球的半径为R,那么其面积就是
R^2*(S-(n-2)*Pi)