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 { 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) } 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" } 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, false) 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 }