门店配送面积计算(根据经纬度计算多边形球面积-凸多边形和凹多边形)

business/jxutils/jxutils.go

@@ -192,6 +192,175 @@ func EarthDistance(lng1, lat1, lng2, lat2 float64) float64 {
 	return dist * radius
 }
 
+//计算多边形平面面积
+func CalcPolygonArea(points [][2]float64) (area float64) {
+	count := len(points)
+	for i := 0; i < count-1; i++ {
+		area += (points[i][0] - points[i+1][0]) * (points[i][1] + points[i+1][1])
+	}
+	area += (points[count-1][0] - points[0][0]) * (points[count-1][1] + points[0][1])
+	area = math.Abs(area) / 2
+
+	return area
+}
+
+func ConvertToRadian(value float64) float64 {
+	return value * math.Pi / 180
+}
+
+//计算凹多边形球面积根据经纬度(KM)
+func CalcConcavePolygonAreaByLngAndLatFix(points [][2]float64) (area float64) {
+	radius := 6378.137
+	for i := 0; i < len(points)-1; i++ {
+		p1 := points[i]
+		p2 := points[i+1]
+		area += ConvertToRadian(p2[0]-p1[0]) * (2 + math.Sin(ConvertToRadian(p1[1])) + math.Sin(ConvertToRadian(p2[1])))
+	}
+
+	area = area * radius * radius / 2
+
+	return math.Abs(area)
+}
+
+//计算凸多边形球面积根据经纬度(KM)
+func CalcConvexPolygonAreaByLngAndLat(points [][2]float64) (area float64) {
+	Count := len(points)
+	LowX := 0.0
+	LowY := 0.0
+	MiddleX := 0.0
+	MiddleY := 0.0
+	HighX := 0.0
+	HighY := 0.0
+
+	AM := 0.0
+	BM := 0.0
+	CM := 0.0
+
+	AL := 0.0
+	BL := 0.0
+	CL := 0.0
+
+	AH := 0.0
+	BH := 0.0
+	CH := 0.0
+
+	CoefficientL := 0.0
+	CoefficientH := 0.0
+
+	ALtangent := 0.0
+	BLtangent := 0.0
+	CLtangent := 0.0
+
+	AHtangent := 0.0
+	BHtangent := 0.0
+	CHtangent := 0.0
+
+	ANormalLine := 0.0
+	BNormalLine := 0.0
+	CNormalLine := 0.0
+
+	OrientationValue := 0.0
+
+	AngleCos := 0.0
+
+	Sum1 := 0.0
+	Sum2 := 0.0
+	Count2 := 0
+	Count1 := 0
+
+	Sum := 0.0
+	Radius := 6378.137
+
+	for i := 0; i < Count; i++ {
+		if i == 0 {
+			LowX = ConvertToRadian(points[Count-1][0])
+			LowY = ConvertToRadian(points[Count-1][1])
+			MiddleX = ConvertToRadian(points[0][0])
+			MiddleY = ConvertToRadian(points[0][1])
+			HighX = ConvertToRadian(points[1][0])
+			HighY = ConvertToRadian(points[1][1])
+		} else if i == Count-1 {
+			LowX = ConvertToRadian(points[Count-2][0])
+			LowY = ConvertToRadian(points[Count-2][1])
+			MiddleX = ConvertToRadian(points[Count-1][0])
+			MiddleY = ConvertToRadian(points[Count-1][1])
+			HighX = ConvertToRadian(points[0][0])
+			HighY = ConvertToRadian(points[0][1])
+		} else {
+			LowX = ConvertToRadian(points[i-1][0])
+			LowY = ConvertToRadian(points[i-1][1])
+			MiddleX = ConvertToRadian(points[i][0])
+			MiddleY = ConvertToRadian(points[i][1])
+			HighX = ConvertToRadian(points[i+1][0])
+			HighY = ConvertToRadian(points[i+1][1])
+		}
+
+		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*float64(Count2) - Sum2
+	} else {
+		Sum = 2*math.Pi*float64(Count1) - Sum1 + Sum2
+	}
+
+	area = (Sum - float64(Count-2)*math.Pi) * Radius * Radius
+
+	return area
+}
+
+//计算多边形球面积根据经纬度(KM)
+func CalcPolygonAreaByLngAndLat(points [][2]float64) (area float64) {
+	area = CalcConvexPolygonAreaByLngAndLat(points)
+	if math.IsNaN(area) {
+		area = CalcConcavePolygonAreaByLngAndLatFix(points)
+	}
+
+	return area
+}
+
 func StandardCoordinate2Int(value float64) int {
 	return int(math.Round(value * 1000000))
 }