Search for precip in local area, translate in to text warning.

$ ./nexrad_annunciator ../test-data/example.radar 38.32 -88.86 47
moderate precip 1 o'clock, 5.2 nm.
$ ./nexrad_annunciator ../test-data/example.radar 38.38 -88.41 140
heavy precip 2 o'clock, 8.9 nm.
$ ./nexrad_annunciator ../test-data/example.radar 38.22 -88.25 140
very heavy precip 11 o'clock, 6.7 nm.

test/nexrad_annunciator.go

@@ -6,6 +6,8 @@ import (
 	"strconv"
 	"os"
 	"bufio"
+	"github.com/kellydunn/golang-geo"
+	"math"
 	)
 
 
@@ -17,6 +19,8 @@ const (
 	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 {
@@ -62,10 +66,10 @@ func block_location(block_num int, ns_flag bool, scale_factor int) (float64, flo
 	} 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
 
@@ -174,10 +178,83 @@ func parseInput(buf string) []NEXRADFrame {
 	return ret
 }
 
-func main() {
+// 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
+}
 
-	if len(os.Args) < 2 {
-		fmt.Printf("%s <uat log>\n", os.Args[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
 	}
 
@@ -188,26 +265,89 @@ func main() {
 		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 {
-			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")
+		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")
+	}
+
 }