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

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

一般说来，经纬度坐标多边形面积指的是球面多边形面积。我曾经在作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

{

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

{

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;

}

----------------------------------------------------------------------------------------------------------------------------

到此结束，敬请批评指正