sjy01-image-proc/pkg/producer/warp.go

package producer

import (
	"os"
	"os/exec"
	"path/filepath"
	"strconv"
	"strings"

	"github.com/airbusgeo/godal"
	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 L1AtoL2(in, out, meta, rpb, demtif string) error {
	var m *ProductMeta
	var err error

	sensor := "MSS"
	reverse := true
	if meta == "" {
		sensor = "FUS"
		reverse = true
		m = &ProductMeta{}
	} else {
		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(
		m.Corners.Extend(),
	)
	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)
	}

	// gdalwarp  -overwrite -t_srs EPSG:4326 -r average -of GTiff -rpc
	cmd = exec.Command("gdalwarp", "-overwrite", "-t_srs", "EPSG:4326", "-r", "average", "-of", "GTiff", "-rpc", in, out)
	log.Info("Running command: ", cmd.String())

	dir := filepath.Dir(out)
	os.MkdirAll(dir, 0755)

	err = cmd.Run()
	if err != nil {
		return err
	}

	if strings.Contains(m.ProductID, "PAN") {
		sensor = "PAN"
		reverse = false
	} else if strings.Contains(m.ProductID, "FUS") {
		sensor = "FUS"
		reverse = false
	}

	corners, width, height, 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
	m.Width = width
	m.Height = height
	xmlfile := filepath.Join(dir, id+".meta.xml")
	writeProductMeta(m, xmlfile)

	GTiffToJPG(out, strings.Replace(out, ".tiff", ".jpg", 1), sensor, reverse)

	return nil
}

func computeBound(tif string) (*Corners, int, int, error) {
	ds, err := godal.Open(tif)
	if err != nil {
		log.Printf("Error opening TIFF file %s: %v", tif, err)
		return nil, -1, -1, err
	}
	defer ds.Close()

	geotransform, err := ds.GeoTransform()
	if err != nil {
		return nil, -1, -1, err
	}

	// Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
	// Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];

	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)

	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, -1, -1, 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, width, height, nil
}