From ea23839bf2dada9022e4544452cee9b34a33604a Mon Sep 17 00:00:00 2001
From: nuknal <linktolane@gmail.com>
Date: Mon, 19 Aug 2024 18:05:55 +0800
Subject: [PATCH] add dem

---
 .gitignore           |  2 ++
 cmd/proc.go          |  4 +++
 config/config.yaml   |  5 +++-
 pkg/config/config.go | 12 ++++++++
 pkg/dem/dem.go       | 70 ++++++++++++++++++++++++++++++++++++++++++++
 pkg/utils/tiff.go    | 24 +++++++++++++++
 6 files changed, 116 insertions(+), 1 deletion(-)
 create mode 100644 pkg/dem/dem.go

diff --git a/.gitignore b/.gitignore
index bc7e2b0..ebc359a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,5 @@ bin
 log
 deployment
 build/go
+*.tif
+*.tiff
diff --git a/cmd/proc.go b/cmd/proc.go
index 6c119ee..a0c97af 100644
--- a/cmd/proc.go
+++ b/cmd/proc.go
@@ -11,6 +11,7 @@ import (
 	"github.com/spf13/cobra"
 	"starwiz.cn/sjy01/image-proc/pkg/calculator"
 	"starwiz.cn/sjy01/image-proc/pkg/config"
+	"starwiz.cn/sjy01/image-proc/pkg/dem"
 	producer "starwiz.cn/sjy01/image-proc/pkg/producer"
 )
 
@@ -35,6 +36,9 @@ var procCmd = &cobra.Command{
 
 		logrus.SetLevel(config.GCONFIG.Log.LogLevel)
 
+		// dem
+		dem.Dem1KmLT = dem.NewDem1Km(config.GCONFIG.Dem.Dem1Km)
+
 		calculator.EOP = calculator.NewEOPTable()
 		calculator.EOP.Load(eopData, eopp5Line)
 
diff --git a/config/config.yaml b/config/config.yaml
index 4f1f6a2..d0089a3 100644
--- a/config/config.yaml
+++ b/config/config.yaml
@@ -17,4 +17,7 @@ radiation:
   min_hist_level: 0.3
   max_hist_level: 0.6
   hf_band_stop_width: 24 # 带阻滤波器宽度
-  scene_moment_matching: false
\ No newline at end of file
+  scene_moment_matching: false
+
+dem:
+  dem_1km: "dem/gdlebm.tif"
\ No newline at end of file
diff --git a/pkg/config/config.go b/pkg/config/config.go
index aa4a210..3b417d9 100644
--- a/pkg/config/config.go
+++ b/pkg/config/config.go
@@ -10,6 +10,7 @@ type Config struct {
 	Radiation      RadiationConfig      `yaml:"radiation" mapstructure:"radiation"`
 	BrowserImg     BrowserImgConfig     `yaml:"browser_img" mapstructure:"browser_img"`
 	Log            LogConfig            `yaml:"log" mapstructure:"log"`
+	Dem            DemConfig            `yaml:"dem" mapstructure:"dem"`
 }
 
 type LogConfig struct {
@@ -39,6 +40,12 @@ type RadiationConfig struct {
 	HFBandStopWidth     int     `yaml:"hf_band_stop_width" mapstructure:"hf_band_stop_width"`
 }
 
+type DemConfig struct {
+	DemType string `yaml:"dem_type" mapstructure:"dem_type"`
+	DemDir  string `yaml:"dem_dir" mapstructure:"dem_dir"`
+	Dem1Km  string `yaml:"dem_1km" mapstructure:"dem_1km"`
+}
+
 type BrowserImgConfig struct {
 	Enhancement        map[string]string `yaml:"enhancement" mapstructure:"enhancement"`
 	CumulativeCutLower float64           `yaml:"cumulative_cut_lower" mapstructure:"cumulative_cut_lower"`
@@ -82,5 +89,10 @@ func init() {
 			CumulativeCutLower: 0.01,
 			CumulativeCutUpper: 0.99,
 		},
+		Dem: DemConfig{
+			DemType: "SRTM",
+			DemDir:  "dem/SRTM",
+			Dem1Km:  "dem/gdlebm.tif",
+		},
 	}
 }
diff --git a/pkg/dem/dem.go b/pkg/dem/dem.go
new file mode 100644
index 0000000..7384c80
--- /dev/null
+++ b/pkg/dem/dem.go
@@ -0,0 +1,70 @@
+package dem
+
+import (
+	"math"
+
+	"github.com/airbusgeo/godal"
+	log "github.com/sirupsen/logrus"
+)
+
+type Dem1Km struct {
+	File          string
+	LT            [][]float32
+	width, height int
+	wUnit         float32 // degree per pixel
+	hUnit         float32 // degree per pixel
+}
+
+var Dem1KmLT *Dem1Km
+
+func NewDem1Km(file string) *Dem1Km {
+	d := &Dem1Km{
+		File: file,
+	}
+
+	if err := d.Load(); err != nil {
+		log.Error("cannot load DEM file", err)
+	}
+	return d
+}
+
+func (d *Dem1Km) Load() error {
+	hDataset, err := godal.Open(d.File)
+	if err != nil {
+		return err
+	}
+
+	band := hDataset.Bands()[0]
+
+	structure := band.Structure()
+
+	for y := 1; y <= structure.SizeY; y++ {
+		scanline := make([]float32, structure.SizeX)
+		err = band.Read(0, y, scanline, structure.SizeX, 1)
+		if err != nil {
+			return err
+		}
+		d.LT = append(d.LT, scanline)
+	}
+
+	d.wUnit = 360.0 / float32(structure.SizeX)
+	d.hUnit = 180.0 / float32(structure.SizeY)
+	d.width = structure.SizeX
+	d.height = structure.SizeY
+
+	return nil
+}
+
+func (d *Dem1Km) Elevation(lng, lat float64) float32 {
+	lng = math.Abs(lng + 180.0)
+	lat = math.Abs(lat - 90.0)
+	w := int(lng / float64(d.wUnit))
+	h := int(lat / float64(d.hUnit))
+
+	if w > d.width-1 || h > d.height-1 {
+		return 0.0
+	}
+
+	return d.LT[h][w]
+
+}
diff --git a/pkg/utils/tiff.go b/pkg/utils/tiff.go
index 4737199..f62ef87 100644
--- a/pkg/utils/tiff.go
+++ b/pkg/utils/tiff.go
@@ -150,3 +150,27 @@ func setGeoTransform(ds *godal.Dataset, topLeftLng, topLeftLat, resolution float
 	ds.SetMetadata("TIFFTAG_DATETIME", time.Now().String())
 	ds.SetMetadata("TIFFTAG_SOFTWARE", "StarWiz-SJY01-IMAGE-PROC")
 }
+
+func ReadTiff(tifFile string) ([][][]float32, error) {
+	hDataset, err := godal.Open(tifFile)
+	if err != nil {
+		return nil, err
+	}
+
+	var data [][][]float32
+	for _, band := range hDataset.Bands() {
+		structure := band.Structure()
+		var dataBand [][]float32
+		for y := 1; y <= structure.SizeY; y++ {
+			scanline := make([]float32, structure.SizeX)
+			err = band.Read(0, y, scanline, structure.SizeX, 1)
+			if err != nil {
+				return nil, err
+			}
+			dataBand = append(dataBand, scanline)
+		}
+
+		data = append(data, dataBand)
+	}
+	return data, nil
+}