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

经纬度坐标下的球面多边形面积计算公式前段时间，想做一个根据地球经纬度坐标计算地球表面面积的软件，查阅大量资料，找到如下方法，仅供参考。一般说来，经纬度坐标多边形面积指的是球面多边形面积。我曾经在作ArcIMS项目时写了一个Javascript函数，特贴出来，大家需要时可以参考。为方便大家直接调用，我做了简单修改，如果有问题，请批评指正。还需要注意的是，该函数不适用于自交叉多边形。不太好注释，具体原理请参考前人的定理：球面多边形计算面积的关键在于计算多边形所有角的度数.对于球面n边形，所有角的和为S，球的半径为R，那么其面积就是球面面积=R^2*(S-(n-2)*Pi)---------------------------------------------------------------------------------------------------------------------------------CODE://calculateAreafunctioncalcArea(PointX,PointY,MapUnits){varCount=PointX.lengthif(Count3){//至少3个点varmtotalArea=0;if((PointX[0]!=PointX[Count-1])||(PointY[0]!=PointY[Count-1]))//第1个点与最后1个点不重合{return;}if(MapUnits==DEGREES)//经纬度坐标下的球面多边形//////////////////degrees度数{varLowX=0.0;varLowY=0.0;varMiddleX=0.0;varMiddleY=0.0;varHighX=0.0;varHighY=0.0;varAM=0.0;varBM=0.0;varCM=0.0;varAL=0.0;varBL=0.0;varCL=0.0;varAH=0.0;varBH=0.0;varCH=0.0;varCoefficientL=0.0;//Coefficient系数varCoefficientH=0.0;varALtangent=0.0;//tangent切线varBLtangent=0.0;varCLtangent=0.0;varAHtangent=0.0;varBHtangent=0.0;varCHtangent=0.0;varANormalLine=0.0;//NormalLine法线varBNormalLine=0.0;varCNormalLine=0.0;varOrientationValue=0.0;//OrientationValue方向值varAngleCos=0.0;//余弦角varSum1=0.0;varSum2=0.0;varCount2=0;varCount1=0;varSum=0.0;varRadius=6378000;//半径for(i=0;iCount;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;}elseif(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;elseif(BM!=0)OrientationValue=BNormalLine/BM;elseOrientationValue=CNormalLine/CM;if(OrientationValue0){Sum1+=AngleCos;Count1++;}else{Sum2+=AngleCos;Count2++;//Sum+=2*Math.PI-AngleCos;}}if(Sum1Sum2){Sum=Sum1+(2*Math.PI*Count2-Sum2);}else{Sum=(2*Math.PI*Count1-Sum1)+Sum2;}//平方米mtotalArea=(Sum-(Count-2)*Math.PI)*Radius*Radius;}else{//非经纬度坐标下的平面多边形vari,j;varj;varp1x,p1y;varp2x,p2y;for(i=Count-1,j=0;jCount;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/=2.0;}returnmtotalArea;}return;}----------------------------------------------------------------------------------------------------------------------------到此结束，敬请批评指正