From f1ea6bd1a6c76505456da0326e7ca522731d5dd2 Mon Sep 17 00:00:00 2001 From: Christopher Young Date: Tue, 22 Sep 2015 19:40:29 -0400 Subject: [PATCH] Working NEXRAD decode. --- test/nexrad_annunciator.go | 213 +++++++++++++++++++++++++++++++++++++ test/uatparse/uatparse.go | 28 ++--- 2 files changed, 229 insertions(+), 12 deletions(-) create mode 100644 test/nexrad_annunciator.go diff --git a/test/nexrad_annunciator.go b/test/nexrad_annunciator.go new file mode 100644 index 00000000..34319c5e --- /dev/null +++ b/test/nexrad_annunciator.go @@ -0,0 +1,213 @@ +package main + +import ( + "fmt" + "./uatparse" + "strconv" + "os" + "bufio" + ) + + +// 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 +) + +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 +} + +func main() { + + if len(os.Args) < 2 { + fmt.Printf("%s \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 + } + + reader := bufio.NewReader(fd) + + for { + buf, err := reader.ReadString('\n') + if err != nil { + fmt.Printf("lost stdin.\n") + break + } + z := parseInput(buf) + for _, zz := range z { + n := "Regional" + if zz.radar_type == 64 { + n = "CONUS" + } + fmt.Printf("NEXRAD %s %s %d %.0f %.0f %.0f %.0f ", n, zz.ts, zz.scale, zz.latNorth * 60, zz.lonWest * 60, zz.height * 60, zz.width * 60) + for _, intens := range zz.intensity { + fmt.Printf("%d", intens) + } + fmt.Printf("\n") + } + } + +} \ No newline at end of file diff --git a/test/uatparse/uatparse.go b/test/uatparse/uatparse.go index d06ecf76..0645411a 100644 --- a/test/uatparse/uatparse.go +++ b/test/uatparse/uatparse.go @@ -27,7 +27,10 @@ const ( ) type UATFrame struct { - raw_data []byte + Raw_data []byte + FISB_data []byte + FISB_hours uint32 + FISB_minutes uint32 frame_length uint32 Frame_type uint32 @@ -81,11 +84,7 @@ func dlac_decode(data []byte, data_len uint32) string { func (f *UATFrame) decodeInfoFrame() { - f.Product_id = ((uint32(f.raw_data[0]) & 0x1f) << 6) | (uint32(f.raw_data[1]) >> 2) - - if f.Product_id != 413 { // FIXME. - return - } + f.Product_id = ((uint32(f.Raw_data[0]) & 0x1f) << 6) | (uint32(f.Raw_data[1]) >> 2) if f.Frame_type != 0 { return // Not FIS-B. @@ -94,7 +93,7 @@ func (f *UATFrame) decodeInfoFrame() { return // Too short for FIS-B. } - t_opt := ((uint32(f.raw_data[1]) & 0x01) << 1) | (uint32(f.raw_data[2]) >> 7) + t_opt := ((uint32(f.Raw_data[1]) & 0x01) << 1) | (uint32(f.Raw_data[2]) >> 7) if t_opt != 0 { //FIXME. // fmt.Printf("don't know time format %d\n", t_opt) @@ -102,11 +101,16 @@ func (f *UATFrame) decodeInfoFrame() { return } - /* fisb_hours := (uint32(data[2]) & 0x7c) >> 2 - fisb_minutes := ((uint32(data[2]) & 0x03) << 4) | (uint32(data[3]) >> 4) - */ + f.FISB_hours = (uint32(f.Raw_data[2]) & 0x7c) >> 2 + f.FISB_minutes = ((uint32(f.Raw_data[2]) & 0x03) << 4) | (uint32(f.Raw_data[3]) >> 4) + fisb_length := f.frame_length - 4 - fisb_data := f.raw_data[4:] + fisb_data := f.Raw_data[4:] + f.FISB_data = fisb_data + + if f.Product_id != 413 { // Doesn't have text data that we'd be interested in, so we're done. + return + } p := dlac_decode(fisb_data, fisb_length) ret := make([]string, 0) @@ -178,7 +182,7 @@ func (u *UATMsg) DecodeUplink() error { data = data[2 : frame_length+2] thisFrame := new(UATFrame) - thisFrame.raw_data = data + thisFrame.Raw_data = data thisFrame.frame_length = frame_length thisFrame.Frame_type = frame_type