phase 0-2 complete
This commit is contained in:
Blake Ridgway
230
internal/radar/parse.go
Normal file
230
internal/radar/parse.go
Normal file
@@ -0,0 +1,230 @@
|
||||
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
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user