L1A -> L2

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
 }