diff --git a/cmd/proc.go b/cmd/proc.go index df56fb8..0938f45 100644 --- a/cmd/proc.go +++ b/cmd/proc.go @@ -90,7 +90,7 @@ var procCmd = &cobra.Command{ logrus.Error(err) } - reg.SaveScenesToTiff(panScenes, mssScenes) + reg.OutputL1A(panScenes, mssScenes) producer.CleanScenes(panScenes) producer.CleanScenes(mssScenes) diff --git a/cmd/warp.go b/cmd/warp.go index a4ca6bb..76ac864 100644 --- a/cmd/warp.go +++ b/cmd/warp.go @@ -22,15 +22,15 @@ var ( var cmdWarp = &cobra.Command{ Use: "warp", - Short: "transform the input image using rpc", + Short: "L1A to L2: transform the input image using rpc", Run: func(cmd *cobra.Command, args []string) { godal.RegisterAll() dem.Dem1KmLT = dem.NewDem1Km(config.GCONFIG.Dem.Dem1Km) warpDEM, _ = filepath.Abs(config.GCONFIG.Dem.Dem1Km) - err := producer.Warp(warpIn, warpOut, warpMeta, warpRPC, warpDEM) + err := producer.L1AtoL2(warpIn, warpOut, warpMeta, warpRPC, warpDEM) if err != nil { logrus.Error("Generate L2 product failed:", err) - os.Exit(10) + os.Exit(1) } }, } diff --git a/pkg/producer/meta.go b/pkg/producer/meta.go index 16f3af5..4cb70d1 100644 --- a/pkg/producer/meta.go +++ b/pkg/producer/meta.go @@ -2,6 +2,7 @@ package producer import ( "encoding/xml" + "math" "os" "path/filepath" "strings" @@ -58,6 +59,26 @@ type Corners struct { LowerLeft Location `xml:"LowerLeft"` } +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) + + return +} + func (r *Registrator) makeProductMeta(scene *Scene) *ProductMeta { meta := &ProductMeta{ Satellite: "SJY01", @@ -89,7 +110,7 @@ func (r *Registrator) makeProductMeta(scene *Scene) *ProductMeta { return meta } -func (r *Registrator) writeProductMeta(productMeta *ProductMeta, filename string) error { +func writeProductMeta(productMeta *ProductMeta, filename string) error { output, err := xml.MarshalIndent(productMeta, "", " ") if err != nil { return err diff --git a/pkg/producer/scenes.go b/pkg/producer/scenes.go index 451c785..cf362e6 100644 --- a/pkg/producer/scenes.go +++ b/pkg/producer/scenes.go @@ -159,7 +159,7 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e return panScenes, mssScenes, err } -func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) error { +func (r *Registrator) OutputL1A(panScenes []*Scene, mssScenes []*Scene) error { var fc geojson.FeatureCollection for i, scene := range panScenes { dir := filepath.Join(r.Params.OutputDir, scene.SceneIndex, "PAN") @@ -182,7 +182,7 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e scene.Meta.DataSize = sizeOfFile(scene.Tiff) metaFile := strings.Replace(scene.Tiff, ".tiff", ".meta.xml", 1) - r.writeProductMeta(scene.Meta, metaFile) + writeProductMeta(scene.Meta, metaFile) jpg := strings.Replace(scene.Tiff, ".tiff", ".jpg", 1) GTiffToJPG(scene.Tiff, jpg, "PAN", false) r.report.Scenes = append(r.report.Scenes, ReportScene{ @@ -233,7 +233,7 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e scene.Meta.DataSize = sizeOfFile(scene.Tiff) metaFile := strings.Replace(scene.Tiff, ".tiff", ".meta.xml", 1) - r.writeProductMeta(scene.Meta, metaFile) + writeProductMeta(scene.Meta, metaFile) jpg := strings.Replace(scene.Tiff, ".tiff", ".jpg", 1) GTiffToJPG(scene.Tiff, jpg, "MSS", false) r.report.Scenes = append(r.report.Scenes, ReportScene{ diff --git a/pkg/producer/warp.go b/pkg/producer/warp.go index b607b5d..3ac46a5 100644 --- a/pkg/producer/warp.go +++ b/pkg/producer/warp.go @@ -1,13 +1,13 @@ package producer import ( - "math" "os" "os/exec" "path/filepath" "strconv" "strings" + "github.com/airbusgeo/godal" log "github.com/sirupsen/logrus" "starwiz.cn/sjy01/image-proc/pkg/dem" ) @@ -15,7 +15,7 @@ import ( // gdalwarp -rpc -to "RPC_HEIGHT=2375" in.tif out.tif // gdalwarp -rpc -to "RPC_DEM=/path/to/dem/gdlebm.tif" in.tif out.tif -func Warp(in, out, meta, rpb, demtif string) error { +func L1AtoL2(in, out, meta, rpb, demtif string) error { meta, _ = filepath.Abs(meta) m, err := ParseProductMeta(meta) if err != nil { @@ -27,7 +27,7 @@ func Warp(in, out, meta, rpb, demtif string) error { log.Infof("Warping %s to %s", in, out) minEv, maxEv := dem.Dem1KmLT.MinMaxElevationInRect( - computeExtend(&m.Corners), + m.Corners.Extend(), ) log.Infof("DEM min/max elevation: %f/%f", minEv, maxEv) var cmd *exec.Cmd @@ -49,27 +49,106 @@ func Warp(in, out, meta, rpb, demtif string) error { if strings.Contains(m.ProductID, "PAN") { sensor = "PAN" } + + corners, err := computeBound(out) + if err != nil { + return err + } + + id := filepath.Base(out) + id = strings.Replace(id, ".tiff", "", 1) + m.ProductID = id + m.ProductLevel = "L2" + m.Corners = *corners + xmlfile := filepath.Join(dir, id+".meta.xml") + writeProductMeta(m, xmlfile) + GTiffToJPG(out, strings.Replace(out, ".tiff", ".jpg", 1), sensor, false) return nil } -func computeExtend(corners *Corners) (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) +func computeBound(tif string) (*Corners, error) { + ds, err := godal.Open(tif) + if err != nil { + log.Printf("Error opening TIFF file %s: %v", tif, err) + return nil, err + } + defer ds.Close() - lat1 = math.Max(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude) - lat1 = math.Max(lat1, corners.UpperRight.Latitude) - lat1 = math.Max(lat1, corners.LowerRight.Latitude) + geotransform, err := ds.GeoTransform() + if err != nil { + return nil, err + } - lng2 = math.Max(corners.UpperLeft.Longitude, corners.LowerLeft.Longitude) - lng2 = math.Max(lng2, corners.UpperRight.Longitude) - lng2 = math.Max(lng2, corners.LowerRight.Longitude) + // Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2]; + // Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5]; - lat2 = math.Min(corners.UpperLeft.Latitude, corners.LowerLeft.Latitude) - lat2 = math.Min(lat2, corners.UpperRight.Latitude) - lat2 = math.Min(lat2, corners.LowerRight.Latitude) + structure := ds.Structure() + xmin := geotransform[0] + ymax := geotransform[3] + xmax := xmin + geotransform[1]*float64(structure.SizeX) + ymin := ymax + geotransform[5]*float64(structure.SizeY) + log.Infof("extend: %f,%f,%f,%f", xmin, ymin, xmax, ymax) + log.Infof("geotransform: %v", geotransform) - return + band := ds.Bands()[0] + width := structure.SizeX + height := structure.SizeY + data := make([]uint16, width*height) + err = band.Read(0, 0, data, width, height) + if err != nil { + return nil, err + } + + nodata, _ := band.NoData() + + index1 := 0 + for x := 0; x < width; x++ { + value := data[0*width+x] + if value != uint16(nodata) { + index1 = x + break + } + } + + index2 := 0 + for x := 0; x < width; x++ { + value := data[(height-1)*width+x] + if value != uint16(nodata) { + index2 = x + break + } + } + + index3 := 0 + for y := 0; y < height; y++ { + value := data[y*width] + if value != uint16(nodata) { + index3 = y + break + } + } + + index4 := 0 + for y := 0; y < height; y++ { + value := data[y*width+width-1] + if value != uint16(nodata) { + index4 = y + break + } + } + + var corners Corners + + corners.UpperLeft.Longitude = xmin + float64(index1)*geotransform[1] + corners.UpperLeft.Latitude = ymax + corners.UpperRight.Longitude = xmax + corners.UpperRight.Latitude = ymax + float64(index4)*geotransform[5] + corners.LowerLeft.Longitude = xmin + corners.LowerLeft.Latitude = ymax + float64(index3)*geotransform[5] + corners.LowerRight.Longitude = xmin + float64(index2)*geotransform[1] + corners.LowerRight.Latitude = ymin + + return &corners, nil }