使用 gdalwarp 生产 L2

cmd/warp.go

@@ -0,0 +1,46 @@
+package main
+
+import (
+	"os"
+	"path/filepath"
+
+	"github.com/airbusgeo/godal"
+	"github.com/sirupsen/logrus"
+	"github.com/spf13/cobra"
+	"starwiz.cn/sjy01/image-proc/pkg/config"
+	"starwiz.cn/sjy01/image-proc/pkg/dem"
+	"starwiz.cn/sjy01/image-proc/pkg/producer"
+)
+
+var (
+	warpIn   string
+	warpOut  string
+	warpMeta string
+	warpRPC  string
+	warpDEM  string
+)
+
+var cmdWarp = &cobra.Command{
+	Use:   "warp",
+	Short: "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)
+		if err != nil {
+			logrus.Error("Generate L2 product failed:", err)
+			os.Exit(10)
+		}
+	},
+}
+
+func init() {
+	rootCmd.AddCommand(cmdWarp)
+
+	cmdWarp.Flags().StringVarP(&warpIn, "input", "i", "", "input image file")
+	cmdWarp.Flags().StringVarP(&warpOut, "output", "o", "", "output image file")
+	cmdWarp.Flags().StringVarP(&warpMeta, "meta", "m", "", "rpc meta file")
+	cmdWarp.Flags().StringVarP(&warpRPC, "rpb", "r", "", "rpb file")
+	cmdWarp.Flags().StringVarP(&warpDEM, "dem", "d", "", "dem file")
+}

pkg/producer/meta.go

@@ -102,3 +102,19 @@ func (r *Registrator) writeProductMeta(productMeta *ProductMeta, filename string
 
 	return nil
 }
+
+func ParseProductMeta(filename string) (*ProductMeta, error) {
+	f, err := os.Open(filename)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+
+	var meta ProductMeta
+	err = xml.NewDecoder(f).Decode(&meta)
+	if err != nil {
+		return nil, err
+	}
+
+	return &meta, nil
+}

pkg/producer/warp.go

@@ -0,0 +1,75 @@
+package producer
+
+import (
+	"math"
+	"os"
+	"os/exec"
+	"path/filepath"
+	"strconv"
+	"strings"
+
+	log "github.com/sirupsen/logrus"
+	"starwiz.cn/sjy01/image-proc/pkg/dem"
+)
+
+// 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 {
+	meta, _ = filepath.Abs(meta)
+	m, err := ParseProductMeta(meta)
+	if err != nil {
+		return err
+	}
+
+	in, _ = filepath.Abs(in)
+	out, _ = filepath.Abs(out)
+
+	log.Infof("Warping %s to %s", in, out)
+	minEv, maxEv := dem.Dem1KmLT.MinMaxElevationInRect(
+		computeExtend(&m.Corners),
+	)
+	log.Infof("DEM min/max elevation: %f/%f", minEv, maxEv)
+	var cmd *exec.Cmd
+	if demtif == "" {
+		cmd = exec.Command("gdalwarp", "-rpc", "-to", "RPC_HEIGHT="+strconv.Itoa(int(maxEv)), in, out)
+	} else {
+		cmd = exec.Command("gdalwarp", "-rpc", "-to", "RPC_DEM="+demtif, in, out)
+	}
+
+	dir := filepath.Dir(out)
+	os.MkdirAll(dir, 0755)
+
+	err = cmd.Run()
+	if err != nil {
+		return err
+	}
+
+	sensor := "MSS"
+	if strings.Contains(m.ProductID, "PAN") {
+		sensor = "PAN"
+	}
+	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)
+
+	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
+}