方位角和高度角

go.mod

@@ -50,6 +50,7 @@ require (
 	github.com/spf13/afero v1.11.0 // indirect
 	github.com/spf13/cast v1.6.0 // indirect
 	github.com/spf13/pflag v1.0.5 // indirect
+	github.com/starainrt/astro v0.0.0-20240209133137-7fd9deb71470 // indirect
 	github.com/subosito/gotenv v1.6.0 // indirect
 	go.mongodb.org/mongo-driver v1.11.4 // indirect
 	go.opentelemetry.io/contrib/instrumentation/net/http/otelhttp v0.49.0 // indirect

go.sum

@@ -1226,6 +1226,8 @@ github.com/spf13/pflag v1.0.5/go.mod h1:McXfInJRrz4CZXVZOBLb0bTZqETkiAhM9Iw0y3An
 github.com/spf13/viper v1.8.1/go.mod h1:o0Pch8wJ9BVSWGQMbra6iw0oQ5oktSIBaujf1rJH9Ns=
 github.com/spf13/viper v1.18.2 h1:LUXCnvUvSM6FXAsj6nnfc8Q2tp1dIgUfY9Kc8GsSOiQ=
 github.com/spf13/viper v1.18.2/go.mod h1:EKmWIqdnk5lOcmR72yw6hS+8OPYcwD0jteitLMVB+yk=
+github.com/starainrt/astro v0.0.0-20240209133137-7fd9deb71470 h1:agRZxDh63dxmxhSRdUHHWVH11MjrdbQq0XJy0dYYrDs=
+github.com/starainrt/astro v0.0.0-20240209133137-7fd9deb71470/go.mod h1:TN9w1RvRiEuEX3tgaMmw/3nrd7EDjzcP3wYfbK7oRgc=
 github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME=
 github.com/stretchr/objx v0.4.0/go.mod h1:YvHI0jy2hoMjB+UWwv71VJQ9isScKT/TqJzVSSt89Yw=
 github.com/stretchr/objx v0.5.0/go.mod h1:Yh+to48EsGEfYuaHDzXPcE3xhTkx73EhmCGUpEOglKo=

pkg/producer/angle.go

@@ -0,0 +1,42 @@
+package producer
+
+import (
+	"math"
+	"time"
+
+	"github.com/starainrt/astro/sun"
+)
+
+// Azimuth 太阳方位角
+// 返回给定经纬度、对应date时区date时刻的太阳方位角（正北为0，向东增加）
+func SunAzimuth(date time.Time, lon, lat float64) float64 {
+	return sun.Azimuth(date, lon, lat)
+}
+
+// Zenith 太阳高度角
+// 返回给定经纬度、对应date时区date时刻的太阳高度角
+func SunZenith(date time.Time, lon, lat float64) float64 {
+	return sun.Zenith(date, lon, lat)
+}
+
+func SatAzimuth(lon, lat, lonSat, latSat float64) float64 {
+	lon = lon * math.Pi / 180
+	lat = lat * math.Pi / 180
+	lonSat = lonSat * math.Pi / 180
+	latSat = latSat * math.Pi / 180
+
+	beta := math.Acos(math.Cos(lat-latSat) * math.Cos(lon-lonSat))
+	phi := math.Asin(math.Sin(lonSat-lon) / math.Sin(beta))
+	return phi * 180 / math.Pi
+}
+
+func SatZenith(lon, lat, lonSat, latSat float64) float64 {
+	lon = lon * math.Pi / 180
+	lat = lat * math.Pi / 180
+	lonSat = lonSat * math.Pi / 180
+	latSat = latSat * math.Pi / 180
+
+	beta := math.Acos(math.Cos(lat-latSat) * math.Cos(lon-lonSat))
+	theta := math.Asin(42164.0 * math.Sin(beta) / math.Sqrt(1.8084*1e9-5.3725*1e8*math.Cos(beta)))
+	return theta * 180 / math.Pi
+}

pkg/producer/aux.go

@@ -140,6 +140,25 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P
 	scene.Meta.Pitch = aux0.Eular2 * 180 / math.Pi
 	scene.Meta.Roll = aux0.Eular1 * 180 / math.Pi
 
+	_, centerT, _ := r.SceneImageTime(scene)
+	scene.Meta.SunAzimuth = SunAzimuth(centerT,
+		scene.Meta.CentreLocation.Longitude,
+		scene.Meta.CentreLocation.Latitude)
+	scene.Meta.SunElevation = SunZenith(centerT,
+		scene.Meta.CentreLocation.Longitude,
+		scene.Meta.CentreLocation.Latitude)
+
+	scene.Meta.SatAzimuth = SatAzimuth(
+		scene.Meta.CentreLocation.Longitude,
+		scene.Meta.CentreLocation.Latitude,
+		scene.Meta.SatPosY,
+		scene.Meta.SatPosX)
+	scene.Meta.SatElevation = SatZenith(
+		scene.Meta.CentreLocation.Longitude,
+		scene.Meta.CentreLocation.Latitude,
+		scene.Meta.SatPosY,
+		scene.Meta.SatPosX)
+
 	// 计算RPC
 	rpc := NewRPC(r, scene)
 	if err := rpc.RPC(); err != nil {

pkg/producer/meta.go

@@ -2,11 +2,12 @@ package producer
 
 import (
 	"encoding/xml"
-	"math"
 	"os"
 	"path/filepath"
 	"strings"
 	"time"
+
+	"github.com/duke-git/lancet/v2/mathutil"
 )
 
 // 定义与XML结构对应的Go结构体
@@ -60,22 +61,14 @@ type Corners struct {
 }
 
 func (corners Corners) Extend() (lng1, lat1, lng2, lat2 float64) {
-	lng1 = math.Min(corners.UpperLeft.Longitude, corners.LowerLeft.Longitude)
-	lng1 = math.Min(lng1, corners.UpperRight.Longitude)
-	lng1 = math.Min(lng1, corners.LowerRight.Longitude)
-
-	lat1 = math.Max(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude)
-	lat1 = math.Max(lat1, corners.UpperRight.Latitude)
-	lat1 = math.Max(lat1, corners.LowerRight.Latitude)
-
-	lng2 = math.Max(corners.UpperLeft.Longitude, corners.LowerLeft.Longitude)
-	lng2 = math.Max(lng2, corners.UpperRight.Longitude)
-	lng2 = math.Max(lng2, corners.LowerRight.Longitude)
-
-	lat2 = math.Min(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude)
-	lat2 = math.Min(lat2, corners.UpperRight.Latitude)
-	lat2 = math.Min(lat2, corners.LowerRight.Latitude)
-
+	lng1 = mathutil.Min(corners.UpperLeft.Longitude, corners.LowerLeft.Longitude,
+		corners.UpperRight.Longitude, corners.LowerRight.Longitude)
+	lng2 = mathutil.Max(corners.UpperLeft.Longitude, corners.LowerLeft.Longitude,
+		corners.UpperRight.Longitude, corners.LowerRight.Longitude)
+	lat1 = mathutil.Max(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude,
+		corners.UpperRight.Latitude, corners.LowerRight.Latitude)
+	lat2 = mathutil.Min(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude,
+		corners.UpperRight.Latitude, corners.LowerRight.Latitude)
 	return
 }