This commit is contained in:
nuknal
2024-05-20 11:08:24 +08:00
parent f609b2b023
commit efd824cc8a
14 changed files with 332 additions and 69 deletions

10
extract/attitude.go Normal file
View File

@@ -0,0 +1,10 @@
package extract
/*
将欧拉角 (yaw, pitch, roll) 转换为旋转矩阵
参数:
- yaw: 偏航角绕Z轴旋转
- pitch: 俯仰角绕Y轴旋转
- roll: 滚转角绕X轴旋转
*/

View File

@@ -142,6 +142,9 @@ func (e *Extractor) SeprateAuxAndImgData(dataFile string, segmentIndex int) erro
var afh AuxFrameHead
msdata := make([][]byte, 4)
var header []byte
var ebAux []byte
var platAux []byte
for i := 0; i < dataLen; {
if i+4 > dataLen {
@@ -182,11 +185,17 @@ func (e *Extractor) SeprateAuxAndImgData(dataFile string, segmentIndex int) erro
break
}
auxTime := binary.BigEndian.Uint32(data[i+32 : i+36])
if math.Abs(float64(auxTime)-float64(afh.TimeSec)) < 30 {
header = append(header, data[i:i+24]...)
}
// 存储辅助数据到临时文件
dataIndex := i + 24
lw.ws[EB_AUX].w.Write(data[dataIndex : dataIndex+8]) // 8字节焦面电箱辅助数据
// lw.ws[EB_AUX].w.Write(data[dataIndex : dataIndex+8]) // 8字节焦面电箱辅助数据
ebAux = append(ebAux, data[dataIndex:dataIndex+8]...)
dataIndex += 8
lw.ws[PLAT_AUX].w.Write(data[dataIndex : dataIndex+32]) // 32字节中心机辅助数据
// lw.ws[PLAT_AUX].w.Write(data[dataIndex : dataIndex+32]) // 32字节中心机辅助数据
platAux = append(platAux, data[dataIndex:dataIndex+32]...)
dataIndex += 32
// 存储图像数据到临时文件 - 以 ENVI BSQ 格式存储,同时提供 HDR 描述文件
@@ -221,15 +230,41 @@ func (e *Extractor) SeprateAuxAndImgData(dataFile string, segmentIndex int) erro
i = dataIndex // 完成一行数据解析
}
var bands = 0
for i := 0; i < 4; i++ {
if len(msdata[i]) > 0 {
log.Println("write mss data of band B", i+1)
_, err := write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[i])
if err != nil {
log.Error("write mss data error:", err.Error())
}
bands += 1
// var bands = 0
// for i := 0; i < 4; i++ {
// if len(msdata[i]) > 0 {
// log.Println("write mss data of band B", i+1)
// _, err := write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[i])
// if err != nil {
// log.Error("write mss data error:", err.Error())
// }
// bands += 1
// }
// }
if len(msdata[0]) != len(msdata[1]) || len(msdata[2]) != len(msdata[3]) {
log.Error("mss data of bands B1-B4 are not equal")
return errors.New("mss data of bands B1-B4 are not equal")
}
mssRowLen := 1192 * 4
for i := 0; i < len(msdata[0]); i += mssRowLen {
var err error
_, err = write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[0][i:i+mssRowLen])
if err != nil {
log.Error("write mss 1 data error:", err.Error())
}
_, err = write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[1][i:i+mssRowLen])
if err != nil {
log.Error("write mss 2 data error:", err.Error())
}
_, err = write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[2][i:i+mssRowLen])
if err != nil {
log.Error("write mss 3 data error:", err.Error())
}
_, err = write16bPixelLittleEndian(lw.ws[MSS_RAW].w, msdata[3][i:i+mssRowLen])
if err != nil {
log.Error("write mss 4 data error:", err.Error())
}
}
@@ -238,16 +273,29 @@ func (e *Extractor) SeprateAuxAndImgData(dataFile string, segmentIndex int) erro
lw.ws[PAN_HDR].w.Write([]byte(panEnviHdr.String()))
mssEnviHdr.Lines = mssEnviHdr.Lines / 4 // 多光谱波段分别在 4 行中传输
mssEnviHdr.Samples = 2384
mssEnviHdr.Bands = bands
mssEnviHdr.Samples = 2384 * 4
mssEnviHdr.Bands = 1
lw.ws[MSS_HDR].w.Write([]byte(mssEnviHdr.String()))
// 帧头+辅助数据
if len(header)/24 < len(ebAux)/128 || len(ebAux)/128 != len(platAux)/512 {
fmt.Println("aux data length:", len(header)/24, len(ebAux)/128, len(platAux)/512)
return errors.New("aux data length is not equal")
}
auxRows := len(header) / 24
for i := 0; i < auxRows; i++ {
lw.ws[AUX].w.Write(header[i*24 : (i+1)*24])
lw.ws[AUX].w.Write(ebAux[i*128 : (i+1)*128])
lw.ws[AUX].w.Write(platAux[i*512 : (i+1)*512])
}
return nil
}
func (e *Extractor) quanlityAnalysis(data []byte) {
// 数据质量分析
fimg, _ := os.Create("demo/temp/" + e.params.DataId + "_aux_img.txt")
fimg, _ := os.Create("demo/temp/" + e.params.DataId + "/aux.txt")
defer fimg.Close()
fimg.WriteString("字节数 帧头流水号 文件号 帧头时间 中心辅助数据前4字节(行33-36)\n")

View File

@@ -346,14 +346,22 @@ func (e *Extractor) ParseAuxPlatform(auxfile string) ([]*AuxPlatform, error) {
}
var aps []*AuxPlatform
for i := 0; i < len(data); i += AuxPlatformFrmSize {
for i := 0; i < len(data); i += 24 + 128 + 512 {
ap := AuxPlatform{}
if err := ap.Parse(data[i : i+AuxPlatformFrmSize]); err != nil {
if err := ap.Parse(data[i+24+128 : i+24+128+512]); err != nil {
return nil, err
}
aps = append(aps, &ap)
ap.Print()
// Calculate(
// Vector3{ap.J2000Pos_X, ap.J2000Pos_Y, ap.J2000Pos_Y},
// Vector3{ap.J2000Vel_X, ap.J2000Vel_Y, ap.J2000Vel_Z},
// Quaternion{ap.SS1_Q1, ap.SS1_Q2, ap.SS1_Q3, ap.SS1_Q4},
// )
// pp.Println(calculator.WGS84XYZtoLngLat(ap.W84PosX, ap.W84PosY, ap.W84PosZ))
// pp.Println(calculator.WGS84XYZtoLngLat(ap.WGS84PosX, ap.WGS84PosY, ap.WGS84PosZ))
break

114
extract/calculator.go Normal file
View File

@@ -0,0 +1,114 @@
package extract
import (
"fmt"
"math"
)
// 定义常数
const (
GM = 3.986004418e14 // 地球引力常数m^3/s^2
Re = 6378137.0 // 地球半径m
)
// 四元数
type Quaternion struct {
w, x, y, z float64
}
// 矢量
type Vector3 struct {
x, y, z float64
}
// 计算轨道参数
func calculateOrbitalParameters(position, velocity Vector3) (float64, float64, float64, float64, float64) {
r := math.Sqrt(position.x*position.x + position.y*position.y + position.z*position.z)
v := math.Sqrt(velocity.x*velocity.x + velocity.y*velocity.y + velocity.z*velocity.z)
// 计算轨道参数
energy := 0.5*v*v - GM/r // 比能
a := -GM / (2 * energy) // 轨道半长轴
e := math.Sqrt(1 - (math.Pow(math.Sqrt(math.Pow(position.x*velocity.y-position.y*velocity.x, 2)+
math.Pow(position.y*velocity.z-position.z*velocity.y, 2)+
math.Pow(position.z*velocity.x-position.x*velocity.z, 2)), 2) / (GM * a)))
i := math.Acos(position.z / r) // 轨道倾角
Omega := math.Atan2(position.x, -position.y) // 升交点赤经
omega := math.Atan2(position.z*velocity.x-position.x*velocity.z, position.x*velocity.y-position.y*velocity.x) // 升交点赤纬
return a, e, i * 180 / math.Pi, Omega * 180 / math.Pi, omega * 180 / math.Pi
}
// 四元数到旋转矩阵
func quaternionToRotationMatrix(q Quaternion) [3][3]float64 {
w, x, y, z := q.w, q.x, q.y, q.z
return [3][3]float64{
{1 - 2*y*y - 2*z*z, 2*x*y - 2*z*w, 2*x*z + 2*y*w},
{2*x*y + 2*z*w, 1 - 2*x*x - 2*z*z, 2*y*z - 2*x*w},
{2*x*z - 2*y*w, 2*y*z + 2*x*w, 1 - 2*x*x - 2*y*y},
}
}
// 向量加法
func add(v1, v2 Vector3) Vector3 {
return Vector3{v1.x + v2.x, v1.y + v2.y, v1.z + v2.z}
}
// 向量数乘
func scale(v Vector3, c float64) Vector3 {
return Vector3{v.x * c, v.y * c, v.z * c}
}
// 矩阵乘法
func matrixMult(m [3][3]float64, v Vector3) Vector3 {
return Vector3{
m[0][0]*v.x + m[0][1]*v.y + m[0][2]*v.z,
m[1][0]*v.x + m[1][1]*v.y + m[1][2]*v.z,
m[2][0]*v.x + m[2][1]*v.y + m[2][2]*v.z,
}
}
// 计算图像位置
func calculateImagePosition(roll, pitch, yaw float64, quat Quaternion, satellitePos, satelliteVel Vector3) (float64, float64) {
// 将角度转换为弧度
roll = roll * math.Pi / 180
pitch = pitch * math.Pi / 180
yaw = yaw * math.Pi / 180
// 构造旋转矩阵
R := quaternionToRotationMatrix(quat)
// 传感器坐标系中的Z轴向量
sensorZ := Vector3{0, 0, 1}
// 将Z轴向量通过旋转矩阵转换到卫星坐标系中
lineOfSight := matrixMult(R, sensorZ)
// 计算交点经纬度
lat := math.Asin(lineOfSight.z) * 180 / math.Pi
lon := math.Atan2(lineOfSight.y, lineOfSight.x) * 180 / math.Pi
return lat, lon
}
func Calculate(satellitePos, satelliteVel Vector3, quat Quaternion) {
// 假设一些已知数据
// satellitePos := Vector3{x: 7000e3, y: 0, z: 0} // 假设一个简单的地心坐标位置
// satelliteVel := Vector3{x: 0, y: 7.5e3, z: 0} // 假设一个简单的速度
// quat := Quaternion{w: 0.7071, x: 0.7071, y: 0, z: 0} // 假设一个简单的四元数
roll := 0.0
pitch := 0.0
yaw := 0.0
// 计算轨道参数
a, e, i, Omega, omega := calculateOrbitalParameters(satellitePos, satelliteVel)
// 打印轨道参数
fmt.Println("轨道参数:")
fmt.Printf("轨道半长轴 (a): %.2f m\n", a)
fmt.Printf("偏心率 (e): %.6f\n", e)
fmt.Printf("轨道倾角 (i): %.2f 度\n", i)
fmt.Printf("升交点赤经 (Ω): %.2f 度\n", Omega)
fmt.Printf("近地点幅角 (ω): %.2f 度\n", omega)
// 计算图像位置
lat, lon := calculateImagePosition(roll, pitch, yaw, quat, satellitePos, satelliteVel)
fmt.Printf("\n图像位置\n纬度: %.6f\n经度: %.6f\n", lat, lon)
}

1
extract/orbit.go Normal file
View File

@@ -0,0 +1 @@
package extract

51
extract/util.go Normal file
View File

@@ -0,0 +1,51 @@
package extract
import (
"bytes"
"encoding/binary"
)
// int32ToBigEndian converts an int32 to a big-endian byte array.
func int32ToBigEndian(num int32) ([]byte, error) {
buf := new(bytes.Buffer)
err := binary.Write(buf, binary.BigEndian, num)
if err != nil {
return nil, err
}
return buf.Bytes(), nil
}
// bigEndianToInt32 converts a big-endian byte array to an int32.
func bigEndianToInt32(data []byte) int32 {
var num int32
buf := bytes.NewReader(data)
err := binary.Read(buf, binary.BigEndian, &num)
if err != nil {
return 0
}
return num
}
// bigEndianToInt24 converts a big-endian 3-byte array to an int32 (simulating int24).
func bigEndianToInt24(data []byte) int32 {
if len(data) != 3 {
return 0
}
num := int32(data[0])<<16 | int32(data[1])<<8 | int32(data[2])
// Adjust for negative values
if num&0x800000 != 0 {
num |= ^0xFFFFFF
}
return num
}
// bigEndianToInt16 converts a big-endian byte array to an int16.
func bigEndianToInt16(data []byte) int16 {
var num int16
buf := bytes.NewReader(data)
err := binary.Read(buf, binary.BigEndian, &num)
if err != nil {
return 0
}
return num
}

View File

@@ -5,13 +5,15 @@ import (
"io"
"os"
"path/filepath"
"strings"
"github.com/sirupsen/logrus"
)
const (
EB_AUX = "EB.AUX"
PLAT_AUX = "PLAT.AUX"
AUX = ".AUX"
EB_AUX = ".EB.AUX"
PLAT_AUX = ".PLAT.AUX"
PAN_RAW = "PAN.RAW"
MSS_RAW = "MSS.RAW"
PAN_HDR = "PAN.HDR"
@@ -32,8 +34,8 @@ type L0Writer struct {
func newL0Writer(outputDir, name string) *L0Writer {
var err error
suffix := []string{EB_AUX, PLAT_AUX, PAN_RAW, MSS_RAW, PAN_HDR, MSS_HDR}
suffix := []string{AUX, PAN_RAW, MSS_RAW, PAN_HDR, MSS_HDR}
nameEles := strings.Split(name, "_")
lw := &L0Writer{
ws: make(map[string]*l0w),
mssStoredType: MSS_ALL_IN_ONE,
@@ -41,7 +43,14 @@ func newL0Writer(outputDir, name string) *L0Writer {
for _, s := range suffix {
w := &l0w{}
w.name = filepath.Join(outputDir, name+"_"+s)
fname := name + s
if s == MSS_HDR || s == MSS_RAW || s == PAN_HDR || s == PAN_RAW {
eles := strings.Split(s, ".")
nameEles[1] = eles[0]
fname = strings.Join(nameEles, "_") + "." + eles[1]
}
w.name = filepath.Join(outputDir, fname)
w.f, err = os.OpenFile(w.name, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
if err != nil {
logrus.Panic("create file failed: ", w.name, err)