关于经纬度得到的多边形面积。

（两个点）半正矢公式计算（Haversine formula）：

关于经纬度得到的多边形面积。

public static class AreaOfTwoPoint

{

private static double EARTH_RADIUS = 6370996.81;

public static double AreaOfPoint(List<Point> list)

{

double result;

if (list.Count<Point>() == 2)

{

double lat = list[0].lat;

double lat2 = list[0].lat;

double lon = list[1].lon;

double lon2 = list[1].lon;

double num = AreaOfTwoPoint.Distance(lat, lon , lat2, lon2);

result = 3.1415926535897931 * Math.Pow(num / 2.0, 2.0);

}

else

{

result = 0.0;

}

return result;

}

private static double HvaerSin(double thera)

{

double num = Math.Sin(thera / 2.0);

return num * num;

}

private static double Distance(double lat1, double lon1, double lat2, double lon2)

{

lat1 = AreaOfTwoPoint.ConvertDegreesToRadians(lat1);

lon1 = AreaOfTwoPoint.ConvertDegreesToRadians(lon1);

lat2 = AreaOfTwoPoint.ConvertDegreesToRadians(lat2);

lon2 = AreaOfTwoPoint.ConvertDegreesToRadians(lon2);

double thera = Math.Abs(lon1 - lon2);

double thera2 = Math.Abs(lat1 - lat2);

double d = AreaOfTwoPoint.HvaerSin(thera2) + Math.Cos(lat1) * Math.Cos(lat2) * AreaOfTwoPoint.HvaerSin(thera);

return 2.0 * AreaOfTwoPoint.EARTH_RADIUS * Math.Asin(Math.Sqrt(d));

}

private static double ConvertDegreesToRadians(double degrees)

{

return degrees * 3.1415926535897931 / 180.0;

}

　　具体请参考：https://en.wikipedia.org/wiki/Haversine_formula

　　多个点，宝宝没本事，没有办法，只能就说按照多边形的投影走了，就是得到多边形，然后切割成N-2个小三角。然后计算面积。

　　在这里必须要提醒的是，在处理的时候，有时候当三个点在经纬度上面显示出来约等于一条直线的时候，因为小数问题，可能会被计算机识别为一条直线，那么，此时，得到的面积就会无穷大，所以，在实际运用过程中，请注意，每个点之间的对应的投影点，及时做好处理。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

private double CalcArea(List<Points> list)

{

var count = list.Count;

if (count > 2)

{

//数组中的元素值

double mtotalArea = 0;

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, BM = 0.0, CM = 0.0;

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

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

double CoefficientL = 0.0, CoefficientH = 0.0;

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

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

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

//定位置

double OrientationValue = 0.0;

//余弦函数

double AngleCos = 0.0;

double Sum1 = 0.0, Sum2 = 0.0;

double Count1 = 0,Count2 = 0;

double Sum = 0.0;

double Radius = 6370996.81;//地球半径

for (int i = 0; i < count; i++)

{

//坐标系中，一般X代表纬度(Lon)，Y代表经度(Lat)

if (i == 0)

{

LowX = (list[count - 1].Lon) * Math.PI / 180;

LowY = (list[count - 1].Lat) * Math.PI / 180;

MiddleX = (list[0].Lon) * Math.PI / 180;

MiddleY = (list[0].Lat) * Math.PI / 180;

HighX = (list[1].Lon) * Math.PI / 180;

HighY = (list[1].Lat) * Math.PI / 180;

}

else if (i == count - 1)

{

LowX = (list[count-2].Lon) * Math.PI / 180;

LowY = (list[count-2].Lat) * Math.PI / 180;

MiddleX = (list[count - 1].Lon) * Math.PI / 180;

MiddleY = (list[count - 1].Lat) * Math.PI / 180;

HighX = (list[0].Lon) * Math.PI / 180;

HighY = (list[0].Lat) * Math.PI / 180;

}

else

{

LowX = (list[i-1].Lon) * Math.PI / 180;

LowY = (list[i - 1].Lat) * Math.PI / 180;

MiddleX = (list[i].Lon) * Math.PI / 180;

MiddleY = (list[i].Lat) * Math.PI / 180;

HighX = (list[i+1].Lon) * Math.PI / 180;

HighY = (list[i+1].Lat) * 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++;

}

if (Sum1 > Sum2)

{

Sum = Sum1 + (2 * Math.PI * Count2 - Sum2);

}

else

{

Sum = (2 * Math.PI * Count1 - Sum1) + Sum2;

}

return Math.Abs((Sum - (count - 2) * Math.PI) * Radius * Radius);

}

return 0;

}

关于经纬度得到的多边形面积。

相关推荐