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