add dem

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",
+		},
 	}
 }

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]
+
+}

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
+}