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

// calculate Area
public static double calcArea(ArrayList PointX, ArrayList PointY, string MapUnits)
{

int Count = PointX.Count;
if (Count > 2)
{
double mtotalArea = 0;

if (MapUnits == "DEGREES")//经纬度坐标下的球面多边形
{
double LowX = 0.0;
double LowY = 0.0;
double MiddleX = 0.0;
double MiddleY = 0.0;
double HighX = 0.0;
double HighY = 0.0;

double AM = 0.0;
double BM = 0.0;
double CM = 0.0;

double AL = 0.0;
double BL = 0.0;
double CL = 0.0;

double AH = 0.0;
double BH = 0.0;
double CH = 0.0;

double CoefficientL = 0.0;
double CoefficientH = 0.0;

double ALtangent = 0.0;
double BLtangent = 0.0;
double CLtangent = 0.0;

double AHtangent = 0.0;
double BHtangent = 0.0;
double CHtangent = 0.0;

double ANormalLine = 0.0;
double BNormalLine = 0.0;
double CNormalLine = 0.0;

double OrientationValue = 0.0;

double AngleCos = 0.0;

double Sum1 = 0.0;
double Sum2 = 0.0;
double Count2 = 0;
double Count1 = 0;

double Sum = 0.0;
double Radius = 6378000;

for (int i = 0; i < Count; i++)
{
if (i == 0)
{
LowX =(double) PointX[Count - 1] * Math.PI / 180;
LowY = (double)PointY[Count - 1] * Math.PI / 180;
MiddleX = (double)PointX[0] * Math.PI / 180;
MiddleY = (double)PointY[0] * Math.PI / 180;
HighX = (double)PointX[1] * Math.PI / 180;
HighY = (double)PointY[1] * Math.PI / 180;
}
else if (i == Count - 1)
{
LowX = (double)PointX[Count - 2] * Math.PI / 180;
LowY = (double)PointY[Count - 2] * Math.PI / 180;
MiddleX = (double)PointX[Count - 1] * Math.PI / 180;
MiddleY = (double)PointY[Count - 1] * Math.PI / 180;
HighX = (double)PointX[0] * Math.PI / 180;
HighY = (double)PointY[0] * Math.PI / 180;
}
else
{
LowX = (double)PointX[i - 1] * Math.PI / 180;
LowY = (double)PointY[i - 1] * Math.PI / 180;
MiddleX = (double)PointX[i] * Math.PI / 180;
MiddleY = (double)PointY[i] * Math.PI / 180;
HighX = (double)PointX[i + 1] * Math.PI / 180;
HighY = (double)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
{ //非经纬度坐标下的平面多边形

int i, j;
//double j;
double p1x, p1y;
double p2x, p2y;
for ( i = Count - 1, j = 0; j < Count; i = j, j++)
{

p1x = (double)PointX[i];
p1y = (double)PointY[i];

p2x = (double)PointX[j];
p2y = (double)PointY[j];

mtotalArea += p1x * p2y - p2x * p1y;
}
mtotalArea /= 2.0;
}
return mtotalArea;
}
return 0;
}

from：http://blog.csdn.net/lockepeak/article/details/6079206