From 8f2b297a02d8a01b6aadc2c95936bcc5faab2baf Mon Sep 17 00:00:00 2001 From: nuknal Date: Thu, 30 May 2024 18:11:42 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9A=82=E6=97=B6=E4=BD=BF=E7=94=A8=E6=98=9F?= =?UTF-8?q?=E4=B8=8B=E7=82=B9=E5=9D=90=E6=A0=87=E4=BD=9C=E4=B8=BA=E5=9B=BE?= =?UTF-8?q?=E5=83=8F=E5=B7=A6=E4=B8=8A=E8=A7=92=E5=9D=90=E6=A0=87?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cmd/fus.go | 38 +++ cmd/proc.go | 14 +- cmd/root.go | 1 + go.mod | 7 +- go.sum | 24 ++ pkg/auxilary/aux.go | 46 ++++ pkg/auxilary/aux_ebox.go | 134 +++++++++ pkg/auxilary/aux_elements.go | 420 +++++++++++++++++++++++++++++ pkg/auxilary/aux_head.go | 111 ++++++++ pkg/auxilary/aux_platform.go | 390 +++++++++++++++++++++++++++ pkg/auxilary/util.go | 51 ++++ pkg/calculator/const.go | 10 + pkg/calculator/j2000.go | 78 ++++++ pkg/calculator/quaternion.go | 70 +++++ pkg/calculator/wgs84.go | 77 ++++++ pkg/config/config.go | 4 + pkg/config/viper.go | 4 +- pkg/producer/aux.go | 109 ++++++++ pkg/producer/image_registration.go | 5 + pkg/producer/jpg.go | 32 +-- pkg/producer/meta.go | 36 ++- pkg/producer/output.go | 31 +-- pkg/producer/proj.go | 4 +- pkg/producer/scenes.go | 53 +++- pkg/producer/util.go | 45 +++- 25 files changed, 1710 insertions(+), 84 deletions(-) create mode 100644 cmd/fus.go create mode 100644 pkg/auxilary/aux.go create mode 100644 pkg/auxilary/aux_ebox.go create mode 100644 pkg/auxilary/aux_elements.go create mode 100644 pkg/auxilary/aux_head.go create mode 100644 pkg/auxilary/aux_platform.go create mode 100644 pkg/auxilary/util.go create mode 100644 pkg/calculator/const.go create mode 100644 pkg/calculator/j2000.go create mode 100644 pkg/calculator/quaternion.go create mode 100644 pkg/calculator/wgs84.go create mode 100644 pkg/producer/aux.go diff --git a/cmd/fus.go b/cmd/fus.go new file mode 100644 index 0000000..7828a7a --- /dev/null +++ b/cmd/fus.go @@ -0,0 +1,38 @@ +package main + +import ( + "path/filepath" + "strings" + + "github.com/sirupsen/logrus" + "github.com/spf13/cobra" + producer "starwiz.cn/sjy01/image-proc/pkg/producer" +) + +var ( + panImage string + mssImage string + outputDir string +) + +var fusCmd = &cobra.Command{ + Use: "fus", + Short: "use GDAL_Pansharpen to merge PAN and MSS images", + Run: func(cmd *cobra.Command, args []string) { + fusedTiff := filepath.Base(mssImage) + fusedTiff = strings.Replace(fusedTiff, "MSS", "FUS", -1) + err := producer.GDALPansharpen(panImage, mssImage, filepath.Join(outputDir, fusedTiff)) + if err != nil { + logrus.Fatal(err) + } + }, +} + +func init() { + rootCmd.AddCommand(fusCmd) + + fusCmd.Flags().StringVarP(&panImage, "pan", "p", "", "path to the PAN image") + fusCmd.Flags().StringVarP(&mssImage, "mss", "m", "", "path to the MSS image") + fusCmd.Flags().StringVarP(&outputDir, "output", "o", "", "path to the output directory") + fusCmd.Flags().StringVarP(¶msXML, "params", "x", "", "path to the XML file containing parameters for GDAL_Pansharpen") +} diff --git a/cmd/proc.go b/cmd/proc.go index 4a8687c..260f36e 100644 --- a/cmd/proc.go +++ b/cmd/proc.go @@ -22,7 +22,8 @@ var procCmd = &cobra.Command{ Short: "process images", Long: `process images`, Run: func(cmd *cobra.Command, args []string) { - _ = config.InitViper(configFile) + config.GViper = config.InitViper(configFile) + logrus.SetLevel(config.GCONFIG.LogLevel) reg := producer.NewRegistrator(producer.DownSampled) reg.Params = params @@ -30,19 +31,23 @@ var procCmd = &cobra.Command{ reg.Params.PanTiffFile = filepath.Join(params.OutputDir, strings.TrimSuffix(filepath.Base(params.PanRawFile), filepath.Ext(params.PanRawFile))+".tiff") reg.Params.FusTIffFile = strings.Replace(reg.Params.MssTiffFile, ".tiff", "_FUS.tiff", 1) + if err := reg.LoadAuxData(); err != nil { + logrus.Fatal(err) + } + if err := reg.LoadMssRaw(); err != nil { - panic(err) + logrus.Fatal(err) } if err := reg.LoadPanRaw(); err != nil { - panic(err) + logrus.Fatal(err) } godal.RegisterAll() os.MkdirAll(params.OutputDir, 0755) if err := reg.DoPhaseCorrelation(); err != nil { - panic(err) + logrus.Fatal(err) } reg.DoCoRegistration() @@ -85,4 +90,5 @@ func init() { procCmd.Flags().StringVarP(¶ms.OutputDir, "output-dir", "o", "data", "output directory") procCmd.Flags().BoolVarP(¶ms.SubScenes, "sub-scenes", "", false, "process sub-scenes") procCmd.Flags().BoolVarP(&saveStrip, "save-strip", "", false, "save original and registered images as GDAL GTiff") + procCmd.Flags().StringVarP(¶msXML, "params", "x", "", "params xml file path") } diff --git a/cmd/root.go b/cmd/root.go index d6fb05b..8ad5473 100644 --- a/cmd/root.go +++ b/cmd/root.go @@ -8,6 +8,7 @@ import ( var ( configFile string + paramsXML string ) var rootCmd = &cobra.Command{ diff --git a/go.mod b/go.mod index b8bc872..6f7e895 100644 --- a/go.mod +++ b/go.mod @@ -10,6 +10,7 @@ require ( require ( github.com/fsnotify/fsnotify v1.7.0 // indirect + github.com/gogo/protobuf v1.3.2 // indirect github.com/hashicorp/hcl v1.0.0 // indirect github.com/inconshreveable/mousetrap v1.1.0 // indirect github.com/lestrrat-go/strftime v1.0.6 // indirect @@ -18,6 +19,8 @@ require ( github.com/mattn/go-isatty v0.0.17 // indirect github.com/mgutz/ansi v0.0.0-20200706080929-d51e80ef957d // indirect github.com/mitchellh/mapstructure v1.5.0 // indirect + github.com/paulmach/orb v0.11.1 // indirect + github.com/paulmach/protoscan v0.2.1 // indirect github.com/pelletier/go-toml/v2 v2.1.0 // indirect github.com/pkg/errors v0.9.1 // indirect github.com/sagikazarmark/locafero v0.4.0 // indirect @@ -28,10 +31,12 @@ require ( github.com/spf13/pflag v1.0.5 // indirect github.com/spf13/viper v1.18.2 // indirect github.com/subosito/gotenv v1.6.0 // indirect + github.com/twpayne/go-geom v1.5.4 // indirect + go.mongodb.org/mongo-driver v1.11.4 // indirect go.uber.org/atomic v1.9.0 // indirect go.uber.org/multierr v1.9.0 // indirect golang.org/x/crypto v0.23.0 // indirect - golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa // indirect + golang.org/x/exp v0.0.0-20240222234643-814bf88cf225 // indirect golang.org/x/term v0.20.0 // indirect golang.org/x/text v0.15.0 // indirect gopkg.in/ini.v1 v1.67.0 // indirect diff --git a/go.sum b/go.sum index 4d4efec..ffa353d 100644 --- a/go.sum +++ b/go.sum @@ -887,6 +887,7 @@ github.com/go-pdf/fpdf v0.5.0/go.mod h1:HzcnA+A23uwogo0tp9yU+l3V+KXhiESpt1PMayhO github.com/go-pdf/fpdf v0.6.0/go.mod h1:HzcnA+A23uwogo0tp9yU+l3V+KXhiESpt1PMayhOh5M= github.com/goccy/go-json v0.9.11/go.mod h1:6MelG93GURQebXPDq3khkgXZkazVtN9CRI+MGFi0w8I= github.com/godbus/dbus/v5 v5.0.4/go.mod h1:xhWf0FNVPg57R7Z0UbKHbJfkEywrmjJnf7w5xrFpKfA= +github.com/gogo/protobuf v1.3.2 h1:Ov1cvc58UF3b5XjBnZv7+opcTcQFZebYjWzi34vdm4Q= github.com/gogo/protobuf v1.3.2/go.mod h1:P1XiOD3dCwIKUDQYPy72D8LYyHL2YPYrpS2s69NZV8Q= github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0/go.mod h1:E/TSTwGwJL78qG/PmXZO1EjYhfJinVAhrmmHX6Z8B9k= github.com/golang/glog v0.0.0-20160126235308-23def4e6c14b/go.mod h1:SBH7ygxi8pfUlaOkMMuAQtPIUF8ecWP5IEl/CR7VP2Q= @@ -925,6 +926,8 @@ github.com/golang/protobuf v1.5.1/go.mod h1:DopwsBzvsk0Fs44TXzsVbJyPhcCPeIwnvohx github.com/golang/protobuf v1.5.2/go.mod h1:XVQd3VNwM+JqD3oG2Ue2ip4fOMUkwXdXDdiuN0vRsmY= github.com/golang/protobuf v1.5.3 h1:KhyjKVUg7Usr/dYsdSqoFveMYd5ko72D+zANwlG1mmg= github.com/golang/protobuf v1.5.3/go.mod h1:XVQd3VNwM+JqD3oG2Ue2ip4fOMUkwXdXDdiuN0vRsmY= +github.com/golang/protobuf v1.5.4 h1:i7eJL8qZTpSEXOPTxNKhASYpMn+8e5Q6AdndVa1dWek= +github.com/golang/snappy v0.0.1/go.mod h1:/XxbfmMg8lxefKM7IXC3fBNl/7bRcc72aCRzEWrmP2Q= github.com/golang/snappy v0.0.3/go.mod h1:/XxbfmMg8lxefKM7IXC3fBNl/7bRcc72aCRzEWrmP2Q= github.com/golang/snappy v0.0.4/go.mod h1:/XxbfmMg8lxefKM7IXC3fBNl/7bRcc72aCRzEWrmP2Q= github.com/google/btree v0.0.0-20180813153112-4030bb1f1f0c/go.mod h1:lNA+9X1NB3Zf8V7Ke586lFgjr2dZNuvo3lPJSGZ5JPQ= @@ -1057,6 +1060,7 @@ github.com/kballard/go-shellquote v0.0.0-20180428030007-95032a82bc51/go.mod h1:C github.com/kisielk/errcheck v1.5.0/go.mod h1:pFxgyoBC7bSaBwPgfKdkLd5X25qrDl4LWUI2bnpBCr8= github.com/kisielk/gotool v1.0.0/go.mod h1:XhKaO+MFFWcvkIS/tQcRk01m1F5IRFswLeQ+oQHNcck= github.com/klauspost/asmfmt v1.3.2/go.mod h1:AG8TuvYojzulgDAMCnYn50l/5QV3Bs/tp6j0HLHbNSE= +github.com/klauspost/compress v1.13.6/go.mod h1:/3/Vjq9QcHkK5uEr5lBEmyoZ1iFhe47etQ6QUkpK6sk= github.com/klauspost/compress v1.15.9/go.mod h1:PhcZ0MbTNciWF3rruxRgKxI5NkcHHrHUDtV4Yw2GlzU= github.com/klauspost/cpuid/v2 v2.0.9/go.mod h1:FInQzS24/EEf25PyTYn52gqo7WaD8xa0213Md/qVLRg= github.com/kr/fs v0.1.0/go.mod h1:FFnZGqtBN9Gxj7eW1uZ42v5BccTP0vu6NEaFoC2HwRg= @@ -1106,7 +1110,12 @@ github.com/mitchellh/mapstructure v1.5.0/go.mod h1:bFUtVrKA4DC2yAKiSyO/QUcy7e+RR github.com/modern-go/concurrent v0.0.0-20180228061459-e0a39a4cb421/go.mod h1:6dJC0mAP4ikYIbvyc7fijjWJddQyLn8Ig3JB5CqoB9Q= github.com/modern-go/reflect2 v0.0.0-20180701023420-4b7aa43c6742/go.mod h1:bx2lNnkwVCuqBIxFjflWJWanXIb3RllmbCylyMrvgv0= github.com/modern-go/reflect2 v1.0.1/go.mod h1:bx2lNnkwVCuqBIxFjflWJWanXIb3RllmbCylyMrvgv0= +github.com/montanaflynn/stats v0.0.0-20171201202039-1bf9dbcd8cbe/go.mod h1:wL8QJuTMNUDYhXwkmfOly8iTdp5TEcJFWZD2D7SIkUc= github.com/pascaldekloe/goe v0.0.0-20180627143212-57f6aae5913c/go.mod h1:lzWF7FIEvWOWxwDKqyGYQf6ZUaNfKdP144TG7ZOy1lc= +github.com/paulmach/orb v0.11.1 h1:3koVegMC4X/WeiXYz9iswopaTwMem53NzTJuTF20JzU= +github.com/paulmach/orb v0.11.1/go.mod h1:5mULz1xQfs3bmQm63QEJA6lNGujuRafwA5S/EnuLaLU= +github.com/paulmach/protoscan v0.2.1 h1:rM0FpcTjUMvPUNk2BhPJrreDKetq43ChnL+x1sRg8O8= +github.com/paulmach/protoscan v0.2.1/go.mod h1:SpcSwydNLrxUGSDvXvO0P7g7AuhJ7lcKfDlhJCDw2gY= github.com/pelletier/go-toml v1.9.3 h1:zeC5b1GviRUyKYd6OJPvBU/mcVDVoL1OhT17FCt5dSQ= github.com/pelletier/go-toml v1.9.3/go.mod h1:u1nR/EPcESfeI/szUZKdtJ0xRNbUoANCkoOuaOx1Y+c= github.com/pelletier/go-toml/v2 v2.1.0 h1:FnwAJ4oYMvbT/34k9zzHuZNrhlz48GB3/s6at6/MHO4= @@ -1191,8 +1200,15 @@ github.com/stretchr/testify v1.8.4/go.mod h1:sz/lmYIOXD/1dqDmKjjqLyZ2RngseejIcXl github.com/subosito/gotenv v1.2.0/go.mod h1:N0PQaV/YGNqwC0u51sEeR/aUtSLEXKX9iv69rRypqCw= github.com/subosito/gotenv v1.6.0 h1:9NlTDc1FTs4qu0DDq7AEtTPNw6SVm7uBMsUCUjABIf8= github.com/subosito/gotenv v1.6.0/go.mod h1:Dk4QP5c2W3ibzajGcXpNraDfq2IrhjMIvMSWPKKo0FU= +github.com/tidwall/pretty v1.0.0/go.mod h1:XNkn88O1ChpSDQmQeStsy+sBenx6DDtFZJxhVysOjyk= +github.com/twpayne/go-geom v1.5.4 h1:b8fiZd0SsEmQEeUdz2atT6KggF1KHiaZIi3DGi5p+sI= +github.com/twpayne/go-geom v1.5.4/go.mod h1:Hw8RszQ2/d9Y/KfOm9CvUJo78BOoIA5g0e4P7JCVKvo= github.com/x-cray/logrus-prefixed-formatter v0.5.2 h1:00txxvfBM9muc0jiLIEAkAcIMJzfthRT6usrui8uGmg= github.com/x-cray/logrus-prefixed-formatter v0.5.2/go.mod h1:2duySbKsL6M18s5GU7VPsoEPHyzalCE06qoARUCeBBE= +github.com/xdg-go/pbkdf2 v1.0.0/go.mod h1:jrpuAogTd400dnrH08LKmI/xc1MbPOebTwRqcT5RDeI= +github.com/xdg-go/scram v1.1.1/go.mod h1:RaEWvsqvNKKvBPvcKeFjrG2cJqOkHTiyTpzz23ni57g= +github.com/xdg-go/stringprep v1.0.3/go.mod h1:W3f5j4i+9rC0kuIEJL0ky1VpHXQU3ocBgklLGvcBnW8= +github.com/youmark/pkcs8 v0.0.0-20181117223130-1be2e3e5546d/go.mod h1:rHwXgn7JulP+udvsHwJoVG1YGAP6VLg4y9I5dyZdqmA= github.com/yuin/goldmark v1.1.25/go.mod h1:3hX8gzYuyVAZsxl0MRgGTJEmQBFcNTphYh9decYSb74= github.com/yuin/goldmark v1.1.27/go.mod h1:3hX8gzYuyVAZsxl0MRgGTJEmQBFcNTphYh9decYSb74= github.com/yuin/goldmark v1.1.32/go.mod h1:3hX8gzYuyVAZsxl0MRgGTJEmQBFcNTphYh9decYSb74= @@ -1205,6 +1221,8 @@ github.com/zeebo/xxh3 v1.0.2/go.mod h1:5NWz9Sef7zIDm2JHfFlcQvNekmcEl9ekUZQQKCYaD go.etcd.io/etcd/api/v3 v3.5.0/go.mod h1:cbVKeC6lCfl7j/8jBhAK6aIYO9XOjdptoxU/nLQcPvs= go.etcd.io/etcd/client/pkg/v3 v3.5.0/go.mod h1:IJHfcCEKxYu1Os13ZdwCwIUTUVGYTSAM3YSwc9/Ac1g= go.etcd.io/etcd/client/v2 v2.305.0/go.mod h1:h9puh54ZTgAKtEbut2oe9P4L/oqKCVB6xsXlzd7alYQ= +go.mongodb.org/mongo-driver v1.11.4 h1:4ayjakA013OdpGyL2K3ZqylTac/rMjrJOMZ1EHizXas= +go.mongodb.org/mongo-driver v1.11.4/go.mod h1:PTSz5yu21bkT/wXpkS7WR5f0ddqw5quethTUn9WM+2g= go.opencensus.io v0.21.0/go.mod h1:mSImk1erAIZhrmZN+AvHh14ztQfjbGwt4TtuofqLduU= go.opencensus.io v0.22.0/go.mod h1:+kGneAE2xo2IficOXnaByMWTGM9T73dGwxeWcUqIpI8= go.opencensus.io v0.22.2/go.mod h1:yxeiOL68Rb0Xd1ddK5vPZ/oVn4vY4Ynel7k9FzqtOIw= @@ -1237,6 +1255,7 @@ golang.org/x/crypto v0.0.0-20210421170649-83a5a9bb288b/go.mod h1:T9bdIzuCu7OtxOm golang.org/x/crypto v0.0.0-20210921155107-089bfa567519/go.mod h1:GvvjBRRGRdwPK5ydBHafDWAxML/pGHZbMvKqRZ5+Abc= golang.org/x/crypto v0.0.0-20211108221036-ceb1ce70b4fa/go.mod h1:GvvjBRRGRdwPK5ydBHafDWAxML/pGHZbMvKqRZ5+Abc= golang.org/x/crypto v0.0.0-20220314234659-1baeb1ce4c0b/go.mod h1:IxCIyHEi3zRg3s0A5j5BB6A9Jmi73HwBIUl50j+osU4= +golang.org/x/crypto v0.0.0-20220622213112-05595931fe9d/go.mod h1:IxCIyHEi3zRg3s0A5j5BB6A9Jmi73HwBIUl50j+osU4= golang.org/x/crypto v0.1.0/go.mod h1:RecgLatLF4+eUMCP1PoPZQb+cVrJcOPbHkTkbkB9sbw= golang.org/x/crypto v0.7.0/go.mod h1:pYwdfH91IfpZVANVyUOhSIPZaFoJGxTFbZhFTx+dXZU= golang.org/x/crypto v0.9.0/go.mod h1:yrmDGqONDYtNj3tH8X9dzUun2m2lzPa9ngI6/RUPGR0= @@ -1263,6 +1282,8 @@ golang.org/x/exp v0.0.0-20200224162631-6cc2880d07d6/go.mod h1:3jZMyOhIsHpP37uCMk golang.org/x/exp v0.0.0-20220827204233-334a2380cb91/go.mod h1:cyybsKvd6eL0RnXn6p/Grxp8F5bW7iYuBgsNCOHpMYE= golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa h1:FRnLl4eNAQl8hwxVVC17teOw8kdjVDVAiFMtgUdTSRQ= golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa/go.mod h1:zk2irFbV9DP96SEBUUAy67IdHUaZuSnrz1n472HUCLE= +golang.org/x/exp v0.0.0-20240222234643-814bf88cf225 h1:LfspQV/FYTatPTr/3HzIcmiUFH7PGP+OQ6mgDYo3yuQ= +golang.org/x/exp v0.0.0-20240222234643-814bf88cf225/go.mod h1:CxmFvTBINI24O/j8iY7H1xHzx2i4OsyguNBmN/uPtqc= golang.org/x/image v0.0.0-20180708004352-c73c2afc3b81/go.mod h1:ux5Hcp/YLpHSI86hEcLt0YII63i6oz57MZXIpbrjZUs= golang.org/x/image v0.0.0-20190227222117-0694c2d4d067/go.mod h1:kZ7UVZpmo3dzQBMxlp+ypCbDeSB+sBbTgSJuh5dn5js= golang.org/x/image v0.0.0-20190802002840-cff245a6509b/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0= @@ -1910,6 +1931,7 @@ google.golang.org/genproto/googleapis/rpc v0.0.0-20230807174057-1744710a1577/go. google.golang.org/genproto/googleapis/rpc v0.0.0-20230822172742-b8732ec3820d h1:uvYuEyMHKNt+lT4K3bN6fGswmK8qSvcreM3BwjDh+y4= google.golang.org/genproto/googleapis/rpc v0.0.0-20230822172742-b8732ec3820d/go.mod h1:+Bk1OCOj40wS2hwAMA+aCW9ypzm63QTBBHp6lQ3p+9M= google.golang.org/genproto/googleapis/rpc v0.0.0-20231120223509-83a465c0220f h1:ultW7fxlIvee4HYrtnaRPon9HpEgFk5zYpmfMgtKB5I= +google.golang.org/genproto/googleapis/rpc v0.0.0-20240311173647-c811ad7063a7 h1:8EeVk1VKMD+GD/neyEHGmz7pFblqPjHoi+PGQIlLx2s= google.golang.org/grpc v1.19.0/go.mod h1:mqu4LbDTu4XGKhr4mRzUsmM4RtVoemTSY81AxZiDr8c= google.golang.org/grpc v1.20.1/go.mod h1:10oTOabMzJvdu6/UiuZezV6QK5dSlG84ov/aaiqXj38= google.golang.org/grpc v1.21.1/go.mod h1:oYelfM1adQP15Ek0mdvEgi9Df8B9CZIaU1084ijfRaM= @@ -1958,6 +1980,7 @@ google.golang.org/grpc v1.57.0/go.mod h1:Sd+9RMTACXwmub0zcNY2c4arhtrbBYD1AUHI/dt google.golang.org/grpc v1.58.0 h1:32JY8YpPMSR45K+c3o6b8VL73V+rR8k+DeMIr4vRH8o= google.golang.org/grpc v1.58.0/go.mod h1:tgX3ZQDlNJGU96V6yHh1T/JeoBQ2TXdr43YbYSsCJk0= google.golang.org/grpc v1.59.0 h1:Z5Iec2pjwb+LEOqzpB2MR12/eKFhDPhuqW91O+4bwUk= +google.golang.org/grpc v1.62.1 h1:B4n+nfKzOICUXMgyrNd19h/I9oH0L1pizfk1d4zSgTk= google.golang.org/grpc/cmd/protoc-gen-go-grpc v1.1.0/go.mod h1:6Kw0yEErY5E/yWrBtf03jp27GLLJujG4z/JK95pnjjw= google.golang.org/protobuf v0.0.0-20200109180630-ec00e32a8dfd/go.mod h1:DFci5gLYBciE7Vtevhsrf46CRTquxDuWsQurQQe4oz8= google.golang.org/protobuf v0.0.0-20200221191635-4d8936d0db64/go.mod h1:kwYJMbMJ01Woi6D6+Kah6886xMZcty6N08ah7+eCXa0= @@ -1978,6 +2001,7 @@ google.golang.org/protobuf v1.29.1/go.mod h1:HV8QOd/L58Z+nl8r43ehVNZIU/HEI6OcFqw google.golang.org/protobuf v1.30.0/go.mod h1:HV8QOd/L58Z+nl8r43ehVNZIU/HEI6OcFqwMG9pJV4I= google.golang.org/protobuf v1.31.0 h1:g0LDEJHgrBl9N9r17Ru3sqWhkIx2NB67okBHPwC7hs8= google.golang.org/protobuf v1.31.0/go.mod h1:HV8QOd/L58Z+nl8r43ehVNZIU/HEI6OcFqwMG9pJV4I= +google.golang.org/protobuf v1.33.0 h1:uNO2rsAINq/JlFpSdYEKIZ0uKD/R9cpdv0T+yoGwGmI= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= gopkg.in/check.v1 v1.0.0-20180628173108-788fd7840127/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= gopkg.in/check.v1 v1.0.0-20201130134442-10cb98267c6c/go.mod h1:JHkPIbrfpd72SG/EVd6muEfDQjcINNoR0C8j2r3qZ4Q= diff --git a/pkg/auxilary/aux.go b/pkg/auxilary/aux.go new file mode 100644 index 0000000..85b727b --- /dev/null +++ b/pkg/auxilary/aux.go @@ -0,0 +1,46 @@ +package auxilary + +import ( + "os" + + log "github.com/sirupsen/logrus" +) + +// 卫星时间起点 北京时间 2000-01-01 20:00:00 +var ( + ReferenceTime2000 = 946728000 +) + +func ExtractAux(auxfile string) ([]*AuxFrameHead, []*AuxFocalBox, []*AuxPlatform, error) { + data, err := os.ReadFile(auxfile) + if err != nil { + log.Println("read aux data from", auxfile, "error:", err.Error()) + return nil, nil, nil, err + } + + var afh []*AuxFrameHead + var afb []*AuxFocalBox + var aps []*AuxPlatform + + for i := 0; i < len(data); i += 24 + 128 + 512 { + if i+24+128+512 > len(data) { + break + } + var head AuxFrameHead + head.Decode(data[i : i+24]) + afh = append(afh, &head) + + var box AuxFocalBox + box.Decode(data[i+24 : i+24+128]) + afb = append(afb, &box) + + var plat AuxPlatform + plat.Decode(data[i+24+128 : i+24+128+512]) + aps = append(aps, &plat) + } + + log.Printf("extract %d aux frame heads, %d aux focal boxes, %d aux platforms from %s", + len(afh), len(afb), len(aps), auxfile) + + return afh, afb, aps, err +} diff --git a/pkg/auxilary/aux_ebox.go b/pkg/auxilary/aux_ebox.go new file mode 100644 index 0000000..4fdb4e5 --- /dev/null +++ b/pkg/auxilary/aux_ebox.go @@ -0,0 +1,134 @@ +package auxilary + +import ( + "fmt" + + "github.com/k0kubun/pp/v3" +) + +// 焦面电箱辅助数据 128 字节 (每 16 行原始图像数据为一组) +type AuxFocalBox struct { + TransferTime float32 // 执行转移时间 0.5us/bit + TrainingDone byte // 1B成功；0B失败 + WorkMode byte // 工作模式 0B正常；1B故障 + IntegralDirection byte // 积分方向 1B正向；0B反向 + PGAGain byte // PGA增益 + PIntegrationLevel uint8 // [5]P波积分级数 + B1IntegrationLevel uint8 // [6.(7~6)]B1积分级数 + B2IntegrationLevel uint8 // [6.(5~4)]B2积分级数 + B3IntegrationLevel uint8 // [6.(3~2)]B3积分级数 + B4IntegrationLevel uint8 // [6.(1~0)]B4积分级数 + SecPluseState byte // [7.(7)]秒脉冲状态 0x1 成功；0x0 失败 + DarkFieldBias uint16 // [7.(5~0)-8] 暗场偏置 1DN/bit，十六进制无符号整型数 + PWinAddr uint16 // [9-10]全色开窗地址 + B1WinAddr uint16 // [11-12]B1开窗地址 + B2WinAddr uint16 // [13-14]B2开窗地址 + B3WinAddr uint16 // [15-16]B3开窗地址 + B4WinAddr uint16 // [17-18]B4开窗地址 + CCDWorkMode byte // [38]工作模式 00H=TDI推扫;11H=框幅; + RawDiskAvailableCap uint16 // [39-40]原始磁盘可用存储容量 量纲：2MBytes/bit + ZipDiskAvailableCap uint16 // [41-42]压缩盘可用存储容量 + CommandCount uint8 // [53]指令计数器 + LastCommandCode uint8 // [54]最后指令代码 +} + +func (ab *AuxFocalBox) Decode(data []byte) error { + if len(data) < 128 { + return fmt.Errorf("eletric box aux data length error") + } + byte0 := data[0] & 0b00111111 + ab.TransferTime = float32(uint32(byte0)<<16|uint32(data[1])<<8|uint32(data[2])) * 0.5 + ab.TrainingDone = data[3] & 0b10000000 + ab.WorkMode = data[3] & 0b01000000 + ab.IntegralDirection = data[3] & 0b00100000 + ab.PGAGain = data[3] & 0b00011111 + ab.PIntegrationLevel = data[4] + ab.B1IntegrationLevel = data[5] >> 6 & 0x03 + ab.B2IntegrationLevel = data[5] >> 4 & 0x03 + ab.B3IntegrationLevel = data[5] >> 2 & 0x03 + ab.B4IntegrationLevel = data[5] & 0x03 + ab.SecPluseState = data[6] >> 7 & 0x01 + ab.DarkFieldBias = uint16(data[6]&0x3f)<<8 | uint16(data[7]) + ab.PWinAddr = uint16(data[8])<<8 | uint16(data[9]) + ab.B1WinAddr = uint16(data[10])<<8 | uint16(data[11]) + ab.B2WinAddr = uint16(data[12])<<8 | uint16(data[13]) + ab.B3WinAddr = uint16(data[14])<<8 | uint16(data[15]) + ab.B4WinAddr = uint16(data[16])<<8 | uint16(data[17]) + ab.CCDWorkMode = data[37] + ab.RawDiskAvailableCap = (uint16(data[38])<<8 | uint16(data[39])) * 2 + ab.ZipDiskAvailableCap = (uint16(data[40])<<8 | uint16(data[41])) * 2 + ab.CommandCount = data[52] + ab.LastCommandCode = data[53] + return nil +} + +func (ab AuxFocalBox) Print() { + pp.Println(ab) +} + +func (ab AuxFocalBox) String() string { + return fmt.Sprintf( + `执行转移时间：%f us\n + TrainingDone：%x\n + 工作模式：%x\n + 积分方向：%x\n + PGA增益：%s + `, + ab.TransferTime, + ab.TrainingDone, + ab.WorkMode, + ab.IntegralDirection, + ab.PGAGainValue(), + ) +} + +func (ab AuxFocalBox) PGAGainValue() string { + switch ab.PGAGain { + case 0b0000: + return "0.75x" + case 0b0001: + return "1.25x" + case 0b0010: + return "1.75x" + case 0b0011: + return "2.25x" + case 0b0100: + return "2.75x" + case 0b0101: + return "3.25x" + case 0b0110: + return "3.75x" + case 0b0111: + return "4.25x" + case 0b1000: + return "4.75x" + case 0b1001: + return "5.25x" + case 0b1010: + return "5.75x" + case 0b10000: + return "1x" + case 0b10001: + return "1.5x" + case 0b10010: + return "2x" + case 0b10011: + return "2.5x" + case 0b10100: + return "3x" + case 0b10101: + return "3.5x" + case 0b10110: + return "4x" + case 0b10111: + return "4.5x" + case 0b11000: + return "5x" + case 0b11001: + return "5.5x" + case 0b11010: + return "6x" + default: + return "unknown" + } +} diff --git a/pkg/auxilary/aux_elements.go b/pkg/auxilary/aux_elements.go new file mode 100644 index 0000000..d9dc07a --- /dev/null +++ b/pkg/auxilary/aux_elements.go @@ -0,0 +1,420 @@ +package auxilary + +var AuxHeader = []interface{}{ + "包头信息", + "填充位", + "帧头信息", + "当前图像模式", + "B4谱状态", + "B3谱状态", + "B2谱状态", + "B1谱状态", + "全色状态", + "当前行流水号", + "时间秒", + "秒小数", + "文件号", + "执行行转移时间", + "Trainingdone", + "工作模式", + "积分方向", + "PGA增益", + "P谱段积分级数", + "B1积分级数", + "B2积分级数", + "B3积分级数", + "B4积分级数", + "秒脉冲状态", + "保留", + "暗场偏置", + "全色开窗地址", + "多光谱1开窗地址", + "多光谱2开窗地址", + "多光谱3开窗地址", + "多光谱4开窗地址", + "面阵模式Linetime时钟周期数", + "面阵模式开窗地址", + "面阵模式开窗大小", + "面阵曝光时间粗调EXP_C", + "面阵曝光时间精调EXP_F", + "面阵模式最小读出行", + "硬盘1温度", + "硬盘2温度", + "保留", + "传感器温度", + "FPGA逻辑版本号", + "工作模式", + "原始盘可用存储容量", + "压缩盘可用存储容量", + "原始盘状态", + "原始盘Host初始化状态", + "原始盘SATA控制器状态", + "原始盘SATA错误计数", + "压缩盘状态", + "压缩盘Host初始化状态", + "压缩盘SATA控制器状态", + "压缩盘SATA错误计数", + "保留", + "DDR初始化状态", + "原始图像硬盘状态", + "压缩数据硬盘状态", + "硬盘1读写状态", + "硬盘2读写状态", + "硬盘1初始化状态", + "硬盘2初始化状态", + "保留", + "保留", + "硬盘1禁用标志", + "硬盘2禁用标志", + "保留", + "B2数据移位", + "B1数据移位", + "B4数据移位", + "B3数据移位", + "指令计数", + "最后一条指令编码", + "指令接收状态", + "错误指令计数", + "错误指令帧编号", + "保留", + "传感器数字电路温度", + "卫星UTC时间秒", + "卫星复波道标志位", + "卫星秒小数", + "定姿四元数(J2000)q0", + "定姿四元数(J2000)q1", + "定姿四元数(J2000)q2", + "定姿四元数(J2000)q3", + "本体相对轨道四元数矢部q1", + "本体相对轨道四元数矢部q2", + "本体相对轨道四元数矢部q3", + "轨道相对惯性系四元数矢部q1", + "轨道相对惯性系四元数矢部q2", + "轨道相对惯性系四元数矢部q3", + "本体相对轨道姿态角1", + "本体相对轨道姿态角2", + "本体相对轨道姿态角3", + "本体相对轨道角速度1", + "本体相对轨道角速度2", + "本体相对轨道角速度3", + "模式运行时间", + "姿控调用周期", + "姿轨控部件使用标志", + "姿轨控算法执行标记", + "偏差四元数q1", + "偏差四元数q2", + "偏差四元数q3", + "偏差角速度1", + "偏差角速度2", + "偏差角速度3", + "X飞轮估计摩擦力矩", + "Y飞轮估计摩擦力矩", + "Z飞轮估计摩擦力矩", + "三轴陀螺X轴角速度漂移估计", + "三轴陀螺Y轴角速度漂移估计", + "三轴陀螺Z轴角速度漂移估计", + "三轴角速度卡尔曼漂移X", + "三轴角速度卡尔曼漂移Y", + "三轴角速度卡尔曼漂移Z", + "X轴估计环境干扰力矩", + "Y轴估计环境干扰力矩", + "Z轴估计环境干扰力矩", + "X轴计算飞轮控制力矩", + "Y轴计算飞轮控制力矩", + "Z轴计算飞轮控制力矩", + "期望四元数矢部1", + "期望四元数矢部2", + "期望四元数矢部3", + "期望角速度1", + "期望角速度2", + "期望角速度3", + "计算当前J2000位置X", + "计算当前J2000位置Y", + "计算当前J2000位置Z", + "计算当前J2000速度X", + "计算当前J2000速度Y", + "计算当前J2000速度Z", + "计算当前84位置X", + "计算当前84位置Y", + "计算当前84位置Z", + "计算当前84速度X", + "计算当前84速度Y", + "计算当前84速度Z", + "偏流角", + "数传点经度", + "数传点纬度", + "数传点地程高", + "定姿方式标志1", + "定姿方式标志2", + "定姿方式标志3", + "X飞轮期望转速", + "Y飞轮期望转速", + "Z飞轮期望转速", + "保留", + "保留", + "保留", + "保留", + "锂电池温度", + "相机主镜温度", + "相机次镜温度", + "推进模块温度", + "负X侧相机桁架杆温度", + "正X正Y侧相机桁架杆温度", + "正X负Y侧相机桁架杆温度", + "正X侧相机支撑背板温度", + "负X侧相机支撑背板温度", + "星敏支架温度", + "成像电箱温度", + "正Y帆板温度", + "电源下位机温度", + "配电热控驱动温度", + "电源控制器温度", + "数字太阳敏矢量数据有效位", + "数字太阳敏矢量数据X", + "数字太阳敏矢量数据Y", + "数字太阳敏位置X1", + "数字太阳敏位置X2", + "数字太阳敏位置Y1", + "数字太阳敏位置Y2", + "太阳敏温度", + "数字太阳敏当前正在使用的阈值（源码）", + "数字太阳敏当前正在使用的增益", + "保留", + "保留", + "保留", + "错误码计数", + "单粒子错误计数", + "配电错误码1", + "配电错误码2", + "配电错误码3", + "配电错误码4", + "配电错误码5", + "GPS天内秒", + "GPSUTC时间累计秒", + "太阳阵电流", + "母线电流", + "负载电流", + "蓄电池电压", + "电源母线电压", + "CPU5.2V电压遥测值", + "5.2V配电电压遥测值", + "保留", + "蓄电池组当前电量", + "模式运行时间秒", + "卫星现运行模式", + "组合业务标识", + "当前业务状态", + "中心机指令接收总计数", + "中心机错误指令计数", + "执行指令所在分系统", + "执行指令的指令代码", + "执行延时指令总计数", + "当前延时指令计数", + "执行延时指令所在分系统", + "执行延时指令的指令代码", + "当前延时业务计数", + "成功执行业务计数", + "异常中止业务计数", + "指令执行状态", + "业务异常中止状态", + "测控数传一体机通信状态", + "保留", + "北斗短报文通信状态", + "GPS接收机通信状态", + "数字太阳敏通信状态", + "星敏1通信状态", + "星敏2通信状态", + "光纤陀螺通信状态", + "MEMS陀螺通信状态", + "飞轮1通信状态", + "飞轮2通信状态", + "飞轮3通信状态", + "飞轮4通信状态", + "智能载荷通信状态", + "保留", + "保留", + "电磁阀开关状态", + "业务保存状态", + "卫星类型标识", + "卫星序号标识", + "保留", + "保留", + "保留", + "保留", + "锂电池加热器通断状态", + "相机主镜加热器通断状态", + "相机次镜加热器通断状态", + "推进模块加热器通断状态", + "相机负X侧桁架杆加热器通断状态", + "成像电箱加热器通断状态", + "相机正X正Y侧桁架杆加热器通断状态", + "相机正X负Y侧桁架杆加热器通断状态", + "相机正X侧支撑背板加热器通断状态", + "相机负X侧支撑背板加热器通断状态", + "星敏支架加热器通断状态", + "保留", + "温度修正系数校验状态", + "电源下位机广播帧监视功能", + "当前控温码表", + "默认控温码表", + "飞轮1电源供电状态", + "飞轮2电源供电状态", + "飞轮3电源供电状态", + "SADA1电源供电状态", + "SADA2电源供电状态", + "测控数传电源供电状态", + "保留", + "保留", + "北斗短报文供电状态", + "推进电源供电状态", + "焦面电源供电状态", + "飞轮4电源供电状态", + "星敏1电源供电状态", + "星敏2电源供电状态", + "数字太阳敏电源供电状态", + "导航电源供电状态", + "三轴光纤陀螺电源供电状态", + "MEMS陀螺电源供电状态", + "热控正线状态", + "热控1状态", + "热控2状态", + "保留", + "保留", + "保留", + "锂电池加热器控温模式", + "相机主镜加热器控温模式", + "相机次镜加热器控温模式", + "推进模块加热器控温模式", + "相机负X侧桁架杆加热器控温模式", + "相机正X正Y侧桁架杆加热器控温模式", + "相机正X负Y侧桁架杆加热器控温模式", + "相机正X侧支撑背板加热器控温模式", + "相机负X侧支撑背板加热器控温模式", + "星敏支架加热器控温模式", + "成像电箱加热器控温模式", + "保留", + "保留", + "保留", + "保留", + "保留", + "保留", + "保留", + "接收机时间来源", + "定位模式", + "轨道数据有效标示", + "主备机标志", + "PPS状态", + "GPS最高信噪比", + "BD最高信噪比", + "参与定位的GPS导航星数", + "参与定位的BD导航星数", + "GPS几何精度因子", + "GPS连续工作时间", + "保留", + "保留", + "WGS-84系X位置", + "WGS-84系Y位置", + "WGS-84系Z位置", + "WGS-84系X速度", + "WGS-84系Y速度", + "WGS-84系Z速度", + "J2000系X位置", + "J2000系Y位置", + "J2000系Z位置", + "J2000系X速度", + "J2000系Y速度", + "J2000系Z速度", + "三轴光纤陀螺X轴角速度", + "三轴光纤陀螺Y轴角速度", + "三轴光纤陀螺Z轴角速度", + "三轴光纤陀螺X轴角速度（1）", + "三轴光纤陀螺Y轴角速度（1）", + "三轴光纤陀螺Z轴角速度（1）", + "星敏AUTC时间", + "星敏AUTC秒小数", + "星敏A四元数q1", + "星敏A四元数q2", + "星敏A四元数q3", + "星敏A四元数q4", + "星敏A曝光时间", + "星敏A阈值", + "星敏A背景值", + "星敏A上电进入boot标志", + "星敏AEDAC打开标志", + "星敏A程序版本", + "星敏A四元数滤波标志", + "星敏A系统内部工作进程代号", + "星敏A系统工作模式", + "星敏A提取星数", + "星敏A四元数有效时导航星数", + "星敏A图像增益", + "星敏A识别星数", + "星敏A打开或者关断外部图像", + "星敏A姿态数据有效标志", + "星敏A内部软件版本号低3位", + "星敏A产品设备编号低5位", + "星敏A成像传感器温度", + "星敏A在轨EDAC错误计数", + "星敏A图像帧号", + "星敏A四星寻找阈值", + "星敏A跟踪阈值", + "星敏ASAA阈值", + "星敏ASAA工作模式", + "星敏A动态模式标志位", + "星敏AX方向角速度", + "星敏AY方向角速度", + "星敏AZ方向角速度", + "星敏A星点阈值自适应功能", + "保留", + "星敏A速率质量", + "星敏BUTC时间", + "星敏BUTC秒小数", + "星敏B四元数q1", + "星敏B四元数q2", + "星敏B四元数q3", + "星敏B四元数q4", + "星敏B曝光时间", + "星敏B阈值", + "星敏B背景值", + "星敏B上电进入boot标志", + "星敏BEDAC打开标志", + "星敏B程序版本", + "星敏B四元数滤波标志", + "星敏B系统内部工作进程代号", + "星敏B系统工作模式", + "星敏B提取星数", + "星敏B四元数有效时导航星数", + "星敏B图像增益", + "星敏B识别星数", + "星敏B打开或者关断外部图像", + "星敏B姿态数据有效标志", + "星敏B内部软件版本号低3位", + "星敏B产品设备编号低5位", + "星敏B成像传感器温度", + "星敏B在轨EDAC错误计数", + "星敏B图像帧号", + "星敏B四星寻找阈值", + "星敏B跟踪阈值", + "星敏BSAA阈值", + "星敏BSAA工作模式", + "星敏B动态模式标志位", + "保留", + "飞轮1转速", + "飞轮1当前电流", + "飞轮2转速", + "飞轮2当前电流", + "飞轮3转速", + "飞轮3当前电流", + "飞轮4转速", + "飞轮4当前电流", + "MEMS陀螺X方向角速度", + "MEMS陀螺Y方向角速度", + "MEMS陀螺Z方向角速度", + "X方向磁场强度", + "Y方向磁场强度", + "Z方向磁场强度", + "输入姿轨控数据UTC时间秒", + "输入姿轨控数据时间秒小数", + "保留", + "校验和", +} diff --git a/pkg/auxilary/aux_head.go b/pkg/auxilary/aux_head.go new file mode 100644 index 0000000..fa800ef --- /dev/null +++ b/pkg/auxilary/aux_head.go @@ -0,0 +1,111 @@ +package auxilary + +import ( + "encoding/binary" + "errors" +) + +// 每行传感器数据帧头信息长度为24字节 +type AuxFrameHead struct { + PkgHead [4]byte // 包头头标识符，固定为0xD15BD15B + FillByte0 byte // 填充字节 0x00 + FrmHead [6]byte // 帧头标识符，固定为0xEB90EB90EB90 + // 图像模式 + // B7 模式标识 000：为线阵模式 + // B6 + // B5 + // B4 B4谱状态 0：不输出1：输出 + // B3 B3谱状态 0：不输出1：输出 + // B2 B2谱状态 0：不输出1：输出 + // B1 B1谱状态 0：不输出1：输出 + // B0 全色状态 0：不输出1：输出 + ImgMode byte + SerialNo uint32 // 当前流水号 + TimeSec uint32 // 时间（秒） + TimeSecFrac uint32 // 秒小数 量纲：1us/bit，十六进制无符号整型数 + FileNo uint8 // 文件号 + + IsValidFrmHead bool + B0 bool + B1 bool + B2 bool + B3 bool + B4 bool + IsLinerMatrix bool + RowLength int +} + +func (afh *AuxFrameHead) Decode(data []byte) error { + if len(data) < 24 { + return errors.New("length of AuxFrameHead is not 24") + } + + afh.PkgHead = [4]byte{data[0], data[1], data[2], data[3]} + afh.FillByte0 = data[4] + afh.FrmHead = [6]byte{data[5], data[6], data[7], data[8], data[9], data[10]} + afh.ImgMode = data[11] + afh.B0 = (afh.ImgMode&(1<<0) != 0x0) + afh.B1 = (afh.ImgMode&(1<<1) != 0x0) + afh.B2 = (afh.ImgMode&(1<<2) != 0x0) + afh.B3 = (afh.ImgMode&(1<<3) != 0x0) + afh.B4 = (afh.ImgMode&(1<<4) != 0x0) + afh.IsLinerMatrix = (afh.ImgMode&(1<<5) == 0x0 && afh.ImgMode&(1<<6) == 0x0 && afh.ImgMode&(1<<7) == 0x0) + afh.SerialNo = binary.BigEndian.Uint32(data[12:16]) + afh.TimeSec = binary.BigEndian.Uint32(data[16:20]) + afh.TimeSecFrac = uint32(uint32(data[20])<<16 | uint32(data[21])<<8 | uint32(data[22])) + afh.FileNo = data[23] + + afh.RowLength = afh.LengthOfRow() + afh.IsValidFrmHead = afh.CheckFrmHead() + + return nil +} + +// 计算一行数据长度：帧头+辅助数据+图像数据 +func (afh AuxFrameHead) LengthOfRow() int { + if !afh.IsLinerMatrix { + return 14344 + } + + length := 64 // 帧头+辅助数据长度 + if afh.B0 { + length += 19040 + } + if afh.B1 { + length += 1192 + } + if afh.B2 { + length += 1192 + } + if afh.B3 { + length += 1192 + } + if afh.B4 { + length += 1192 + } + + return length +} + +func (afh AuxFrameHead) CheckFrmHead() bool { + if afh.FrmHead[0] == 0xEB && afh.FrmHead[1] == 0x90 && + afh.FrmHead[2] == 0xEB && afh.FrmHead[3] == 0x90 && + afh.FrmHead[4] == 0xEB && afh.FrmHead[5] == 0x90 { + return true + } + return false +} + +func (afh AuxFrameHead) ImageMode() string { + if afh.IsLinerMatrix { + return "线阵模式" + } + return "面阵模式" +} + +func (afh AuxFrameHead) BandStatus(b bool) string { + if b { + return "输出" + } + return "不输出" +} diff --git a/pkg/auxilary/aux_platform.go b/pkg/auxilary/aux_platform.go new file mode 100644 index 0000000..0f720e6 --- /dev/null +++ b/pkg/auxilary/aux_platform.go @@ -0,0 +1,390 @@ +package auxilary + +import ( + "encoding/binary" + "fmt" + "os" + "time" + + "github.com/k0kubun/pp/v3" + "starwiz.cn/sjy01/image-proc/pkg/calculator" + + log "github.com/sirupsen/logrus" +) + +const AuxPlatformFrmSize = 512 // 512 bytes + +// 中心机辅助数据 512 字节 (每 16 行原始图像数据为一组) +type AuxPlatform struct { + // 结构中标注的字节序号从 1 开始计数 + UTCTimeSec uint32 // [1-4] 卫星 UTC 时间戳（秒） + Waveway uint8 // [5] 波道 0x00：波道1（241.0~311.0）0x01：波道2（241.1~311.1） + Microsecond uint32 // [6-8]卫星秒小数（微秒） + QuatAttstarQ0 float64 // [9-12]定姿四元数(J2000) 的 Q0 值，量纲 1/0x40000000 + QuatAttstarQ1 float64 // [13-16]Q1 值 + QuatAttstarQ2 float64 // [17-20]Q2 值 + QuatAttstarQ3 float64 // [21-24]Q3 值 + QuatOrbitQ1 float64 // [25-28]本体相对轨道四元数矢部 的 Q1 值，量纲 1/0x40000000 + QuatOrbitQ2 float64 // [29-32]Q2 值 + QuatOrbitQ3 float64 // [33-36]Q3 值Q3 值 + QuatOrbJQ1 float64 // [37-40]轨道相对惯性系四元数矢部 的 Q1 值，量纲 1/0x40000000 + QuatOrbJQ2 float64 // [41-44]Q2 值 + QuatOrbJQ3 float64 // [45-48]Q3 值 + Eular1 float64 // [49-52]本体相对轨道姿态角，[Z]Yaw 测摆角，量纲：1/1000000 单位：rad + Eular2 float64 // [53-56] Pitch [Y] + Eular3 float64 // [57-60] Roll [X] + DotEular1 float64 // [61-62]本体相对轨道角速度，量纲：1/100000 单位：rad + DotEular2 float64 // [63-64] + DotEular3 float64 // [65-66] + ModTime uint32 // [67-70]模式运行时间 秒 + DTime float64 // [71-72]姿控调用周期，量纲：1/1000 单位：秒 + AutoState [3]uint8 // [73-75]姿轨控部件使用标志 + ProTrack [16]uint8 // [76-91]姿轨控算法执行标记 + QeQ1 float64 // [92-95]偏差四元数 double 量纲：1/0x40000000 /bit + QeQ2 float64 // [96-99] + QeQ3 float64 // [100-103] + We1 float64 // [104-105]偏差角速度 量纲：1/0x100000 单位：rad + We2 float64 // [106-107] + We3 float64 // [108-109] + WTFX float64 // [110-111]X飞轮估计摩擦力矩 量纲：1/0x400000 + WTFY float64 // [112-113] + WTFZ float64 // [114-115] + FbdriftX float64 // [116-117]三轴陀螺X轴角速度漂移估计 1/0x4000000 rad/s + FbdriftY float64 // [118-119] + FbdriftZ float64 // [120-121] + KalbX float64 // [122-123]三轴角速度卡尔曼漂移X 量纲：1/0x4000000 单位：rad/s + KalbY float64 // [124-125] + KalbZ float64 // [126-127] + HTDX float64 // [128-129]X轴估计环境干扰力矩 量纲：1/0x200000 + HTDY float64 // [130-131] + HTDZ float64 // [132-133] + CommandWheelX float64 // [134-135]X轴计算飞轮控制力矩 量纲：1/0x40000 + CommandWheelY float64 // [136-137] + CommandWheelZ float64 // [138-139] + QuatG1 float64 // [140-143]期望四元数矢部1 量纲：1/0x40000000 /bit + QuatG2 float64 // [144-147] + QuatG3 float64 // [148-151] + WG1 float64 // [152-153]期望角速度1 量纲：1/0x40000 + WG2 float64 // [154-155] + WG3 float64 // [156-157] + J2000PosX float64 // [158-161]计算当前J2000位置X 1/0x100 单位：m + J2000PosY float64 // [162-165] + J2000PosZ float64 // [166-169] + J2000VelX float64 // [170-173]计算当前J2000速度X 量纲：1/0x100 单位：m/s + J2000VelY float64 // [174-177] + J2000VelZ float64 // [178-181] + W84PosX float64 // [182-185]计算当前WGS 84位置X 1/0x100 单位：m + W84PosY float64 // [186-189] + W84PosZ float64 // [190-193] + W84VelX float64 // [194-197]计算当前WGS 84速度X 量纲：1/0x100 单位：m/s + W84VelY float64 // [198-201] + W84VelZ float64 // [202-205] + AngleDraft float64 // [206-209]偏流角 量纲：1/10000000 单位：rad + DataTransLong float64 // [210-211]数传点经度 量纲：1/1000 + DataTransLat float64 // [212-213]数传点纬度 量纲：1/1000 + DataTransH float64 // [214-215]数传点地程高 量纲：1 + StarPrio [3]byte // [216-218] 星敏1定姿方式标志1 0x00：星敏1定姿；0xA0:星敏1滤波;0x01星敏2定姿；0xA1:星敏2滤波;0xff:未使用星敏定姿 + WheelspeedcalX float64 // [219-222]X飞轮期望转速 量纲：1/10000 单位：rpm + WheelspeedcalY float64 // [223-226] + WheelspeedcalZ float64 // [227-230] + Reserved0 [2]byte // 预留字节 + Reserved1 [4]byte // 预留字节 + Reserved2 [2]byte // 预留字节 + Reserved3 [2]byte // 预留字节 + // 241.0 - 311.0 波道1参数 + // 241.1 - 311.1 波道2参数 + WGS84PosX float64 // [312-315]计算当前WGS 84位置X 单位：0.01米 + WGS84PosY float64 // [316-319] + WGS84PosZ float64 // [320-323] + WGS84VelX float64 // [324-327]计算当前WGS 84速度X 单位：0.01米/秒 + WGS84VelY float64 // [328-331] + WGS84VelZ float64 // [332-335] + J2000Pos_X float64 // [336-339]计算当前J2000位置X 单位：0.01米 + J2000Pos_Y float64 // [340-343] + J2000Pos_Z float64 // [344-347] + J2000Vel_X float64 // [348-351]计算当前J2000速度X 单位：0.01米/秒 + J2000Vel_Y float64 // [352-355] + J2000Vel_Z float64 // [356-359] + FOGyroXAV float64 // [360-362]三轴光纤陀螺X轴角速度 陀螺速度共3字节单位°/s量纲：1/72359 + FOGyroYAV float64 // [363-365] + FOGyroZAV float64 // [366-368] + FOGyroXAV1 float64 // [369-371] + FOGyroYAV1 float64 // [372-374] + FOGyroZAV1 float64 // [375-377] + SS1_UTCTime uint32 // [378-381] 星敏1 UTC时间 + SS1_UTCTimeFrac float32 // [382-384] 星敏1 UTC秒小数 单位为40.96 us，高字节在前 + SS1_Q1 float64 // [385-388] 星敏1四元数q1 当量：1/2147483647 （星敏坐标系相对于J2000惯性坐标系） + SS1_Q2 float64 // [389-392] 星敏1四元数q2 + SS1_Q3 float64 // [393-396] 星敏1四元数q3 + SS1_Q4 float64 // [397-400] 星敏1四元数q4 + SS1_ExposureTime uint8 // [401] 星敏1曝光时间 无符号数，单位ms + SS1_Threshold uint8 // [402] 星敏1阈值 无符号数 + SS1_BackgroudV uint8 // [403] 星敏1背景值 无符号数 + SS1_Flags byte // [404] 星敏1标志位 0x01：曝光时间有效；0x02：阈值有效；0x04：背景值有效 + SS1_WorkMode uint8 // [405.(7~6)] 星敏1工作模式: 0正常工作模式；1保留；2固定阈值模式；3测试模式（保留） + SS1_ExtractStars uint8 // [405.(5~0)] + SS1_NavStars uint8 // [406.(7~5)] 星敏1导航星数 + SS1_Gain uint8 // [406.(4~0)] 星敏1增益 + SS1_RegonizedStars uint8 // [407.(7~2)] 星敏1识别星数 + SS1_ExtenalImage bool // [407.(1)] 星敏1外部图像标志位 0x01：外部图像；0x00：关闭 + SS1_AttitudeActive bool // [407.(0)] 星敏1姿态有效标志位 0x01：激活；0x00：关闭 + SS1_ImgFrmNo uint32 // [411-413] 星敏1图像帧号 + SS1_XAV float64 // [417-418] 星敏1X方向角速度 单位：2e-11 °/s + SS1_YAV float64 // [419-420] + SS1_ZAV float64 // [421-422] + SS2_UTCTime uint32 // [424-427] 星敏2 UTC时间 + SS2_UTCTimeFrac float32 // [428-430] 星敏2 UTC秒小数 单位为40.96 us，高字节在前 + SS2_Q1 float64 // [431-434] 星敏2四元数q1 当量：1/2147483647 （星敏坐标系相对于J2000惯性坐标系） + SS2_Q2 float64 // [435-438] 星敏2四元数q2 + SS2_Q3 float64 // [439-442] 星敏2四元数q3 + SS2_Q4 float64 // [443-446] 星敏2四元数q4 + SS2_ExposureTime uint8 // [447] 星敏2曝光时间 无符号数，单位ms + SS2_Threshold uint8 // [448] 星敏2阈值 无符号数 + SS2_BackgroudV uint8 // [449] 星敏2背景值 无符号数 + SS2_Flags byte // [450] 星敏2标志位 0x01：曝光时间有效；0x02：阈值有效；0x04：背景值有效 + SS2_WorkMode uint8 // [451.(7~6)] 星敏2工作模式: 0正常工作模式；1保留；2固定阈值模式；3测试模式（保留） + SS2_ExtractStars uint8 // [451.(5~0)] + SS2_NavStars uint8 // [452.(7~5)] 星敏2导航星数 + SS2_Gain uint8 // [452.(4~0)] 星敏2增益 + SS2_RegonizedStars uint8 // [453.(7~2)] 星敏2识别星数 + SS2_ExtenalImage bool // [453.(1)] 星敏2外部图像标志位 0x01：外部图像；0x00：关闭 + SS2_AttitudeActive bool // [453.(0)] 星敏2姿态有效标志位 0x01：激活；0x00：关闭 + SS2_ImgFrmNo uint32 // [457-459] 星敏2图像帧号 + FlyWheel1_Vel float64 // [464-466] 飞轮1速度 单位：0.1rpm + FlyWheel1_Amps float64 // [467] 飞轮1电流 单位：0.009A + FlyWheel2_Vel float64 // [468-470] 飞轮2速度 单位：0.1rpm + FlyWheel2_Amps float64 // [471] 飞轮2电流 单位：0.009A + FlyWheel3_Vel float64 // [472-474] 飞轮3速度 单位：0.1rpm + FlyWheel3_Amps float64 // [475] 飞轮3电流 单位：0.009A + FlyWheel4_Vel float64 // [476-478] 飞轮4速度 单位：0.1rpm + FlyWheel4_Amps float64 // [479] 飞轮4电流 单位：0.009A + MemsGyroXAV float64 // [480-483] MEMS陀螺X轴角速度 0.000001°/s + MemsGyroYAV float64 // [484-487] MEMS陀螺Y轴角速度 0.000001°/s + MemsGyroZAV float64 // [488-491] MEMS陀螺Z轴角速度 0.000001°/s + MagnetictrengthX float64 // [492-495] 磁场强度X 量纲19200,单位uT + MagnetictrengthY float64 // [496-499] 磁场强度Y 量纲19200,单位uT + MagnetictrengthZ float64 // [500-503] 磁场强度Z 量纲19200,单位uT + ASSTimeInt uint32 // [504-507] 输入姿轨控数据UTC时间秒 + ASSTimeDec uint32 // [508-510] 输入姿轨控数据UTC时间秒小数 单位为μs + Reserved4 byte // [511] 保留字节 + CheckSum byte // [512] 校验和 +} + +func (ap *AuxPlatform) Decode(data []byte) error { + if len(data) < 512 { + return fmt.Errorf("data length less than 512") + } + + // 按需解析数据 + ap.UTCTimeSec = binary.BigEndian.Uint32(data[0:4]) + ap.Waveway = data[4] + ap.Microsecond = uint32(data[5])<<16 | uint32(data[6])<<8 | uint32(data[7]) + ap.QuatAttstarQ0 = float64(bigEndianToInt32(data[8:12])) / 0x40000000 + ap.QuatAttstarQ1 = float64(bigEndianToInt32(data[12:16])) / 0x40000000 + ap.QuatAttstarQ2 = float64(bigEndianToInt32(data[16:20])) / 0x40000000 + ap.QuatAttstarQ3 = float64(bigEndianToInt32(data[20:24])) / 0x40000000 + ap.QuatOrbitQ1 = float64(bigEndianToInt32(data[24:28])) / 0x40000000 + ap.QuatOrbitQ2 = float64(bigEndianToInt32(data[28:32])) / 0x40000000 + ap.QuatOrbitQ3 = float64(bigEndianToInt32(data[32:36])) / 0x40000000 + ap.QuatOrbJQ1 = float64(bigEndianToInt32(data[36:40])) / 0x40000000 + ap.QuatOrbJQ2 = float64(bigEndianToInt32(data[40:44])) / 0x40000000 + ap.QuatOrbJQ3 = float64(bigEndianToInt32(data[44:48])) / 0x40000000 + ap.Eular1 = float64(bigEndianToInt32(data[48:52])) / 1000000 + ap.Eular2 = float64(bigEndianToInt32(data[52:56])) / 1000000 + ap.Eular3 = float64(bigEndianToInt32(data[56:60])) / 1000000 + ap.DotEular1 = float64(bigEndianToInt16(data[60:62])) / 100000 + ap.DotEular2 = float64(bigEndianToInt16(data[62:64])) / 100000 + ap.DotEular3 = float64(bigEndianToInt16(data[64:66])) / 100000 + ap.ModTime = binary.BigEndian.Uint32(data[66:70]) + ap.DTime = float64(bigEndianToInt16(data[70:72])) / 1000 + ap.AutoState[0] = uint8(data[72]) + ap.AutoState[1] = uint8(data[73]) + ap.AutoState[2] = uint8(data[74]) + for i := 0; i < 16; i++ { + ap.ProTrack[i] = uint8(data[75+i]) + } + ap.QeQ1 = float64(bigEndianToInt32(data[91:95])) / 0x40000000 + ap.QeQ2 = float64(bigEndianToInt32(data[95:99])) / 0x40000000 + ap.QeQ3 = float64(bigEndianToInt32(data[99:103])) / 0x40000000 + ap.We1 = float64(bigEndianToInt16(data[103:105])) / 0x100000 + ap.We2 = float64(bigEndianToInt16(data[105:107])) / 0x100000 + ap.We3 = float64(bigEndianToInt16(data[107:109])) / 0x100000 + ap.WTFX = float64(bigEndianToInt16(data[109:111])) / 0x400000 + ap.WTFY = float64(bigEndianToInt16(data[111:113])) / 0x400000 + ap.WTFZ = float64(bigEndianToInt16(data[113:115])) / 0x400000 + ap.FbdriftX = float64(bigEndianToInt16(data[115:117])) / 0x4000000 + ap.FbdriftY = float64(bigEndianToInt16(data[117:119])) / 0x4000000 + ap.FbdriftZ = float64(bigEndianToInt16(data[119:121])) / 0x4000000 + ap.KalbY = float64(bigEndianToInt16(data[121:123])) / 0x4000000 + ap.KalbZ = float64(bigEndianToInt16(data[123:125])) / 0x4000000 + ap.KalbX = float64(bigEndianToInt16(data[125:127])) / 0x4000000 + ap.HTDX = float64(bigEndianToInt16(data[127:129])) / 0x200000 + ap.HTDY = float64(bigEndianToInt16(data[129:131])) / 0x200000 + ap.HTDZ = float64(bigEndianToInt16(data[131:133])) / 0x200000 + ap.CommandWheelX = float64(bigEndianToInt16(data[133:135])) / 0x40000 + ap.CommandWheelY = float64(bigEndianToInt16(data[135:137])) / 0x40000 + ap.CommandWheelZ = float64(bigEndianToInt16(data[137:139])) / 0x40000 + ap.QuatG1 = float64(bigEndianToInt32(data[139:143])) / 0x40000000 + ap.QuatG2 = float64(bigEndianToInt32(data[143:147])) / 0x40000000 + ap.QuatG3 = float64(bigEndianToInt32(data[147:151])) / 0x40000000 + ap.WG1 = float64(bigEndianToInt16(data[151:153])) / 0x40000 + ap.WG2 = float64(bigEndianToInt16(data[153:155])) / 0x40000 + ap.WG3 = float64(bigEndianToInt16(data[155:157])) / 0x40000 + ap.J2000PosX = float64(bigEndianToInt32(data[157:161])) / 0x100 + ap.J2000PosY = float64(bigEndianToInt32(data[161:165])) / 0x100 + ap.J2000PosZ = float64(bigEndianToInt32(data[165:169])) / 0x100 + ap.J2000VelX = float64(bigEndianToInt32(data[169:173])) / 0x100 + ap.J2000VelY = float64(bigEndianToInt32(data[173:177])) / 0x100 + ap.J2000VelZ = float64(bigEndianToInt32(data[177:181])) / 0x100 + ap.W84PosX = float64(bigEndianToInt32(data[181:185])) / 0x100 + ap.W84PosY = float64(bigEndianToInt32(data[185:189])) / 0x100 + ap.W84PosZ = float64(bigEndianToInt32(data[189:193])) / 0x100 + ap.W84VelX = float64(bigEndianToInt32(data[193:197])) / 0x100 + ap.W84VelY = float64(bigEndianToInt32(data[197:201])) / 0x100 + ap.W84VelZ = float64(bigEndianToInt32(data[201:205])) / 0x100 + ap.AngleDraft = float64(bigEndianToInt32(data[205:209])) / 10000000 + ap.DataTransLong = float64(bigEndianToInt16(data[209:211])) / 1000 + ap.DataTransLat = float64(bigEndianToInt16(data[211:213])) / 1000 + ap.DataTransH = float64(bigEndianToInt16(data[213:215])) + copy(ap.StarPrio[:], data[215:218]) + ap.WheelspeedcalX = float64(bigEndianToInt32(data[218:222])) / 10000 + ap.WheelspeedcalY = float64(bigEndianToInt32(data[222:226])) / 10000 + ap.WheelspeedcalZ = float64(bigEndianToInt32(data[226:230])) / 10000 + if ap.Waveway == 0x0 { + // 241.0 - 311.0 波道1参数 + } else { + // 241.1 - 311.1 波道2参数 + } + ap.WGS84PosX = float64(bigEndianToInt32(data[311:315])) / 100 + ap.WGS84PosY = float64(bigEndianToInt32(data[315:319])) / 100 + ap.WGS84PosZ = float64(bigEndianToInt32(data[319:323])) / 100 + ap.WGS84VelX = float64(bigEndianToInt32(data[323:327])) / 100 + ap.WGS84VelY = float64(bigEndianToInt32(data[327:331])) / 100 + ap.WGS84VelZ = float64(bigEndianToInt32(data[331:335])) / 100 + ap.J2000Pos_X = float64(bigEndianToInt32(data[335:339])) / 100 + ap.J2000Pos_Y = float64(bigEndianToInt32(data[339:343])) / 100 + ap.J2000Pos_Z = float64(bigEndianToInt32(data[343:347])) / 100 + ap.J2000Vel_X = float64(bigEndianToInt32(data[347:351])) / 100 + ap.J2000Vel_Y = float64(bigEndianToInt32(data[351:355])) / 100 + ap.J2000Vel_Z = float64(bigEndianToInt32(data[355:359])) / 100 + ap.FOGyroXAV = float64(bigEndianToInt24(data[359:362])) / 72359 + ap.FOGyroYAV = float64(bigEndianToInt24(data[362:365])) / 72359 + ap.FOGyroZAV = float64(bigEndianToInt24(data[365:368])) / 72359 + ap.FOGyroXAV1 = float64(bigEndianToInt24(data[368:371])) / 72359 + ap.FOGyroYAV1 = float64(bigEndianToInt24(data[371:374])) / 72359 + ap.FOGyroZAV1 = float64(bigEndianToInt24(data[374:377])) / 72359 + ap.SS1_UTCTime = binary.BigEndian.Uint32(data[377:381]) + ap.SS1_UTCTimeFrac = float32(uint32(data[381])<<16|uint32(data[382])<<8|uint32(data[383])) * 40.96 + ap.SS1_Q1 = float64(bigEndianToInt32(data[384:388])) / 2147483647 + ap.SS1_Q2 = float64(bigEndianToInt32(data[388:392])) / 2147483647 + ap.SS1_Q3 = float64(bigEndianToInt32(data[392:396])) / 2147483647 + ap.SS1_Q4 = float64(bigEndianToInt32(data[396:400])) / 2147483647 + ap.SS1_ExposureTime = uint8(data[400]) + ap.SS1_Threshold = uint8(data[401]) + ap.SS1_BackgroudV = uint8(data[402]) + ap.SS1_Flags = data[403] + ap.SS1_WorkMode = uint8(data[404] >> 6) + ap.SS1_ExtractStars = uint8(data[404] & 0x3F) + ap.SS1_NavStars = uint8(data[405] >> 5) + ap.SS1_Gain = uint8(data[405] & 0x1F) + ap.SS1_RegonizedStars = uint8(data[406] >> 2) + ap.SS1_ExtenalImage = data[406]&0x02 == 0x02 + ap.SS1_AttitudeActive = data[406]&0x01 == 0x01 + ap.SS1_ImgFrmNo = uint32(data[410])<<16 | uint32(data[411])<<8 | uint32(data[412]) + ap.SS1_XAV = float64(bigEndianToInt16(data[416:418])) * 2e-11 + ap.SS1_YAV = float64(bigEndianToInt16(data[418:420])) * 2e-11 + ap.SS1_ZAV = float64(bigEndianToInt16(data[420:422])) * 2e-11 + ap.SS2_UTCTime = binary.BigEndian.Uint32(data[423:427]) + ap.SS2_UTCTimeFrac = float32(uint32(data[427])<<16|uint32(data[428])<<8|uint32(data[429])) * 40.96 + ap.SS2_Q1 = float64(bigEndianToInt32(data[430:434])) / 2147483647 + ap.SS2_Q2 = float64(bigEndianToInt32(data[434:438])) / 2147483647 + ap.SS2_Q3 = float64(bigEndianToInt32(data[438:442])) / 2147483647 + ap.SS2_Q4 = float64(bigEndianToInt32(data[442:446])) / 2147483647 + ap.SS2_ExposureTime = uint8(data[446]) + ap.SS2_Threshold = uint8(data[447]) + ap.SS2_BackgroudV = uint8(data[448]) + ap.SS2_Flags = data[449] + ap.SS2_WorkMode = uint8(data[450] >> 6) + ap.SS2_ExtractStars = uint8(data[450] & 0x3F) + ap.SS2_NavStars = uint8(data[451] >> 5) + ap.SS2_Gain = uint8(data[451] & 0x1F) + ap.SS2_RegonizedStars = uint8(data[452] >> 2) + ap.SS2_ExtenalImage = data[452]&0x02 == 0x02 + ap.SS2_AttitudeActive = data[452]&0x01 == 0x01 + ap.SS2_ImgFrmNo = uint32(data[456])<<16 | uint32(data[457])<<8 | uint32(data[458]) + ap.FlyWheel1_Vel = float64(bigEndianToInt16(data[463:466])) * 0.1 + ap.FlyWheel1_Amps = float64(uint8(data[466])) * 0.009 + ap.FlyWheel2_Vel = float64(bigEndianToInt16(data[467:470])) * 0.1 + ap.FlyWheel2_Amps = float64(uint8(data[470])) * 0.009 + ap.FlyWheel3_Vel = float64(bigEndianToInt16(data[471:474])) * 0.1 + ap.FlyWheel3_Amps = float64(uint8(data[474])) * 0.009 + ap.FlyWheel4_Vel = float64(bigEndianToInt16(data[475:478])) * 0.1 + ap.FlyWheel4_Amps = float64(uint8(data[478])) * 0.009 + ap.MemsGyroXAV = float64(bigEndianToInt32(data[479:483])) * 0.000001 + ap.MemsGyroYAV = float64(bigEndianToInt32(data[483:487])) * 0.000001 + ap.MemsGyroZAV = float64(bigEndianToInt32(data[487:491])) * 0.000001 + ap.MagnetictrengthX = float64(bigEndianToInt32(data[491:495])) / 19200 + ap.MagnetictrengthY = float64(bigEndianToInt32(data[495:499])) / 19200 + ap.MagnetictrengthZ = float64(bigEndianToInt32(data[499:503])) / 19200 + ap.ASSTimeInt = binary.BigEndian.Uint32(data[503:507]) + ap.ASSTimeDec = uint32(data[507])<<16 | uint32(data[508])<<8 | uint32(data[509]) + ap.CheckSum = data[511] + + return nil +} + +func (ap AuxPlatform) Print() { + pp.Println(ap) +} + +// Extractor 辅助数据提取器 +// auxfile 辅助数据文件路径 +func ParseAuxPlatform(auxfile string) ([]*AuxPlatform, error) { + data, err := os.ReadFile(auxfile) + if err != nil { + log.Println("read aux data from", auxfile, "error:", err.Error()) + return nil, err + } + + var aps []*AuxPlatform + for i := 0; i < len(data); i += 24 + 128 + 512 { + ap := AuxPlatform{} + if err := ap.Decode(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.WGS84XYZtoLatLngH(ap.W84PosX, ap.W84PosY, ap.W84PosZ)) + pp.Println(calculator.WGS84XYZtoLatLngH(ap.WGS84PosX, ap.WGS84PosY, ap.WGS84PosZ)) + + pp.Println(calculator.J2000ToWGS84( + ap.J2000PosX, + ap.J2000PosY, + ap.J2000PosZ, + time.Unix( + int64(ap.UTCTimeSec+uint32(ReferenceTime2000)), + int64(ap.Microsecond)*1000).UTC())) + + pp.Println(calculator.J2000ToWGS84( + ap.J2000Pos_X, + ap.J2000Pos_Y, + ap.J2000Pos_Z, + time.Unix( + int64(ap.UTCTimeSec+uint32(ReferenceTime2000)), + int64(ap.Microsecond)*1000).UTC())) + break + } + return aps, nil +} + +func Time2000UTCSec() int64 { + t, _ := time.ParseInLocation("2006-01-02 15:04:05", "2000-01-01 12:00:00", time.UTC) + return t.Unix() +} diff --git a/pkg/auxilary/util.go b/pkg/auxilary/util.go new file mode 100644 index 0000000..37b5ad2 --- /dev/null +++ b/pkg/auxilary/util.go @@ -0,0 +1,51 @@ +package auxilary + +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 +} diff --git a/pkg/calculator/const.go b/pkg/calculator/const.go new file mode 100644 index 0000000..05f92d7 --- /dev/null +++ b/pkg/calculator/const.go @@ -0,0 +1,10 @@ +package calculator + +// WGS84 ellipsoid constants +const ( + EarthRadius = 6378137.0 // 地球半径，单位米 + a = 6378137.0 // semi-major axis in meters + f = 1 / 298.257223563 // flattening + e2 = 2*f - f*f // square of eccentricity + J2000Epoch = 2451545.0 // Julian date of J2000 epoch +) diff --git a/pkg/calculator/j2000.go b/pkg/calculator/j2000.go new file mode 100644 index 0000000..fac7b8c --- /dev/null +++ b/pkg/calculator/j2000.go @@ -0,0 +1,78 @@ +package calculator + +import ( + "math" + "time" +) + +// Function to convert current J2000 position to WGS84 +func J2000ToWGS84(j2000X, j2000Y, j2000Z float64, utc time.Time) (float64, float64, float64) { + julianDate := UTCToJulianDate(utc) + gast := CalculateGAST(julianDate, utc) + + // Calculate rotation matrix from J2000 to ITRS + rotationMatrix := CreateRotationMatrix(gast) + + // Rotate J2000 position to ITRS + itrsX := rotationMatrix[0][0]*j2000X + rotationMatrix[0][1]*j2000Y + rotationMatrix[0][2]*j2000Z + itrsY := rotationMatrix[1][0]*j2000X + rotationMatrix[1][1]*j2000Y + rotationMatrix[1][2]*j2000Z + itrsZ := rotationMatrix[2][0]*j2000X + rotationMatrix[2][1]*j2000Y + rotationMatrix[2][2]*j2000Z + + // Convert ITRS to geodetic coordinates (WGS84) + latitude, longitude, height := ECEFToGeodetic(itrsX, itrsY, itrsZ) + return latitude, longitude, height +} + +// Function to convert UTC to Julian Date +func UTCToJulianDate(t time.Time) float64 { + year, month, day := t.Date() + hour, minute, second := t.Clock() + fracOfDay := float64(hour)/24 + float64(minute)/(24*60) + float64(second)/(24*60*60) + + if month <= 2 { + year -= 1 + month += 12 + } + + A := year / 100 + B := 2 - A + A/4 + julianDate := float64(int(365.25*float64(year))) + float64(int(30.6001*float64(month+1))) + float64(day) + 1720994.5 + fracOfDay + float64(B) + return julianDate +} + +// Function to calculate GAST (Greenwich Apparent Sidereal Time) +func CalculateGAST(julianDate float64, utc time.Time) float64 { + // Calculate Julian centuries since J2000.0 + T := (julianDate - J2000Epoch) / 36525.0 + + // Calculate Greenwich Mean Sidereal Time (GMST) in degrees + GMST := 280.46061837 + 360.98564736629*(julianDate-J2000Epoch) + 0.000387933*T*T - T*T*T/38710000 + GMST = math.Mod(GMST, 360.0) + if GMST < 0 { + GMST += 360.0 + } + + // Calculate the equation of the equinoxes (simplified) + omega := 125.04 - 0.052954*T + L := 280.47 + 0.98565*T + epsilon := 23.4393 - 0.0000004*T + deltaPsi := -0.000319*math.Sin(math.Pi/180.0*omega) - 0.000024*math.Sin(math.Pi/180.0*2*L) + + // Greenwich Apparent Sidereal Time (GAST) + GAST := GMST + deltaPsi*math.Cos(math.Pi/180.0*epsilon) + GAST = math.Mod(GAST, 360.0) + if GAST < 0 { + GAST += 360.0 + } + + return GAST * math.Pi / 180.0 // Convert to radians +} + +// Function to create the rotation matrix for the given GAST +func CreateRotationMatrix(gast float64) [3][3]float64 { + return [3][3]float64{ + {math.Cos(gast), math.Sin(gast), 0}, + {-math.Sin(gast), math.Cos(gast), 0}, + {0, 0, 1}, + } +} diff --git a/pkg/calculator/quaternion.go b/pkg/calculator/quaternion.go new file mode 100644 index 0000000..e19b58b --- /dev/null +++ b/pkg/calculator/quaternion.go @@ -0,0 +1,70 @@ +package calculator + +import ( + "fmt" + "math" +) + +// Quaternion represents a quaternion with scalar (w) and vector (x, y, z) parts +type Quaternion struct { + w, x, y, z float64 +} + +// Quaternion multiplication +func (q1 Quaternion) Mul(q2 Quaternion) Quaternion { + return Quaternion{ + w: q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z, + x: q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y, + y: q1.w*q2.y - q1.x*q2.z + q1.y*q2.w + q1.z*q2.x, + z: q1.w*q2.z + q1.x*q2.y - q1.y*q2.x + q1.z*q2.w, + } +} + +// Quaternion conjugate +func (q Quaternion) Conjugate() Quaternion { + return Quaternion{w: q.w, x: -q.x, y: -q.y, z: -q.z} +} + +// Rotate vector by quaternion +func (q Quaternion) Rotate(v [3]float64) [3]float64 { + qv := Quaternion{w: 0, x: v[0], y: v[1], z: v[2]} + qConj := q.Conjugate() + qvRotated := q.Mul(qv).Mul(qConj) + return [3]float64{qvRotated.x, qvRotated.y, qvRotated.z} +} + +func main() { + // 示例数据 + qBI := Quaternion{w: 1, x: 0, y: 0, z: 0} // 本体相对惯性系四元数 + posJ2000 := [3]float64{7000, 0, 0} // J2000位置 + // velJ2000 := [3]float64{0, 7.5, 0} // J2000速度 + + // 相机参数 + const numPixels = 9520 + const fov = 10.0 * math.Pi / 180 // 假设视场角为10度 + + // 逐像素计算地面交点 + for i := 0; i < numPixels; i++ { + // 计算像素点相对光轴的偏角 + pixelOffset := (float64(i) - float64(numPixels)/2) / float64(numPixels) + angle := pixelOffset * fov + + // 假设光轴在本体坐标系中指向-z方向，计算视线方向 + dBody := [3]float64{-math.Sin(angle), 0, -math.Cos(angle)} + + // 转换到惯性系 + dInertial := qBI.Rotate(dBody) + + // 计算地面交点（假设dInertial已经标准化） + k := -posJ2000[2] / dInertial[2] // 简化的交点计算 + groundPoint := [3]float64{ + posJ2000[0] + k*dInertial[0], + posJ2000[1] + k*dInertial[1], + posJ2000[2] + k*dInertial[2], + } + + // 转换到地理坐标 + lat, lon, _ := ECEFToGeodetic(groundPoint[0], groundPoint[1], groundPoint[2]) + fmt.Printf("Pixel %d: Latitude: %f, Longitude: %f\n", i, lat, lon) + } +} diff --git a/pkg/calculator/wgs84.go b/pkg/calculator/wgs84.go new file mode 100644 index 0000000..dcaa793 --- /dev/null +++ b/pkg/calculator/wgs84.go @@ -0,0 +1,77 @@ +package calculator + +import "math" + +func WGS84XYZtoLatLngH(X, Y, Z float64) (float64, float64, float64) { + return ECEFToGeodetic(X, Y, Z) +} + +// Function to convert ECEF (ITRS) coordinates to geodetic coordinates (latitude, longitude, height) +func ECEFToGeodetic(X, Y, Z float64) (float64, float64, float64) { + b := a * (1 - f) + e2Prime := e2 * (a * a) / (b * b) + p := math.Sqrt(X*X + Y*Y) + theta := math.Atan2(Z*a, p*b) + + // Calculate Longitude + longitude := math.Atan2(Y, X) + + // Calculate Latitude + latitude := math.Atan2(Z+e2Prime*b*math.Pow(math.Sin(theta), 3), p-e2*a*math.Pow(math.Cos(theta), 3)) + + // Calculate Height + N := a / math.Sqrt(1-e2*math.Pow(math.Sin(latitude), 2)) + height := p/math.Cos(latitude) - N + + // Convert radians to degrees for latitude and longitude + latitudeDeg := latitude * 180 / math.Pi + longitudeDeg := longitude * 180 / math.Pi + + return latitudeDeg, longitudeDeg, height +} + +// Converts degrees to radians. +func degreesToRadians(degrees float64) float64 { + return degrees * math.Pi / 180 +} + +// Converts radians to degrees. +func radiansToDegrees(radians float64) float64 { + return radians * 180 / math.Pi +} + +// Calculate destination point given a start point (lat1, lng1), distance dx, dy (in meters). +func CalculateDestination(lat1, lng1, dx, dy float64) (float64, float64) { + // Convert latitude and longitude from degrees to radians. + lat1Rad := degreesToRadians(lat1) + lng1Rad := degreesToRadians(lng1) + + // Radius of the Earth at the given latitude. + radius := EarthRadius * math.Cos(lat1Rad) + + // Calculate the new latitude in radians. + newLatRad := lat1Rad + dy/EarthRadius + + // Calculate the new longitude in radians. + newLngRad := lng1Rad + dx/radius + + // Convert the new latitude and longitude back to degrees. + newLat := radiansToDegrees(newLatRad) + newLng := radiansToDegrees(newLngRad) + + return newLat, newLng +} + +const ( + MetersPerDegreeLatitude = 111320.0 // 1度纬度约等于111,320米 +) + +// Converts a distance in meters to degrees longitude at a specific latitude +func MetersToDegreesLongitude(meters float64, latitude float64) float64 { + return meters / (111320.0 * cosLatitude(latitude)) +} + +// Approximates the cosine of the latitude in radians +func cosLatitude(latitude float64) float64 { + return 1.0 / (1.0 + (latitude/90.0)*(latitude/90.0)) +} diff --git a/pkg/config/config.go b/pkg/config/config.go index 19fda63..6a1406a 100644 --- a/pkg/config/config.go +++ b/pkg/config/config.go @@ -1,7 +1,10 @@ package config +import "github.com/sirupsen/logrus" + type Config struct { CRConfig CoRegistrationConfig `yaml:"cr_config" mapstructure:"cr_config"` + LogLevel logrus.Level `yaml:"log_level" mapstructure:"log_level"` } type CoRegistrationConfig struct { @@ -20,6 +23,7 @@ var GCONFIG Config func init() { GCONFIG = Config{ + LogLevel: logrus.DebugLevel, CRConfig: CoRegistrationConfig{ MssBands: 4, PixelBytes: 2, diff --git a/pkg/config/viper.go b/pkg/config/viper.go index 2ce094c..8a017b1 100644 --- a/pkg/config/viper.go +++ b/pkg/config/viper.go @@ -10,7 +10,7 @@ import ( "github.com/spf13/viper" ) -var viperInst *viper.Viper +var GViper *viper.Viper func InitViper(cfg string) *viper.Viper { v := viper.New() @@ -22,7 +22,7 @@ func InitViper(cfg string) *viper.Viper { v.SetConfigType("yaml") err := v.ReadInConfig() if err != nil { - log.Println(fmt.Errorf("fatal error config file: %s", err.Error())) + log.Println(fmt.Errorf("cannot read config file: %s, use default config", err.Error())) return nil } diff --git a/pkg/producer/aux.go b/pkg/producer/aux.go new file mode 100644 index 0000000..8b422d1 --- /dev/null +++ b/pkg/producer/aux.go @@ -0,0 +1,109 @@ +package imageproc + +import ( + "math" + "time" + + "github.com/paulmach/orb" + "github.com/paulmach/orb/planar" + "starwiz.cn/sjy01/image-proc/pkg/auxilary" + "starwiz.cn/sjy01/image-proc/pkg/calculator" +) + +func (r *Registrator) LoadAuxData() error { + var err error + r.auxHeads, r.auxBoxes, r.AuxPlatforms, err = auxilary.ExtractAux(r.Params.AuxRawFile) + return err +} + +func (r *Registrator) SceneImageTime(scene *Scene) (start, center, end time.Time) { + startPosInAux, endPosInAux := r.SceneInAuxIndex(scene) + centerPosInAux := (startPosInAux + endPosInAux) / 2 + + start = time.Unix(int64(auxilary.ReferenceTime2000)+int64(r.AuxPlatforms[startPosInAux].UTCTimeSec), + int64(r.AuxPlatforms[startPosInAux].Microsecond)*1000) + center = time.Unix(int64(auxilary.ReferenceTime2000)+int64(r.AuxPlatforms[centerPosInAux].UTCTimeSec), + int64(r.AuxPlatforms[centerPosInAux].Microsecond)*1000) + end = time.Unix(int64(auxilary.ReferenceTime2000)+int64(r.AuxPlatforms[endPosInAux].UTCTimeSec), + int64(r.AuxPlatforms[endPosInAux].Microsecond)*1000) + + return +} + +// FIXME: 位置像元经纬度计算方法，暂时使用星下点替代 +func (r *Registrator) ScenePosition(scene *Scene) (topLeft, bottomRight orb.Point) { + // startPosInAux, endPosInAux := r.SceneInAuxIndex(scene) + ap := r.AuxPlatforms[0] + // ap1 := r.AuxPlatforms[endPosInAux] + lat0, lng0, _ := calculator.WGS84XYZtoLatLngH(ap.WGS84PosX, ap.WGS84PosY, ap.WGS84PosZ) + + lat, lng := calculator.CalculateDestination(lat0, lng0, + float64(scene.Width)*scene.Meta.Gsd, float64(-scene.Y)*scene.Meta.Gsd) + + // lat0, lng0, h0 := calculator.WGS84XYZtoLatLngH(ap1.WGS84PosX, ap1.WGS84PosY, ap1.WGS84PosZ) + // fmt.Println("Scene Position End:", lat0, lng0, h0) + + lat1, lng1 := calculator.CalculateDestination(lat, lng, + float64(scene.Width)*scene.Meta.Gsd, float64(-scene.Height)*scene.Meta.Gsd) + + poly := orb.Polygon{ + { + {lng, lat}, + {lng1, lat}, + {lng1, lat1}, + {lng, lat1}, + {lng, lat}, + }, + } + + centroid, _ := planar.CentroidArea(poly) + scene.Meta.CentreLocation.Latitude = centroid.Y() + scene.Meta.CentreLocation.Longitude = centroid.X() + scene.Meta.Corners.UpperLeft.Latitude = lat + scene.Meta.Corners.UpperLeft.Longitude = lng + scene.Meta.Corners.UpperRight.Latitude = lat + scene.Meta.Corners.UpperRight.Longitude = lng1 + scene.Meta.Corners.LowerLeft.Latitude = lat1 + scene.Meta.Corners.LowerLeft.Longitude = lng + scene.Meta.Corners.LowerRight.Latitude = lat1 + scene.Meta.Corners.LowerRight.Longitude = lng1 + + scene.Meta.SatPosX = ap.WGS84PosX + scene.Meta.SatPosY = ap.WGS84PosY + scene.Meta.SatPosZ = ap.WGS84PosZ + scene.Meta.Yaw = ap.Eular1 * 180 / math.Pi + scene.Meta.Pitch = ap.Eular2 * 180 / math.Pi + scene.Meta.Roll = ap.Eular3 * 180 / math.Pi + + // feature := geojson.NewFeature(poly) + // fcs.Features = append(fcs.Features, feature) + + // fd, _ := fcs.MarshalJSON() + // fmt.Println(string(fd)) + + return +} + +func (r *Registrator) SceneInAuxIndex(scene *Scene) (int, int) { + var auxForImageRow int + switch scene.Type { + case "MSS": + auxForImageRow = 4 + case "PAN": + auxForImageRow = 16 + case "FUS": + auxForImageRow = 16 + } + + startPosInAux := scene.Y / auxForImageRow + if startPosInAux >= len(r.AuxPlatforms) { + startPosInAux = len(r.AuxPlatforms) - 1 + } + + endPosInAux := (scene.Y + scene.Height) / auxForImageRow + if endPosInAux >= len(r.AuxPlatforms) { + endPosInAux = len(r.AuxPlatforms) - 1 + } + + return startPosInAux, endPosInAux +} diff --git a/pkg/producer/image_registration.go b/pkg/producer/image_registration.go index 96d0a0c..db8cfbb 100644 --- a/pkg/producer/image_registration.go +++ b/pkg/producer/image_registration.go @@ -11,6 +11,7 @@ import ( "github.com/airbusgeo/godal" log "github.com/sirupsen/logrus" "gocv.io/x/gocv" + "starwiz.cn/sjy01/image-proc/pkg/auxilary" ) type Registrate interface{} @@ -49,6 +50,10 @@ type Registrator struct { rgbirImage gocv.Mat resampleMethod ResampleMethod + + auxHeads []*auxilary.AuxFrameHead + auxBoxes []*auxilary.AuxFocalBox + AuxPlatforms []*auxilary.AuxPlatform } func NewRegistrator(rsmethod ResampleMethod) *Registrator { diff --git a/pkg/producer/jpg.go b/pkg/producer/jpg.go index 6693c3e..0106086 100644 --- a/pkg/producer/jpg.go +++ b/pkg/producer/jpg.go @@ -59,9 +59,7 @@ func GTiffToJPG(ftiff, fjpg string, reversed bool) error { channels := gocv.Split(img) for i, ch := range channels { // 2. 计算图像的最小值和最大值 - minVal, maxVal, minLoc, maxLoc := gocv.MinMaxLoc(ch) - log.Printf("Min value: %f, Max value: %f, min location: %v, max location: %v", minVal, maxVal, minLoc, maxLoc) - + minVal, maxVal, _, _ := gocv.MinMaxLoc(ch) // 3. 将16位图像线性拉伸到8位图像 scale := 255.0 / (maxVal - minVal) shift := -minVal * scale @@ -85,10 +83,6 @@ func GTiffToJPG(ftiff, fjpg string, reversed bool) error { 0, 0, gocv.InterpolationCubic) // 7. 应用伽玛校正提升亮度 - // avgBrightness := calculateAverageBrightness(img8bit) - // gamma := determineGammaValue(avgBrightness) - // log.Printf("Average Brightness: %f, Determined Gamma: %f", avgBrightness, gamma) - gammaCorrected := applyGammaCorrection(img8bit, 1.6) // 伽玛值 defer gammaCorrected.Close() ok := gocv.IMWriteWithParams(fjpg, gammaCorrected, []int{gocv.IMWriteJpegOptimize, 1}) @@ -96,9 +90,6 @@ func GTiffToJPG(ftiff, fjpg string, reversed bool) error { err = fmt.Errorf("error saving %s", fjpg) return err } - - // gocv.IMWrite(strings.Replace(fjpg, ".jpg", "_8.jpg", 1), img8bit) - return nil } @@ -121,24 +112,3 @@ func applyGammaCorrection(src gocv.Mat, gamma float64) gocv.Mat { gocv.LUT(src, gammaMat, &dst) return dst } - -// calculateAverageBrightness calculates the average brightness of the image -func calculateAverageBrightness(img gocv.Mat) float64 { - if img.Channels() == 1 { - return img.Mean().Val1 - } - gray := gocv.NewMat() - defer gray.Close() - gocv.CvtColor(img, &gray, gocv.ColorBGRToGray) - return gray.Mean().Val1 -} - -// determineGammaValue determines an appropriate gamma value based on the average brightness -func determineGammaValue(avgBrightness float64) float64 { - targetBrightness := 180.0 - if avgBrightness < 0.1 { - return 1.0 - } - gamma := math.Log(targetBrightness/255.0) / math.Log(avgBrightness/255.0) - return gamma -} diff --git a/pkg/producer/meta.go b/pkg/producer/meta.go index 8b29c10..4435453 100644 --- a/pkg/producer/meta.go +++ b/pkg/producer/meta.go @@ -3,6 +3,9 @@ package imageproc import ( "encoding/xml" "os" + "path/filepath" + "strings" + "time" ) // 定义与XML结构对应的Go结构体 @@ -55,7 +58,38 @@ type Corners struct { LowerLeft Location `xml:"LowerLeft"` } -func WriteProductMeta(productMeta *ProductMeta, filename string) error { +func (r *Registrator) makeProductMeta(scene *Scene) *ProductMeta { + meta := &ProductMeta{ + Satellite: "SJY01", + Sensor: "PMS", + OutputFormat: "GTiff", + ProductGenTime: time.Now().Format("2006-01-02T15:04:05"), + Width: scene.Width, + Height: scene.Height, + MapProjection: "GEOGRAPHIC", + EarthEllipsoid: "WGS_84", + ProductLevel: "L1A", + } + + switch scene.Type { + case "PAN": + meta.Gsd = 1.25 + meta.Bands = 1 + case "MSS": + meta.Gsd = 5 + meta.Bands = 4 + } + + meta.ProductID = filepath.Base(strings.TrimSuffix(scene.Tiff, ".tiff")) + startTime, centerTime, endTime := r.SceneImageTime(scene) + meta.StartTime = startTime.Format("2006-01-02T15:04:05.000") + meta.EndTime = endTime.Format("2006-01-02T15:04:05.000") + meta.CentreTime = centerTime.Format("2006-01-02T15:04:05.000") + + return meta +} + +func (r *Registrator) writeProductMeta(productMeta *ProductMeta, filename string) error { output, err := xml.MarshalIndent(productMeta, "", " ") if err != nil { return err diff --git a/pkg/producer/output.go b/pkg/producer/output.go index 77c7c97..1b733e4 100644 --- a/pkg/producer/output.go +++ b/pkg/producer/output.go @@ -13,7 +13,7 @@ import ( ) func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error { - err := savePanToGDALGTiff(r.PanImage, tiffFile) + err := savePanToGDALGTiff(r.PanImage, 0, 0, tiffFile) if err != nil { return err } @@ -21,13 +21,7 @@ func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error { return nil } -func (r *Registrator) SaveFilteredPanToGDALGTiff(tiffFile string) error { - img := PANFilter(r.PanImage) - img.ConvertTo(&img, gocv.MatTypeCV16U) - return savePanToGDALGTiff(img, tiffFile) -} - -func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error { +func savePanToGDALGTiff(pan gocv.Mat, topLeftX, topLeftY float64, tiffFile string) error { // log.Println("Saving PAN image to TIFF file:", tiffFile) width := pan.Cols() @@ -40,7 +34,7 @@ func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error { } defer ds.Close() - setGeoTransform(ds, 0, 0, float64(width), float64(height), 1.25) + setGeoTransform(ds, topLeftX, topLeftY, 1.25) ds.SetMetadata("NBITS", "16") // 将通道的数据转换为uint16数组 @@ -74,7 +68,7 @@ func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error { gocv.Merge(r.registeredMssImages[:], &r.rgbirImage) err := SaveBGRToGDALGTiff(r.rgbirImage, - 4, 5, + 4, 0, 0, 5, []godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined}, tiffFile) if err != nil { @@ -84,15 +78,11 @@ func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error { return nil } -func (r *Registrator) SavePansharpenedToGDALGTiff(tiffFile string) error { - ihsImage := PansharpenIHS(r.PanImage, r.rgbirImage) - return SaveBGRToGDALGTiff(ihsImage, - 3, 1.25, - []godal.ColorInterp{godal.CIRed, godal.CIGreen, godal.CIBlue}, - tiffFile) -} - -func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, colorInterps []godal.ColorInterp, tiffFile string) error { +func SaveBGRToGDALGTiff(bgr gocv.Mat, + bands int, + topLeftX, topLeftY float64, + resolution float64, + colorInterps []godal.ColorInterp, tiffFile string) error { width := bgr.Cols() height := bgr.Rows() @@ -123,8 +113,7 @@ func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, colorInterp defer ds.Close() // ds.SetMetadata("NBITS", "16") - - setGeoTransform(ds, 0, 0, float64(width), float64(height), resolution) + setGeoTransform(ds, topLeftX, topLeftY, resolution) for b := 0; b < bands; b++ { band := ds.Bands()[b] diff --git a/pkg/producer/proj.go b/pkg/producer/proj.go index aef3b9f..ff3104f 100644 --- a/pkg/producer/proj.go +++ b/pkg/producer/proj.go @@ -75,6 +75,8 @@ func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"], AUTHORITY["EPSG","3857"]]`) + projections = append(projections, `GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]`) + // 输出投影信息 - return projections[0] + return projections[4] } diff --git a/pkg/producer/scenes.go b/pkg/producer/scenes.go index cacf087..0414483 100644 --- a/pkg/producer/scenes.go +++ b/pkg/producer/scenes.go @@ -13,12 +13,15 @@ import ( ) type Scene struct { - Width int - Height int - X int // coordinate of the left top corner of the scene - Y int - Mat []gocv.Mat - Tiff string + Type string + Width int + Height int + X int // coordinate of the left top corner of the scene + Y int + Mat []gocv.Mat + Tiff string + Meta *ProductMeta + SceneId string } func (s *Scene) Cleanup() { @@ -43,12 +46,17 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e } scene := &Scene{ + Type: "PAN", Width: 9344, Height: y1 - i*hPAN, X: 0, Y: i * hPAN, } + name := filepath.Base(r.Params.PanTiffFile) + name = strings.TrimSuffix(name, ".tiff") + scene.SceneId = fmt.Sprintf("%s_%03d", name, i+1) + mat := r.PanImage.Region(image.Rect(0, i*hPAN, 9344, y1)) scene.Mat = append(scene.Mat, mat) panScenes = append(panScenes, scene) @@ -64,6 +72,7 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e } scene := &Scene{ + Type: "MSS", Width: 2336, Height: y1 - i*hMSS, X: 0, @@ -75,6 +84,9 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e scene.Mat = append(scene.Mat, mat) } + name := filepath.Base(r.Params.MssTiffFile) + name = strings.TrimSuffix(name, ".tiff") + scene.SceneId = fmt.Sprintf("%s_%03d", name, i+1) mssScenes = append(mssScenes, scene) } @@ -91,16 +103,23 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e dir := filepath.Join(r.Params.OutputDir, fmt.Sprintf("%03d", i+1), "PAN") os.MkdirAll(dir, 0755) - tiff := filepath.Base(r.Params.PanTiffFile) - tiff = strings.TrimSuffix(tiff, ".tiff") - filename := fmt.Sprintf("%s_%03d.tiff", tiff, i+1) + filename := fmt.Sprintf("%s_L1A.tiff", scene.SceneId) scene.Tiff = filepath.Join(dir, filename) - err := savePanToGDALGTiff(scene.Mat[0], scene.Tiff) + scene.Meta = r.makeProductMeta(scene) + r.ScenePosition(scene) + + err := savePanToGDALGTiff(scene.Mat[0], + scene.Meta.Corners.UpperLeft.Longitude, + scene.Meta.Corners.UpperLeft.Latitude, + scene.Tiff) if err != nil { log.Errorf("save pan scene %d to tiff failed: %v", i+1, err) return err } + scene.Meta.DataSize = sizeOfFile(scene.Tiff) + metaFile := strings.Replace(scene.Tiff, ".tiff", ".meta.xml", 1) + r.writeProductMeta(scene.Meta, metaFile) GTiffToJPG(scene.Tiff, strings.Replace(scene.Tiff, ".tiff", ".jpg", 1), false) } @@ -108,14 +127,17 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e dir := filepath.Join(r.Params.OutputDir, fmt.Sprintf("%03d", i+1), "MSS") os.MkdirAll(dir, 0755) - tiff := filepath.Base(r.Params.MssTiffFile) - tiff = strings.TrimSuffix(tiff, ".tiff") - filename := fmt.Sprintf("%s_%03d.tiff", tiff, i+1) + filename := fmt.Sprintf("%s_L1A.tiff", scene.SceneId) scene.Tiff = filepath.Join(dir, filename) + scene.Meta = r.makeProductMeta(scene) + r.ScenePosition(scene) rgbirImage, _ := r.MergeMSSToBGRNIR(scene.Mat) err := SaveBGRToGDALGTiff(rgbirImage, - 4, 5, + 4, + scene.Meta.Corners.UpperLeft.Longitude, + scene.Meta.Corners.UpperLeft.Latitude, + 5, []godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined}, scene.Tiff) if err != nil { @@ -123,6 +145,9 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e return err } + scene.Meta.DataSize = sizeOfFile(scene.Tiff) + metaFile := strings.Replace(scene.Tiff, ".tiff", ".meta.xml", 1) + r.writeProductMeta(scene.Meta, metaFile) GTiffToJPG(scene.Tiff, strings.Replace(scene.Tiff, ".tiff", ".jpg", 1), false) } diff --git a/pkg/producer/util.go b/pkg/producer/util.go index bacc338..f5e25bc 100644 --- a/pkg/producer/util.go +++ b/pkg/producer/util.go @@ -1,29 +1,47 @@ package imageproc import ( + "os" "time" "github.com/airbusgeo/godal" log "github.com/sirupsen/logrus" + "starwiz.cn/sjy01/image-proc/pkg/calculator" ) -func setGeoTransform(ds *godal.Dataset, topLeftX, topLeftY, width, height, resolution float64) { +func setGeoTransform(ds *godal.Dataset, topLeftLng, topLeftLat, resolution float64) { + // 转换分辨率（米到度） + resLat := resolution / calculator.MetersPerDegreeLatitude + resLng := calculator.MetersToDegreesLongitude(resolution, topLeftLat) // 设置地理变换（假设左上角坐标为(0,0)，PAN每个像素分辨率为1.2米） geotransform := [6]float64{ - 0, // top left x - resolution, // w-e pixel resolution - 0, // rotation, 0 if image is "north up" - 0, // top left y - 0, // rotation, 0 if image is "north up" - -resolution, // n-s pixel resolution (negative value) + topLeftLng, // top left x + resLng, // w-e pixel resolution + 0, // rotation, 0 if image is "north up" + topLeftLat, // top left y + 0, // rotation, 0 if image is "north up" + -resLat, // n-s pixel resolution (negative value) } err := ds.SetGeoTransform(geotransform) if err != nil { log.Errorf("Failed to set GeoTransform: %v", err) + return } - // 设置投影信息（WGS84） - err = ds.SetProjection(calculateProj(topLeftX, topLeftY, width, height, resolution)) + // 设置投影为 WGS84 + srs, err := godal.NewSpatialRefFromEPSG(4326) // gdal.CreateSpatialReference("") + if err != nil { + log.Errorf("Failed to set spatial reference: %v", err) + return + } + + projWKT, err := srs.WKT() + if err != nil { + log.Errorf("Failed to convert spatial reference to WKT: %v", err) + return + } + + err = ds.SetProjection(projWKT) if err != nil { log.Errorf("Failed to set projection: %v", err) } @@ -32,3 +50,12 @@ func setGeoTransform(ds *godal.Dataset, topLeftX, topLeftY, width, height, resol ds.SetMetadata("TIFFTAG_DATETIME", time.Now().String()) ds.SetMetadata("TIFFTAG_SOFTWARE", "StarWiz-SJY01-IMAGE-PROC") } + +func sizeOfFile(file string) int64 { + fileInfo, err := os.Stat(file) + if err != nil { + log.Errorf("Failed to get file size: %v", err) + return 0 + } + return fileInfo.Size() +}