// 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;
}