kopia lustrzana https://github.com/cyoung/stratux
354 wiersze
8.3 KiB
Go
354 wiersze
8.3 KiB
Go
package main
|
|
|
|
import (
|
|
"fmt"
|
|
"../uatparse"
|
|
"strconv"
|
|
"os"
|
|
"bufio"
|
|
"github.com/kellydunn/golang-geo"
|
|
"math"
|
|
)
|
|
|
|
|
|
// Most adapted from extract_nexrad.c
|
|
|
|
const (
|
|
BLOCK_WIDTH = float64(48.0/60.0)
|
|
WIDE_BLOCK_WIDTH = float64(96.0/60.0)
|
|
BLOCK_HEIGHT = float64(4.0/60.0)
|
|
BLOCK_THRESHOLD = 405000
|
|
BLOCKS_PER_RING = 450
|
|
|
|
WARN_DIST = float64(18.52) // kilometers (10 nm).
|
|
)
|
|
|
|
type NEXRADFrame struct {
|
|
radar_type uint32
|
|
ts string
|
|
scale int
|
|
latNorth float64
|
|
lonWest float64
|
|
height float64
|
|
width float64
|
|
intensity []uint8 // Really only 4-bit values.
|
|
}
|
|
|
|
|
|
func block_location(block_num int, ns_flag bool, scale_factor int) (float64, float64, float64, float64) {
|
|
var realScale float64
|
|
if scale_factor == 1 {
|
|
realScale = float64(5.0)
|
|
} else if scale_factor == 2 {
|
|
realScale = float64(9.0)
|
|
} else {
|
|
realScale = float64(1.0)
|
|
}
|
|
|
|
if block_num >= BLOCK_THRESHOLD {
|
|
block_num = block_num & ^1
|
|
}
|
|
|
|
raw_lat := float64(BLOCK_HEIGHT * float64(int(float64(block_num) / float64(BLOCKS_PER_RING))))
|
|
raw_lon := float64(block_num % BLOCKS_PER_RING) * BLOCK_WIDTH
|
|
|
|
var lonSize float64
|
|
if block_num >= BLOCK_THRESHOLD {
|
|
lonSize = WIDE_BLOCK_WIDTH * realScale
|
|
} else {
|
|
lonSize = BLOCK_WIDTH * realScale
|
|
}
|
|
|
|
latSize := BLOCK_HEIGHT * realScale
|
|
|
|
if ns_flag { // Southern hemisphere.
|
|
raw_lat = 0 - raw_lat
|
|
} else {
|
|
raw_lat = raw_lat + BLOCK_HEIGHT
|
|
}
|
|
|
|
if raw_lon > 180.0 {
|
|
raw_lon = raw_lon - 360.0
|
|
}
|
|
|
|
return raw_lat, raw_lon, latSize, lonSize
|
|
|
|
}
|
|
|
|
func decode_nexrad(f *uatparse.UATFrame) []NEXRADFrame {
|
|
ret := make([]NEXRADFrame, 0)
|
|
|
|
rle_flag := (uint32(f.FISB_data[0]) & 0x80) != 0
|
|
ns_flag := (uint32(f.FISB_data[0]) & 0x40) != 0
|
|
block_num := ((int(f.FISB_data[0]) & 0x0f) << 16) | (int(f.FISB_data[1]) << 8) | (int(f.FISB_data[2]))
|
|
scale_factor := (int(f.FISB_data[0]) & 0x30) >> 4
|
|
|
|
if rle_flag { // Single bin, RLE encoded.
|
|
lat, lon, h, w := block_location(block_num, ns_flag, scale_factor)
|
|
var tmp NEXRADFrame
|
|
tmp.radar_type = f.Product_id
|
|
tmp.ts = strconv.Itoa(int(f.FISB_hours)) + ":" + strconv.Itoa(int(f.FISB_minutes))
|
|
tmp.scale = scale_factor
|
|
tmp.latNorth = lat
|
|
tmp.lonWest = lon
|
|
tmp.height = h
|
|
tmp.width = w
|
|
tmp.intensity = make([]uint8, 0)
|
|
|
|
intensityData := f.FISB_data[3:]
|
|
for _, v := range intensityData {
|
|
intensity := uint8(v) & 0x7;
|
|
runlength := (uint8(v) >> 3) + 1
|
|
for runlength > 0 {
|
|
tmp.intensity = append(tmp.intensity, intensity)
|
|
runlength--
|
|
}
|
|
}
|
|
ret = append(ret, tmp)
|
|
} else {
|
|
var row_start int
|
|
var row_size int
|
|
if block_num >= 405000 {
|
|
row_start = block_num - ((block_num - 405000) % 225)
|
|
row_size = 225
|
|
} else {
|
|
row_start = block_num - (block_num % 450)
|
|
row_size = 450
|
|
}
|
|
|
|
row_offset := block_num - row_start
|
|
|
|
L := int(f.FISB_data[3] & 15)
|
|
for i := 0; i < L; i++ {
|
|
var bb int
|
|
if i == 0 {
|
|
bb = (int(f.FISB_data[3]) & 0xF0) | 0x08
|
|
} else {
|
|
bb = int(f.FISB_data[i+3])
|
|
}
|
|
|
|
for j := 0; j < 8; j++ {
|
|
if bb & (1 << uint(j)) != 0 {
|
|
row_x := (row_offset + 8*i + j - 3) % row_size
|
|
bn := row_start + row_x
|
|
lat, lon, h, w := block_location(bn, ns_flag, scale_factor)
|
|
var tmp NEXRADFrame
|
|
tmp.radar_type = f.Product_id
|
|
tmp.ts = strconv.Itoa(int(f.FISB_hours)) + ":" + strconv.Itoa(int(f.FISB_minutes))
|
|
tmp.scale = scale_factor
|
|
tmp.latNorth = lat
|
|
tmp.lonWest = lon
|
|
tmp.height = h
|
|
tmp.width = w
|
|
tmp.intensity = make([]uint8, 0)
|
|
for k := 0; k < 128; k++ {
|
|
z := uint8(0)
|
|
if f.Product_id == 64 {
|
|
z = 1
|
|
}
|
|
tmp.intensity = append(tmp.intensity, z)
|
|
}
|
|
ret = append(ret, tmp)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return ret
|
|
}
|
|
|
|
func parseInput(buf string) []NEXRADFrame {
|
|
ret := make([]NEXRADFrame, 0)
|
|
|
|
uatmsg, err := uatparse.New(buf)
|
|
if err != nil {
|
|
return ret
|
|
}
|
|
|
|
uatmsg.DecodeUplink()
|
|
|
|
for _, uatframe := range uatmsg.Frames {
|
|
if uatframe.Product_id != 63 && uatframe.Product_id != 64 { // It's neither CONUS nor Regional NEXRAD.
|
|
continue
|
|
}
|
|
p := decode_nexrad(uatframe)
|
|
if p != nil {
|
|
ret = append(ret, p...)
|
|
}
|
|
}
|
|
return ret
|
|
}
|
|
|
|
// Range is 0 to 360.
|
|
func oclock(ang float64) uint8 {
|
|
if ang > 345 || ang <= 15 {
|
|
return 12
|
|
} else if ang > 15 && ang <= 45 {
|
|
return 1
|
|
} else if ang > 45 && ang <= 75 {
|
|
return 2
|
|
} else if ang > 75 && ang <= 105 {
|
|
return 3
|
|
} else if ang > 105 && ang <= 135 {
|
|
return 4
|
|
} else if ang > 135 && ang <= 165 {
|
|
return 5
|
|
} else if ang > 165 && ang <= 195 {
|
|
return 6
|
|
} else if ang > 195 && ang <= 225 {
|
|
return 7
|
|
} else if ang > 225 && ang <= 255 {
|
|
return 8
|
|
} else if ang > 255 && ang <= 285 {
|
|
return 9
|
|
} else if ang > 285 && ang <= 315 {
|
|
return 10
|
|
} else if ang > 315 && ang <= 345 {
|
|
return 11
|
|
}
|
|
return 0
|
|
}
|
|
|
|
func intensityToText(intensity uint8) string {
|
|
if intensity >= 0 && intensity < 3 {
|
|
return "light"
|
|
} else if intensity >= 3 && intensity < 6 {
|
|
return "moderate"
|
|
} else if intensity == 6 {
|
|
return "heavy"
|
|
} else if intensity == 7 {
|
|
return "very heavy"
|
|
}
|
|
return ""
|
|
}
|
|
|
|
func fixHeading(hdg float64) float64 {
|
|
if hdg < 0 {
|
|
return float64(hdg + 360)
|
|
}
|
|
if hdg >= 360 {
|
|
return float64(hdg - 360)
|
|
}
|
|
return float64(hdg)
|
|
}
|
|
|
|
func scanNEXRAD(poly *geo.Polygon, frame NEXRADFrame) (*geo.Point, uint8) {
|
|
var retpt *geo.Point
|
|
var maxIntensity uint8
|
|
for y := 0; y < 4; y++ {
|
|
for x := 0; x < 32; x++ {
|
|
intensity := frame.intensity[x + 32*y]
|
|
lat := frame.latNorth - (float64(y) * (frame.height)/float64(4.0))
|
|
lon := frame.lonWest + (float64(x) * (frame.width)/float64(32.0))
|
|
pt := geo.NewPoint(lat, lon)
|
|
if !poly.Contains(pt) { // Doesn't contain this point - skip.
|
|
continue
|
|
}
|
|
if intensity > maxIntensity {
|
|
retpt = pt
|
|
maxIntensity = intensity
|
|
}
|
|
}
|
|
}
|
|
return retpt, maxIntensity
|
|
}
|
|
|
|
func main() {
|
|
if len(os.Args) < 5 {
|
|
fmt.Printf("%s <uat log> <lat> <lon> <hdg>\n", os.Args[0])
|
|
return
|
|
}
|
|
|
|
fd, err := os.Open(os.Args[1])
|
|
|
|
if err != nil {
|
|
fmt.Printf("can't open file: %s\n", err.Error())
|
|
return
|
|
}
|
|
|
|
hdg, err := strconv.Atoi(os.Args[4])
|
|
if err != nil || hdg > 360 || hdg < 0 {
|
|
fmt.Printf("invalid heading: %s\n", os.Args[4])
|
|
return
|
|
}
|
|
lat, err := strconv.ParseFloat(os.Args[2], 64)
|
|
if err != nil {
|
|
fmt.Printf("invalid lat: %s\n", os.Args[2])
|
|
return
|
|
}
|
|
lon, err := strconv.ParseFloat(os.Args[3], 64)
|
|
if err != nil {
|
|
fmt.Printf("invalid lon: %s\n", os.Args[3])
|
|
return
|
|
}
|
|
|
|
hdgFloat := float64(hdg)
|
|
|
|
frames := make([]NEXRADFrame, 0)
|
|
|
|
reader := bufio.NewReader(fd)
|
|
|
|
for {
|
|
buf, err := reader.ReadString('\n')
|
|
if err != nil {
|
|
break
|
|
}
|
|
o := parseInput(buf)
|
|
frames = append(frames, o...)
|
|
}
|
|
|
|
// Do processing.
|
|
myPos := geo.NewPoint(lat, lon)
|
|
// We'll now draw a rectangle 20nm wide by 10nm tall in space, with the aircraft at the center of the bottom edge.
|
|
// This gives 180 degrees of "visibility" for a decent sized area.
|
|
nineOClock := fixHeading(hdgFloat - 90.0)
|
|
threeOClock := fixHeading(hdgFloat + 90.0)
|
|
|
|
// fmt.Printf("myPos=%v\n", myPos)
|
|
|
|
leftBottom := myPos.PointAtDistanceAndBearing(WARN_DIST, nineOClock)
|
|
rightBottom := myPos.PointAtDistanceAndBearing(WARN_DIST, threeOClock)
|
|
|
|
// fmt.Printf("nineOClock=%f [leftBottom=%v], threeOClock=%f [rightBottom=%v]\n", nineOClock, leftBottom, threeOClock, rightBottom)
|
|
|
|
hypDist := math.Sqrt2 * WARN_DIST
|
|
leftTopHdg := fixHeading(hdgFloat - 45.0)
|
|
rightTopHdg := fixHeading(hdgFloat + 45.0)
|
|
|
|
leftTop := myPos.PointAtDistanceAndBearing(hypDist, leftTopHdg)
|
|
rightTop := myPos.PointAtDistanceAndBearing(hypDist, rightTopHdg)
|
|
|
|
// fmt.Printf("leftTopHdg=%f [leftTop=%v], rightTopHdg=%f [rightTop=%v]\n", leftTopHdg, leftTop, rightTopHdg, rightTop)
|
|
|
|
points := []*geo.Point{leftTop, rightTop, rightBottom, leftBottom, leftTop}
|
|
poly := geo.NewPolygon(points)
|
|
|
|
|
|
var maxpt *geo.Point
|
|
var maxIntensity uint8
|
|
|
|
for _, frame := range frames {
|
|
//FIXME: Scans the whole map.
|
|
thisMaxpt, thisMaxIntensity := scanNEXRAD(poly, frame)
|
|
if thisMaxIntensity > maxIntensity {
|
|
maxpt = thisMaxpt
|
|
maxIntensity = thisMaxIntensity
|
|
}
|
|
}
|
|
|
|
// fmt.Printf("maxes: %d %v\n", maxIntensity, maxpt)
|
|
|
|
|
|
if maxIntensity > 0 && maxpt != nil {
|
|
desc := intensityToText(maxIntensity)
|
|
direction := fixHeading(myPos.BearingTo(maxpt))
|
|
relativeDirection := fixHeading(direction - hdgFloat)
|
|
// fmt.Printf("direction=%f, relativeDirection=%f\n", direction, relativeDirection)
|
|
directionDesc := oclock(relativeDirection)
|
|
dist := myPos.GreatCircleDistance(maxpt) * float64(0.539957) // Convert km -> nm.
|
|
fmt.Printf("%s precip %d o'clock, %0.1f nm.\n", desc, directionDesc, dist)
|
|
} else {
|
|
fmt.Printf("no precip.\n")
|
|
}
|
|
|
|
}
|