391 lines
10 KiB
Go
391 lines
10 KiB
Go
package server
|
|
|
|
import (
|
|
"bytes"
|
|
"encoding/json"
|
|
"fmt"
|
|
"image"
|
|
"image/color"
|
|
"image/png"
|
|
"math"
|
|
"net/http"
|
|
"strconv"
|
|
"strings"
|
|
|
|
"github.com/blakeridgway/heloha/internal/radar"
|
|
"github.com/go-chi/chi/v5"
|
|
)
|
|
|
|
const tileSize = 256
|
|
|
|
// TileHandler serves /api/v1/tile/{site}/{z}/{x}/{y}.png (single-site, reflectivity).
|
|
func TileHandler(store *RadarStore, cache *TileCache) http.HandlerFunc {
|
|
return func(w http.ResponseWriter, r *http.Request) {
|
|
site := strings.ToLower(chi.URLParam(r, "site"))
|
|
z, err1 := strconv.Atoi(chi.URLParam(r, "z"))
|
|
x, err2 := strconv.Atoi(chi.URLParam(r, "x"))
|
|
y, err3 := strconv.Atoi(chi.URLParam(r, "y"))
|
|
if err1 != nil || err2 != nil || err3 != nil {
|
|
http.Error(w, "invalid tile coords", http.StatusBadRequest)
|
|
return
|
|
}
|
|
|
|
key := site + "/refl/latest/" + strconv.Itoa(z) + "/" + strconv.Itoa(x) + "/" + strconv.Itoa(y)
|
|
if data, ok := cache.Get(key); ok {
|
|
w.Header().Set("Content-Type", "image/png")
|
|
w.Write(data)
|
|
return
|
|
}
|
|
|
|
product := store.GetLatest(site, "reflectivity")
|
|
if product == nil {
|
|
http.Error(w, "no radar data", http.StatusServiceUnavailable)
|
|
return
|
|
}
|
|
|
|
img := renderTile(product, z, x, y, dbzColor)
|
|
writeTile(w, cache, key, img)
|
|
}
|
|
}
|
|
|
|
// CompositeTileHandler serves:
|
|
//
|
|
// /api/v1/tile/composite/{product}/{frame}/{z}/{x}/{y}.png
|
|
// /api/v1/tile/composite/{z}/{x}/{y}.png (legacy: product=reflectivity, frame=latest)
|
|
func CompositeTileHandler(store *RadarStore, cache *TileCache) http.HandlerFunc {
|
|
return func(w http.ResponseWriter, r *http.Request) {
|
|
product := chi.URLParam(r, "product")
|
|
if product == "" {
|
|
product = "reflectivity"
|
|
}
|
|
frameStr := chi.URLParam(r, "frame")
|
|
|
|
z, err1 := strconv.Atoi(chi.URLParam(r, "z"))
|
|
x, err2 := strconv.Atoi(chi.URLParam(r, "x"))
|
|
y, err3 := strconv.Atoi(chi.URLParam(r, "y"))
|
|
if err1 != nil || err2 != nil || err3 != nil {
|
|
http.Error(w, "invalid tile coords", http.StatusBadRequest)
|
|
return
|
|
}
|
|
|
|
key := fmt.Sprintf("composite/%s/%s/%d/%d/%d", product, frameStr, z, x, y)
|
|
if data, ok := cache.Get(key); ok {
|
|
w.Header().Set("Content-Type", "image/png")
|
|
w.Write(data)
|
|
return
|
|
}
|
|
|
|
var products []*radar.RadarProduct
|
|
if frameStr == "" || frameStr == "latest" {
|
|
products = store.GetAllLatest(product)
|
|
} else {
|
|
frameIdx, err := strconv.Atoi(frameStr)
|
|
if err != nil {
|
|
http.Error(w, "invalid frame", http.StatusBadRequest)
|
|
return
|
|
}
|
|
products = store.GetAllAtFrame(product, frameIdx)
|
|
}
|
|
|
|
if len(products) == 0 {
|
|
http.Error(w, "no radar data", http.StatusServiceUnavailable)
|
|
return
|
|
}
|
|
|
|
colorFn := dbzColor
|
|
if product == "velocity" {
|
|
colorFn = velColor
|
|
}
|
|
|
|
img := renderCompositeTile(products, z, x, y, colorFn)
|
|
writeTile(w, cache, key, img)
|
|
}
|
|
}
|
|
|
|
// FramesHandler serves /api/v1/frames?product=reflectivity|velocity
|
|
func FramesHandler(store *RadarStore) http.HandlerFunc {
|
|
return func(w http.ResponseWriter, r *http.Request) {
|
|
product := r.URL.Query().Get("product")
|
|
if product == "" {
|
|
product = "reflectivity"
|
|
}
|
|
count, frames := store.FrameMeta(product)
|
|
w.Header().Set("Content-Type", "application/json")
|
|
json.NewEncoder(w).Encode(map[string]any{
|
|
"product": product,
|
|
"frame_count": count,
|
|
"frames": frames,
|
|
})
|
|
}
|
|
}
|
|
|
|
// speckleFilter removes pixels with fewer than minNeighbors non-transparent
|
|
// 8-connected neighbors. Running multiple passes progressively erodes small
|
|
// noise clusters while large coherent precipitation areas survive intact.
|
|
func speckleFilter(src *image.RGBA) *image.RGBA {
|
|
const (
|
|
passes = 4
|
|
minNeighbors = 4
|
|
)
|
|
img := src
|
|
for pass := 0; pass < passes; pass++ {
|
|
dst := image.NewRGBA(img.Bounds())
|
|
for y := 0; y < tileSize; y++ {
|
|
for x := 0; x < tileSize; x++ {
|
|
c := img.RGBAAt(x, y)
|
|
if c.A == 0 {
|
|
continue
|
|
}
|
|
n := 0
|
|
for dy := -1; dy <= 1; dy++ {
|
|
for dx := -1; dx <= 1; dx++ {
|
|
if dx == 0 && dy == 0 {
|
|
continue
|
|
}
|
|
nx, ny := x+dx, y+dy
|
|
if nx >= 0 && nx < tileSize && ny >= 0 && ny < tileSize &&
|
|
img.RGBAAt(nx, ny).A > 0 {
|
|
n++
|
|
}
|
|
}
|
|
}
|
|
if n >= minNeighbors {
|
|
dst.SetRGBA(x, y, c)
|
|
}
|
|
}
|
|
}
|
|
img = dst
|
|
}
|
|
return img
|
|
}
|
|
|
|
func writeTile(w http.ResponseWriter, cache *TileCache, key string, img *image.RGBA) {
|
|
var buf bytes.Buffer
|
|
if err := png.Encode(&buf, img); err != nil {
|
|
http.Error(w, "encode error", http.StatusInternalServerError)
|
|
return
|
|
}
|
|
data := buf.Bytes()
|
|
cache.Set(key, data)
|
|
w.Header().Set("Content-Type", "image/png")
|
|
w.Write(data)
|
|
}
|
|
|
|
// renderCompositeTile builds a mosaic tile: each pixel gets a value from the
|
|
// nearest radar site within 230 km, colored by colorFn.
|
|
func renderCompositeTile(products []*radar.RadarProduct, z, x, y int, colorFn func(float32) color.RGBA) *image.RGBA {
|
|
img := image.NewRGBA(image.Rect(0, 0, tileSize, tileSize))
|
|
|
|
north, west := tileToLatLon(z, x, y)
|
|
south, east := tileToLatLon(z, x+1, y+1)
|
|
latSpan := north - south
|
|
lonSpan := east - west
|
|
|
|
const maxRangeKm = 230.0
|
|
|
|
for py := 0; py < tileSize; py++ {
|
|
lat := north - float64(py)/tileSize*latSpan
|
|
for px := 0; px < tileSize; px++ {
|
|
lon := west + float64(px)/tileSize*lonSpan
|
|
|
|
var bestVal float32 = float32(math.NaN())
|
|
bestDist := math.MaxFloat64
|
|
|
|
for _, p := range products {
|
|
bearing, distKm := bearingDist(p.SiteLat, p.SiteLon, lat, lon)
|
|
if distKm >= maxRangeKm || distKm >= bestDist {
|
|
continue
|
|
}
|
|
val := interpolateValue(p.Radials, bearing, distKm)
|
|
bestDist = distKm
|
|
bestVal = val
|
|
}
|
|
|
|
if math.IsNaN(float64(bestVal)) {
|
|
continue
|
|
}
|
|
img.SetRGBA(px, py, colorFn(bestVal))
|
|
}
|
|
}
|
|
return speckleFilter(img)
|
|
}
|
|
|
|
// renderTile renders a single-site tile.
|
|
func renderTile(p *radar.RadarProduct, z, x, y int, colorFn func(float32) color.RGBA) *image.RGBA {
|
|
img := image.NewRGBA(image.Rect(0, 0, tileSize, tileSize))
|
|
|
|
north, west := tileToLatLon(z, x, y)
|
|
south, east := tileToLatLon(z, x+1, y+1)
|
|
latSpan := north - south
|
|
lonSpan := east - west
|
|
|
|
for py := 0; py < tileSize; py++ {
|
|
lat := north - float64(py)/tileSize*latSpan
|
|
for px := 0; px < tileSize; px++ {
|
|
lon := west + float64(px)/tileSize*lonSpan
|
|
bearing, distKm := bearingDist(p.SiteLat, p.SiteLon, lat, lon)
|
|
val := interpolateValue(p.Radials, bearing, distKm)
|
|
if math.IsNaN(float64(val)) {
|
|
continue
|
|
}
|
|
img.SetRGBA(px, py, colorFn(val))
|
|
}
|
|
}
|
|
return speckleFilter(img)
|
|
}
|
|
|
|
// interpolateValue bilinearly interpolates between the two nearest radials and adjacent bins.
|
|
func interpolateValue(radials []radar.Radial, bearing, distKm float64) float32 {
|
|
if len(radials) == 0 {
|
|
return float32(math.NaN())
|
|
}
|
|
|
|
best := 0
|
|
bestDiff := azDiff(radials[0].Azimuth, bearing)
|
|
for i := 1; i < len(radials); i++ {
|
|
if d := azDiff(radials[i].Azimuth, bearing); d < bestDiff {
|
|
bestDiff = d
|
|
best = i
|
|
}
|
|
}
|
|
|
|
diff := bearing - radials[best].Azimuth
|
|
if diff > 180 {
|
|
diff -= 360
|
|
}
|
|
if diff < -180 {
|
|
diff += 360
|
|
}
|
|
|
|
var adj int
|
|
if diff >= 0 {
|
|
adj = (best + 1) % len(radials)
|
|
} else {
|
|
adj = (best - 1 + len(radials)) % len(radials)
|
|
diff = -diff
|
|
}
|
|
|
|
span := azDiff(radials[best].Azimuth, radials[adj].Azimuth)
|
|
t := 0.0
|
|
if span > 0 {
|
|
t = math.Min(diff/span, 1.0)
|
|
}
|
|
|
|
v1 := radialValueAt(&radials[best], distKm)
|
|
v2 := radialValueAt(&radials[adj], distKm)
|
|
|
|
switch {
|
|
case math.IsNaN(float64(v1)) && math.IsNaN(float64(v2)):
|
|
return float32(math.NaN())
|
|
case math.IsNaN(float64(v1)):
|
|
return v2
|
|
case math.IsNaN(float64(v2)):
|
|
return v1
|
|
}
|
|
return float32(float64(v1)*(1-t) + float64(v2)*t)
|
|
}
|
|
|
|
func radialValueAt(r *radar.Radial, distKm float64) float32 {
|
|
if distKm < r.RangeKm || r.BinSizeKm == 0 {
|
|
return float32(math.NaN())
|
|
}
|
|
exactBin := (distKm - r.RangeKm) / r.BinSizeKm
|
|
b0 := int(exactBin)
|
|
if b0 >= len(r.Values) {
|
|
return float32(math.NaN())
|
|
}
|
|
v0 := r.Values[b0]
|
|
b1 := b0 + 1
|
|
if b1 >= len(r.Values) || math.IsNaN(float64(v0)) {
|
|
return v0
|
|
}
|
|
v1 := r.Values[b1]
|
|
if math.IsNaN(float64(v1)) {
|
|
return v0
|
|
}
|
|
frac := float32(exactBin - float64(b0))
|
|
return v0*(1-frac) + v1*frac
|
|
}
|
|
|
|
func bearingDist(lat1, lon1, lat2, lon2 float64) (bearing, distKm float64) {
|
|
const earthR = 6371.0
|
|
lat1r := lat1 * math.Pi / 180
|
|
lat2r := lat2 * math.Pi / 180
|
|
dlat := (lat2 - lat1) * math.Pi / 180
|
|
dlon := (lon2 - lon1) * math.Pi / 180
|
|
a := math.Sin(dlat/2)*math.Sin(dlat/2) +
|
|
math.Cos(lat1r)*math.Cos(lat2r)*math.Sin(dlon/2)*math.Sin(dlon/2)
|
|
distKm = earthR * 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
|
|
y := math.Sin(dlon) * math.Cos(lat2r)
|
|
x := math.Cos(lat1r)*math.Sin(lat2r) - math.Sin(lat1r)*math.Cos(lat2r)*math.Cos(dlon)
|
|
bearing = math.Atan2(y, x) * 180 / math.Pi
|
|
if bearing < 0 {
|
|
bearing += 360
|
|
}
|
|
return
|
|
}
|
|
|
|
func azDiff(a, b float64) float64 {
|
|
d := math.Abs(a - b)
|
|
if d > 180 {
|
|
d = 360 - d
|
|
}
|
|
return d
|
|
}
|
|
|
|
func tileToLatLon(z, x, y int) (lat, lon float64) {
|
|
n := math.Pow(2, float64(z))
|
|
lon = float64(x)/n*360.0 - 180.0
|
|
latRad := math.Atan(math.Sinh(math.Pi * (1 - 2*float64(y)/n)))
|
|
lat = latRad * 180.0 / math.Pi
|
|
return
|
|
}
|
|
|
|
// dbzColor maps dBZ to RGBA. Below 20 dBZ is transparent — removes AP,
|
|
// biological returns, and ground clutter while keeping real precipitation.
|
|
func dbzColor(dbz float32) color.RGBA {
|
|
switch {
|
|
case dbz < 20:
|
|
return color.RGBA{0, 0, 0, 0}
|
|
case dbz < 25:
|
|
return color.RGBA{0x02, 0x8e, 0x00, 210}
|
|
case dbz < 30:
|
|
return color.RGBA{0x01, 0xb4, 0x4c, 210}
|
|
case dbz < 35:
|
|
return color.RGBA{0x9c, 0xe4, 0x00, 220}
|
|
case dbz < 40:
|
|
return color.RGBA{0xd8, 0xd8, 0x00, 220}
|
|
case dbz < 45:
|
|
return color.RGBA{0xff, 0xaa, 0x00, 230}
|
|
case dbz < 50:
|
|
return color.RGBA{0xff, 0x00, 0x00, 230}
|
|
case dbz < 55:
|
|
return color.RGBA{0xd0, 0x00, 0x00, 240}
|
|
case dbz < 60:
|
|
return color.RGBA{0xc0, 0x00, 0x80, 240}
|
|
default:
|
|
return color.RGBA{0xff, 0xf7, 0x00, 255}
|
|
}
|
|
}
|
|
|
|
// velColor maps radial velocity (m/s) to RGBA.
|
|
// Negative = toward radar (green), positive = away (red).
|
|
func velColor(ms float32) color.RGBA {
|
|
switch {
|
|
case ms <= -27:
|
|
return color.RGBA{0xff, 0x00, 0x00, 230} // strong inbound
|
|
case ms <= -14:
|
|
return color.RGBA{0xff, 0x66, 0x66, 210}
|
|
case ms <= -1:
|
|
return color.RGBA{0xff, 0xaa, 0xaa, 180}
|
|
case ms < 1:
|
|
return color.RGBA{0, 0, 0, 0} // near-zero
|
|
case ms < 14:
|
|
return color.RGBA{0xaa, 0xff, 0xaa, 180}
|
|
case ms < 27:
|
|
return color.RGBA{0x00, 0xcc, 0x00, 210}
|
|
default:
|
|
return color.RGBA{0x00, 0x66, 0x00, 230} // strong outbound
|
|
}
|
|
}
|