经纬度坐标下的球面多边形面积计算公式

经纬度坐标下的球面多边形面积计算公式
前段时间，想做一个根据地球经纬度坐标计算地球表面面积的软件，查阅大量资料，找到如下方法，仅供参考。

一般说来，经纬度坐标多边形面积指的是球面多边形面积。

我曾经在作ArcIMS项目时写了一个Javascript函数，特贴出来，大家需要时可以参考。

为方便大家直接调用，我做了简单修改，如果有问题，请批评指正。

还需要注意的是，该函数不适用于自交叉多边形。

不太好注释，具体原理请参考前人的定理：
球面多边形计算面积的关键在于计算多边形所有角的度数.对于球面n边形，所有角的和为S，球的半径为R，那么其面积就是
---------------------------------------------------------------------------------------------------------------------------------
CODE:
// calculate Area
function calcArea(PointX,PointY,MapUnits) {
var Count =
if (Count>3) {//至少3个点
var mtotalArea = 0;
if((PointX[0]!=PointX[Count-1])||(PointY[0]!=PointY[Count-1]))
//第1个点与最后1个点不重合
{
return;
}
if (MapUnits=="DEGREES")
//经纬度坐标下的球面多边形 //////////////////degrees度数
{
var LowY=;
var MiddleX=;
var MiddleY=;
var HighX=;
var HighY=;
var AM = ;
var BM = ;
var CM = ;
var AL = ;
var BL = ;
var CL = ;
var AH = ;
var BH = ;
var CH = ;
var CoefficientL = ;//Coefficient系数
var CoefficientH = ;
var ALtangent = ; //tangent切线
var BLtangent = ;
var CLtangent = ;
var AHtangent = ;
var BHtangent = ;
var CHtangent = ;
var ANormalLine = ; //NormalLine法线
var BNormalLine = ;
var CNormalLine = ;
var OrientationValue = ; //Orientation Value方向值 var AngleCos = ;//余弦角
var Sum1 = ;
var Sum2 = ;
var Count2 = 0;
var Count1 = 0;
var Radius = 6378000; //半径
for(i=0;i<Count;i++)
{
if(i==0)
{
LowX = PointX[Count-1] * / 180;//换算成弧度
LowY = PointY[Count-1] * / 180;
MiddleX = PointX[0] * / 180;
MiddleY = PointY[0] * / 180;
HighX = PointX[1] * / 180;
HighY = PointY[1] * / 180;
}
else if(i==Count-1)
{
LowX = PointX[Count-2] * / 180;
LowY = PointY[Count-2] * / 180;
MiddleX = PointX[Count-1] * / 180;
MiddleY = PointY[Count-1] * / 180;
HighX = PointX[0] * / 180;
HighY = PointY[0] * / 180;
}
else
{
LowX = PointX[i-1] * / 180;
LowY = PointY[i-1] * / 180;
MiddleX = PointX[i] * / 180;
MiddleY = PointY[i] * / 180;
HighX = PointX[i+1] * / 180;
HighY = PointY[i+1] * / 180; }
AM = (MiddleY) * (MiddleX);
BM = (MiddleY) * (MiddleX);
CM = (MiddleY);
AL = (LowY) * (LowX);
BL = (LowY) * (LowX);
CL = (LowY);
AH = (HighY) * (HighX);
BH = (HighY) * (HighX);
CH = (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)/(AHtangent * AHtangent + BHtangent * BHtangent +CHtangent * CHtangent) * (ALtangent * ALtangent + BLtangent * BLtangent +CLtangent * CLtangent));
AngleCos = (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*;
}
}
if(Sum1>Sum2){
Sum = Sum1+(2**Count2-Sum2);
}
else{
Sum = (2**Count1-Sum1)+Sum2;
}
//平方米
mtotalArea = (Sum-(Count-2)**Radius*Radius; }
else { //非经纬度坐标下的平面多边形
var i,j;
var j;
var p1x,p1y;
var p2x,p2y;
for(i=Count-1, j=0; j<Count; i=j, j++)
{
if(i==Count-1)
{
p1x = mX;
p1y = mY;
}
else
{
p1x = PointX[i];
p1y = PointY[i]; }
if(j==Count-1)
{
p2x = mX;
p2y = mY;
}
else
{
p2x = PointX[j];
p2y = PointY[j]; }
mtotalArea +=p1x*p2y-p2x*p1y;
}
mtotalArea /= ;
}
return mtotalArea;
}
return;
}
----------------------------------------------------------------------------------------------------------------------------
到此结束，敬请批评指正。