231 lines
6.2 KiB
Go
231 lines
6.2 KiB
Go
package radar
|
|
|
|
import (
|
|
"bytes"
|
|
"compress/bzip2"
|
|
"encoding/binary"
|
|
"fmt"
|
|
"io"
|
|
"math"
|
|
"time"
|
|
)
|
|
|
|
// RadarProduct holds parsed data from one Level 3 scan.
|
|
type RadarProduct struct {
|
|
Site string
|
|
Product string // "reflectivity" or "velocity"
|
|
SiteLat float64
|
|
SiteLon float64
|
|
Elevation float64
|
|
Time time.Time
|
|
Radials []Radial
|
|
}
|
|
|
|
// Radial is one spoke of radar data.
|
|
type Radial struct {
|
|
Azimuth float64
|
|
DeltaAz float64
|
|
RangeKm float64
|
|
BinSizeKm float64
|
|
Values []float32 // dBZ per bin; NaN = no data
|
|
}
|
|
|
|
// GateLL returns the lat/lon center of range bin i in this radial.
|
|
func (r *Radial) GateLL(siteLat, siteLon float64, i int) (lat, lon float64) {
|
|
distKm := r.RangeKm + float64(i)*r.BinSizeKm + r.BinSizeKm/2
|
|
return rangeBearing(siteLat, siteLon, distKm, r.Azimuth)
|
|
}
|
|
|
|
// Parse decodes a NEXRAD Level 3 base reflectivity product (N0Q / product 94).
|
|
func Parse(site string, data []byte) (*RadarProduct, error) {
|
|
return parseProduct(site, "reflectivity", data, func(v byte) float32 {
|
|
if v <= 1 {
|
|
return float32(math.NaN())
|
|
}
|
|
return float32(v-2)*0.5 - 32.0 // dBZ
|
|
})
|
|
}
|
|
|
|
// ParseVelocity decodes a NEXRAD Level 3 base velocity product (N0U / product 99).
|
|
func ParseVelocity(site string, data []byte) (*RadarProduct, error) {
|
|
return parseProduct(site, "velocity", data, func(v byte) float32 {
|
|
if v <= 1 {
|
|
return float32(math.NaN())
|
|
}
|
|
return float32(v-2)*0.5 - 63.5 // m/s
|
|
})
|
|
}
|
|
|
|
func parseProduct(site, productType string, data []byte, decode func(byte) float32) (*RadarProduct, error) {
|
|
data = stripWMOHeader(data)
|
|
if len(data) < 120 {
|
|
return nil, fmt.Errorf("file too short after header strip (%d bytes)", len(data))
|
|
}
|
|
|
|
// --- Message Header Block (18 bytes) — skip ---
|
|
|
|
// --- Product Description Block (starts at byte 18) ---
|
|
pdb := data[18:]
|
|
if len(pdb) < 102 {
|
|
return nil, fmt.Errorf("product description block truncated")
|
|
}
|
|
|
|
siteLat := float64(int32(binary.BigEndian.Uint32(pdb[2:]))) / 1000.0
|
|
siteLon := float64(int32(binary.BigEndian.Uint32(pdb[6:]))) / 1000.0
|
|
volDate := int16(binary.BigEndian.Uint16(pdb[22:]))
|
|
volTimeSec := int32(binary.BigEndian.Uint32(pdb[24:])) // seconds since midnight
|
|
scanTime := julianToTime(volDate, volTimeSec)
|
|
|
|
// data[120:] is the Symbology Block, bzip2-compressed for digital products (N0Q).
|
|
// 18 bytes Message Header + 102 bytes PDB = 120 bytes prefix.
|
|
symData := data[120:]
|
|
if len(symData) >= 3 && symData[0] == 'B' && symData[1] == 'Z' && symData[2] == 'h' {
|
|
r := bzip2.NewReader(bytes.NewReader(symData))
|
|
decompressed, err := io.ReadAll(r)
|
|
if err != nil {
|
|
return nil, fmt.Errorf("decompress symbology: %w", err)
|
|
}
|
|
symData = decompressed
|
|
}
|
|
|
|
// Find the Symbology Block by scanning for its signature:
|
|
// block divider (FF FF) + block ID (00 01).
|
|
symPos := -1
|
|
for i := 0; i+3 < len(symData); i++ {
|
|
if symData[i] == 0xFF && symData[i+1] == 0xFF && symData[i+2] == 0x00 && symData[i+3] == 0x01 {
|
|
symPos = i
|
|
break
|
|
}
|
|
}
|
|
if symPos < 0 {
|
|
return nil, fmt.Errorf("symbology block not found in decompressed data")
|
|
}
|
|
|
|
// --- Symbology Block header (10 bytes) ---
|
|
sym := symData[symPos:]
|
|
// [0-1]: block divider (-1)
|
|
// [2-3]: block ID (1)
|
|
// [4-7]: block length (bytes)
|
|
// [8-9]: number of layers
|
|
|
|
// --- Layer header (6 bytes) ---
|
|
if 10+6 > len(sym) {
|
|
return nil, fmt.Errorf("symbology block too short for layer")
|
|
}
|
|
layer := sym[10:]
|
|
// [0-1]: layer divider (-1)
|
|
// [2-5]: layer length (bytes)
|
|
|
|
pkt := layer[6:]
|
|
packetCode := int16(binary.BigEndian.Uint16(pkt[0:]))
|
|
if packetCode != 16 {
|
|
return nil, fmt.Errorf("unsupported packet code %d (want 16)", packetCode)
|
|
}
|
|
|
|
// --- Packet Code 16 header (14 bytes) ---
|
|
// [0-1]: code (16)
|
|
// [2-3]: first range bin index
|
|
// [4-5]: number of range bins
|
|
// [6-7]: I center
|
|
// [8-9]: J center
|
|
// [10-11]: scale factor (thousandths of km per pixel)
|
|
// [12-13]: number of radials
|
|
numBins := int(int16(binary.BigEndian.Uint16(pkt[4:])))
|
|
scaleFactor := int(int16(binary.BigEndian.Uint16(pkt[10:])))
|
|
numRadials := int(int16(binary.BigEndian.Uint16(pkt[12:])))
|
|
|
|
binSizeKm := 1.0
|
|
if scaleFactor > 0 {
|
|
binSizeKm = float64(scaleFactor) / 1000.0
|
|
}
|
|
|
|
pos := 14 // offset into pkt
|
|
radials := make([]Radial, 0, numRadials)
|
|
for i := 0; i < numRadials; i++ {
|
|
if pos+6 > len(pkt) {
|
|
break
|
|
}
|
|
runBytes := int(int16(binary.BigEndian.Uint16(pkt[pos:])))
|
|
startAngle := float64(int16(binary.BigEndian.Uint16(pkt[pos+2:]))) / 10.0
|
|
deltaAngle := float64(int16(binary.BigEndian.Uint16(pkt[pos+4:]))) / 10.0
|
|
pos += 6
|
|
|
|
if pos+runBytes > len(pkt) {
|
|
break
|
|
}
|
|
raw := pkt[pos : pos+runBytes]
|
|
pos += runBytes
|
|
// Pad to even-byte boundary.
|
|
if runBytes%2 != 0 {
|
|
pos++
|
|
}
|
|
|
|
values := make([]float32, numBins)
|
|
for j := 0; j < numBins && j < len(raw); j++ {
|
|
values[j] = decode(raw[j])
|
|
}
|
|
|
|
radials = append(radials, Radial{
|
|
Azimuth: startAngle,
|
|
DeltaAz: deltaAngle,
|
|
RangeKm: 0.0,
|
|
BinSizeKm: binSizeKm,
|
|
Values: values,
|
|
})
|
|
}
|
|
|
|
if len(radials) == 0 {
|
|
return nil, fmt.Errorf("no radials parsed")
|
|
}
|
|
|
|
return &RadarProduct{
|
|
Site: site,
|
|
Product: productType,
|
|
SiteLat: siteLat,
|
|
SiteLon: siteLon,
|
|
Elevation: 0.5,
|
|
Time: scanTime,
|
|
Radials: radials,
|
|
}, nil
|
|
}
|
|
|
|
// stripWMOHeader removes the variable-length WMO/AFOS text header that NWS
|
|
// prepends to distributed products. The binary NEXRAD data begins after the
|
|
// final \r\r\n sequence in the first 100 bytes.
|
|
func stripWMOHeader(data []byte) []byte {
|
|
limit := 100
|
|
if len(data) < limit {
|
|
limit = len(data)
|
|
}
|
|
last := 0
|
|
for i := 0; i+2 < limit; i++ {
|
|
if data[i] == '\r' && data[i+1] == '\r' && data[i+2] == '\n' {
|
|
last = i + 3
|
|
}
|
|
}
|
|
return data[last:]
|
|
}
|
|
|
|
func julianToTime(julianDate int16, secs int32) time.Time {
|
|
base := time.Date(1970, 1, 1, 0, 0, 0, 0, time.UTC)
|
|
d := base.AddDate(0, 0, int(julianDate)-1)
|
|
return d.Add(time.Duration(secs) * time.Second)
|
|
}
|
|
|
|
func rangeBearing(lat0, lon0, distKm, bearing float64) (lat, lon float64) {
|
|
const earthR = 6371.0
|
|
lat0r := lat0 * math.Pi / 180
|
|
lon0r := lon0 * math.Pi / 180
|
|
br := bearing * math.Pi / 180
|
|
dr := distKm / earthR
|
|
|
|
latr := math.Asin(math.Sin(lat0r)*math.Cos(dr) +
|
|
math.Cos(lat0r)*math.Sin(dr)*math.Cos(br))
|
|
lonr := lon0r + math.Atan2(
|
|
math.Sin(br)*math.Sin(dr)*math.Cos(lat0r),
|
|
math.Cos(dr)-math.Sin(lat0r)*math.Sin(latr),
|
|
)
|
|
return latr * 180 / math.Pi, lonr * 180 / math.Pi
|
|
}
|
|
|