|
阅读：4476回复：6
经纬度坐标下的球面多边形面积计算公式
经纬度坐标下的球面多边形面积计算公式<BR><BR>
<DIV >这两天看到有人询问WGS84下的多边形面积计算。一般说来，经纬度坐标多边形面积指的是球面多边形面积。我曾经在作ArcIMS项目时写了一个Javascript函数，特贴出来，大家需要时可以参考。为方便大家直接调用，我做了简单修改，如果有问题，请批评指正。还需要注意的是，该函数不适用于自交叉多边形。<BR><BR> <DIV class=msgheader> <DIV class=right><a href="http://bbs.esrichina-bj.cn/ESRI/viewthread.php?tid=8771;extra=page%3D6###" target="_blank" >[Copy to clipboard]</A></DIV>CODE:</DIV> <DIV>// calculate Area<BR>function calcArea(PointX,PointY,MapUnits) {<BR><BR> var Count = PointX.length<BR> if (Count>3) {<BR> var mtotalArea = 0;<BR> if((PointX[0]!=PointX[Count-1])||(PointY[0]!=PointY[Count-1]))<BR> {<BR> return;<BR> }<BR> <BR> if (MapUnits=="DEGREES")//经纬度坐标下的球面多边形<BR> {<BR> var LowX=0.0;<BR> var LowY=0.0;<BR> var MiddleX=0.0;<BR> var MiddleY=0.0; <BR> var HighX=0.0;<BR> var HighY=0.0;<BR><BR> var AM = 0.0; <BR> var BM = 0.0; <BR> var CM = 0.0; <BR><BR> var AL = 0.0; <BR> var BL = 0.0; <BR> var CL = 0.0; <BR> <BR> var AH = 0.0; <BR> var BH = 0.0; <BR> var CH = 0.0; <BR><BR> var CoefficientL = 0.0;<BR> var CoefficientH = 0.0; <BR> <BR> var ALtangent = 0.0; <BR> var BLtangent = 0.0; <BR> var CLtangent = 0.0; <BR><BR> var AHtangent = 0.0; <BR> var BHtangent = 0.0; <BR> var CHtangent = 0.0;<BR> <BR> var ANormalLine = 0.0; <BR> var BNormalLine = 0.0; <BR> var CNormalLine = 0.0;<BR> <BR> var OrientationValue = 0.0; <BR> <BR> var AngleCos = 0.0;<BR><BR> var Sum1 = 0.0; <BR> var Sum2 = 0.0; <BR> var Count2 = 0; <BR> var Count1 = 0; <BR> <BR> <BR> var Sum = 0.0;<BR> var Radius = 6378000; <BR> <BR> for(i=0;i<Count;i++)<BR> {<BR> if(i==0)<BR> {<BR> LowX = PointX[Count-1] * Math.PI / 180;<BR> LowY = PointY[Count-1] * Math.PI / 180; <BR> MiddleX = PointX[0] * Math.PI / 180;<BR> MiddleY = PointY[0] * Math.PI / 180;<BR> HighX = PointX[1] * Math.PI / 180;<BR> HighY = PointY[1] * Math.PI / 180;<BR> }<BR> else if(i==Count-1)<BR> {<BR> LowX = PointX[Count-2] * Math.PI / 180;<BR> LowY = PointY[Count-2] * Math.PI / 180; <BR> MiddleX = PointX[Count-1] * Math.PI / 180;<BR> MiddleY = PointY[Count-1] * Math.PI / 180; <BR> HighX = PointX[0] * Math.PI / 180;<BR> HighY = PointY[0] * Math.PI / 180; <BR> }<BR> else<BR> {<BR> LowX = PointX[i-1] * Math.PI / 180;<BR> LowY = PointY[i-1] * Math.PI / 180; <BR> MiddleX = PointX * Math.PI / 180;<BR> MiddleY = PointY * Math.PI / 180; <BR> HighX = PointX[i+1] * Math.PI / 180;<BR> HighY = PointY[i+1] * Math.PI / 180; <BR> }<BR> <BR> AM = Math.cos(MiddleY) * Math.cos(MiddleX);<BR> BM = Math.cos(MiddleY) * Math.sin(MiddleX);<BR> CM = Math.sin(MiddleY);<BR> AL = Math.cos(LowY) * Math.cos(LowX);<BR> BL = Math.cos(LowY) * Math.sin(LowX);<BR> CL = Math.sin(LowY);<BR> AH = Math.cos(HighY) * Math.cos(HighX);<BR> BH = Math.cos(HighY) * Math.sin(HighX);<BR> CH = Math.sin(HighY); <BR> <BR> <BR> CoefficientL = (AM*AM + BM*BM + CM*CM)/(AM*AL + BM*BL + CM*CL);<BR> CoefficientH = (AM*AM + BM*BM + CM*CM)/(AM*AH + BM*BH + CM*CH);<BR> <BR> ALtangent = CoefficientL * AL - AM;<BR> BLtangent = CoefficientL * BL - BM;<BR> CLtangent = CoefficientL * CL - CM;<BR> AHtangent = CoefficientH * AH - AM;<BR> BHtangent = CoefficientH * BH - BM;<BR> CHtangent = CoefficientH * CH - CM; <BR> <BR> <BR> AngleCos = (AHtangent * ALtangent + BHtangent * BLtangent + CHtangent * CLtangent)/(Math.sqrt(AHtangent * AHtangent + BHtangent * BHtangent +CHtangent * CHtangent) * Math.sqrt(ALtangent * ALtangent + BLtangent * BLtangent +CLtangent * CLtangent));<BR> <BR> AngleCos = Math.acos(AngleCos);<BR> <BR> ANormalLine = BHtangent * CLtangent - CHtangent * BLtangent;<BR> BNormalLine = 0 - (AHtangent * CLtangent - CHtangent * ALtangent); <BR> CNormalLine = AHtangent * BLtangent - BHtangent * ALtangent;<BR> <BR> if(AM!=0) <BR> OrientationValue = ANormalLine/AM;<BR> else if(BM!=0) <BR> OrientationValue = BNormalLine/BM;<BR> else<BR> OrientationValue = CNormalLine/CM;<BR> <BR> if(OrientationValue>0)<BR> {<BR> Sum1 += AngleCos;<BR> Count1 ++;<BR> <BR> }<BR> else<BR> {<BR> Sum2 += AngleCos;<BR> Count2 ++;<BR> //Sum +=2*Math.PI-AngleCos;<BR> }<BR><BR> }<BR> <BR> if(Sum1>Sum2){<BR> Sum = Sum1+(2*Math.PI*Count2-Sum2);<BR> }<BR> else{<BR> Sum = (2*Math.PI*Count1-Sum1)+Sum2;<BR> }<BR> <BR> //平方米<BR> mtotalArea = (Sum-(Count-2)*Math.PI)*Radius*Radius;<BR> }<BR> else { //非经纬度坐标下的平面多边形<BR><BR> var i,j;<BR> var j;<BR> var p1x,p1y;<BR> var p2x,p2y;<BR> for(i=Count-1, j=0; j<Count; i=j, j++)<BR> {<BR> if(i==Count-1)<BR> {<BR> p1x = mX;<BR> p1y = mY; <BR> }<BR> else<BR> {<BR> p1x = PointX;<BR> p1y = PointY; <BR> }<BR> if(j==Count-1)<BR> {<BR> p2x = mX;<BR> p2y = mY; <BR> }<BR> else<BR> {<BR> p2x = PointX[j];<BR> p2y = PointY[j]; <BR> } <BR> mtotalArea +=p1x*p2y-p2x*p1y;<BR> }<BR> mtotalArea /= 2.0;<BR> }<BR> return mtotalArea;<BR> }<BR> return;<BR>}</DIV><BR>不太好注释，具体原理请参考古人的定理：<BR>球面多边形计算面积的关键在于计算多边形所有角的度数.<BR>对于球面n边形,所有角的和为S,球的半径为R,那么其面积就是<BR>R^2*(S-(n-2)*Pi)<BR></DIV> |
|
|
1楼#
发布于：2007-05-21 13:06
<P>谢谢共享，不知道精度如何</P>
<P>球面的是用高斯投影的算法吗</P> |
|
|
|
2楼#
发布于：2007-07-16 16:24
mX,mY指的是什么啊!
|
|
|
|
3楼#
发布于：2007-07-16 16:37
<img src="images/post/smile/dvbbs/em01.gif" />
|
|
|
|
4楼#
发布于：2007-07-16 19:21
<P>不胜感激</P>
|
|
|
5楼#
发布于：2008-02-26 09:47
路过.不错的同志.谢谢你,虽然我暂时用不到.
|
|
|
6楼#
发布于：2009-03-04 17:32
有关原理这块能不能说明一下，另外多边形的边是通过球心的大弧，还是线性插值的弧？
|
|