1// Package terrain samples ground elevation along a path from the AWS Terrain2// Tiles dataset (Terrarium encoding; SRTM, USGS NED and other sources).3// Tiles are cached on disk, so each area is downloaded once.4package terrain56import (7 "bytes"8 "context"9 "fmt"10 "image"11 "image/png"12 "math"13 "net/http"14 "os"15 "path/filepath"16 "strconv"17 "sync"1819 "airwaves/internal/geo"20 "airwaves/internal/web"21)2223const (24 tileURL = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/%d/%d/%d.png"25 // zoom 10 gives roughly 120 m pixels at US latitudes.26 zoom = 1027 tileSize = 25628 spacingKm = 0.1529)3031// Profile is ground elevation sampled at equal steps from the receiver32// (index 0) to the transmitter (last index).33type Profile struct {34 DistanceKm float64 `json:"distanceKm"`35 Elevations []float64 `json:"elevations"`36}3738type tileKey struct{ x, y int }3940type tile struct {41 once sync.Once42 elev []float3243 err error44}4546// Source reads elevations, downloading tiles on demand.47type Source struct {48 HTTP *http.Client49 Dir string // tile cache directory5051 mu sync.Mutex52 tiles map[tileKey]*tile53}5455// NewSource caches tiles under dir.56func NewSource(c *http.Client, dir string) *Source {57 return &Source{HTTP: c, Dir: dir, tiles: make(map[tileKey]*tile)}58}5960// Profile samples the ground from rx to tx.61func (s *Source) Profile(ctx context.Context, rx, tx geo.Point) (Profile, error) {62 d := geo.DistanceKm(rx, tx)63 n := max(int(math.Ceil(d/spacingKm))+1, 16)64 elev := make([]float64, n)65 for i := range n {66 e, err := s.Elevation(ctx, geo.Intermediate(rx, tx, float64(i)/float64(n-1)))67 if err != nil {68 return Profile{}, err69 }70 elev[i] = e71 }72 return Profile{DistanceKm: d, Elevations: elev}, nil73}7475// Elevation returns meters above sea level at p.76func (s *Source) Elevation(ctx context.Context, p geo.Point) (float64, error) {77 scale := float64(int(1)<<zoom) * tileSize78 lat := p.Lat * math.Pi / 18079 px := (p.Lon + 180) / 360 * scale80 py := (1 - math.Log(math.Tan(lat)+1/math.Cos(lat))/math.Pi) / 2 * scale81 ix, iy := int(px), int(py)82 t, err := s.tile(ctx, tileKey{ix / tileSize, iy / tileSize})83 if err != nil {84 return 0, err85 }86 return float64(t[(iy%tileSize)*tileSize+ix%tileSize]), nil87}8889func (s *Source) tile(ctx context.Context, k tileKey) ([]float32, error) {90 s.mu.Lock()91 t, ok := s.tiles[k]92 if !ok {93 t = &tile{}94 s.tiles[k] = t95 }96 s.mu.Unlock()97 t.once.Do(func() { t.elev, t.err = s.load(ctx, k) })98 if t.err != nil {99 // Let a later call retry a failed download.100 s.mu.Lock()101 if s.tiles[k] == t {102 delete(s.tiles, k)103 }104 s.mu.Unlock()105 }106 return t.elev, t.err107}108109func (s *Source) load(ctx context.Context, k tileKey) ([]float32, error) {110 path := filepath.Join(s.Dir, strconv.Itoa(zoom), strconv.Itoa(k.x), strconv.Itoa(k.y)+".png")111 raw, err := os.ReadFile(path)112 if err != nil {113 raw, err = web.Get(ctx, s.HTTP, fmt.Sprintf(tileURL, zoom, k.x, k.y))114 if err != nil {115 return nil, fmt.Errorf("terrain tile %d/%d/%d: %w", zoom, k.x, k.y, err)116 }117 if err := os.MkdirAll(filepath.Dir(path), 0o755); err == nil {118 _ = os.WriteFile(path, raw, 0o644)119 }120 }121 img, err := png.Decode(bytes.NewReader(raw))122 if err != nil {123 return nil, fmt.Errorf("decode terrain tile %d/%d/%d: %w", zoom, k.x, k.y, err)124 }125 return decodeTerrarium(img), nil126}127128// decodeTerrarium converts pixels to meters: (R*256 + G + B/256) - 32768.129func decodeTerrarium(img image.Image) []float32 {130 b := img.Bounds()131 out := make([]float32, tileSize*tileSize)132 for y := range tileSize {133 for x := range tileSize {134 r, g, bl, _ := img.At(b.Min.X+x, b.Min.Y+y).RGBA()135 out[y*tileSize+x] = float32(float64(r>>8)*256+float64(g>>8)+float64(bl>>8)/256) - 32768136 }137 }138 return out139}