106 lines
2.1 KiB
Go
106 lines
2.1 KiB
Go
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 := 0; 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
|
|
hDataset.Close()
|
|
|
|
log.Infof("DEM loaded: %s, %d x %d, wUnit: %.2f, hUnit: %.2f", d.File, d.width, d.height, d.wUnit, d.hUnit)
|
|
|
|
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]
|
|
|
|
}
|
|
|
|
func (d *Dem1Km) MinMaxElevationInRect(lng1, lat1, lng2, lat2 float64) (float64, float64) {
|
|
minElev := math.MaxFloat32
|
|
maxElev := -math.MaxFloat32
|
|
|
|
lng1 = math.Abs(lng1 + 180.0)
|
|
lat1 = math.Abs(lat1 - 90.0)
|
|
w1 := int(lng1 / float64(d.wUnit))
|
|
h1 := int(lat1 / float64(d.hUnit))
|
|
|
|
lng2 = math.Abs(lng2 + 180.0)
|
|
lat2 = math.Abs(lat2 - 90.0)
|
|
w2 := int(lng2 / float64(d.wUnit))
|
|
h2 := int(lat2 / float64(d.hUnit))
|
|
|
|
if w1 > d.width-1 || h1 > d.height-1 || w2 > d.width-1 || h2 > d.height-1 {
|
|
return 0.0, 0.0
|
|
}
|
|
|
|
for w := w1; w <= w2; w++ {
|
|
for h := h1; h <= h2; h++ {
|
|
if d.LT[h][w] < float32(minElev) {
|
|
minElev = float64(d.LT[h][w])
|
|
}
|
|
if d.LT[h][w] > float32(maxElev) {
|
|
maxElev = float64(d.LT[h][w])
|
|
}
|
|
}
|
|
}
|
|
|
|
return minElev, maxElev
|
|
}
|