#!/usr/bin/env python
#usage: rasterradar.py KIWA.vol140.dat
# v0.24, 8 July 2011
import sys
from math import *

### MAKE THE COLOR CODES FOR PPM FORMAT
# The following are the positive colors in the legend for the RAP radar images
# at http://www.rap.ucar.edu/weather/radar/ :
hexc="fff 93c f0f 900 c30 f00 f60 fc0 ff0 090 0c0 0f0 00f 39f 6ff 000".split()
hexc.reverse() # flip order to be from low to high 
levels=range(0,76,5) # the lower bound to the the dBZ color
rgb=[ ] # will contain 3-byte strings of the color 
# colors in hexc need to be expanded, e.g. fc0 => ffcc00
# tip: list(string) puts characters into a list, 2*s doubles a string
# tip: int('cc',16) => 204, the 16 is the base of the number system
# tip: also see "List comprehensions" at http://docs.python.org/tut/node7.html
for abc in hexc: # loop over all color codes, make new codes
	triplet=[int(2*s,16) for s in list(abc)] # e.g, 'fc0' => [255, 204, 0]
	c3="%c%c%c" % tuple(triplet) # 3 bytes, but most of the characters do not print
	print "converted ",abc," to ",triplet," and ",c3
	rgb.append(c3) #store the color as 3 bytes

### MAKE A VERY TINY PPM OF THE RADAR COLOR LEGEND
ouc=open("radarcolors.ppm",'wb')
# write PPM header:
ouc.write('P6\n')
ouc.write('#radar colors\n')
ouc.write("1 %d\n" % len(rgb))
ouc.write("255\n")
rgb.reverse() # flip order to print high values at top of image
ouc.write("".join(rgb)) #this trick concatenates a list of strings
rgb.reverse() # flip back to low to high
ouc.close()
print "wrote radarcolors.ppm" 

### OPEN THE DATA FILE, AND SAMPLE IT	
try: # try to get file name from command line value
	filename=sys.argv[1] 
except: # if not successful, prompt user to enter it
	filename=raw_input("enter file name of radar data=>")
inf=open(filename,'r')
header=inf.readline().split() # read first line, and split at white space
scale=header[20] # this prints 100, but I think we should ignore it
missing=header[21] # missing value, but we simply replace all negatives as 0
datalines=inf.readlines() # read all the rest of the lines into a list, one line per item
w=len(datalines[0].split()) # number of angular measurements per line
h=len(datalines) # number of lines of data, of increasing radius (number of "range gates")
print "sampled "+filename+" :   scale=",scale,"   missing=",missing,"  num of angl=",w,"  num ranges=",h

### PLOT DATA IN THETA,R  PLANE 
# rasterize as rows and columns, ignore polar coordinates
oup=open("radarimage.ppm",'wb')
# write PPM header:
oup.write('P6\n')
oup.write('#radar image\n')
oup.write("%d %d\n" % (w,h)) 
oup.write("255\n")
datart=[] # will be list of lists of data (colorized), for your task 
print "reading data from "+filename
# tip: also see "List comprehensions" at http://docs.python.org/tut/node7.html
for line in datalines:
	v=[max(float(x),0.) for x in line.split()] # convert 360 items to floating point values
	c=[rgb[int(x/5.)] for x in v] # convert floats to integer index, then to color
	datart.append(c) #store list of colorized row of data in datart (use datart in your task)
	oup.write("".join(c))
oup.close()
print "wrote radarimage.ppm" 
 
### TASK:  PLOT THE DATA IN POLAR COORDINATES  
print "plotting in polar coordinates..."
print "sample color in r,theta array:", datart[10][359]
xras=600 # your choice for pixel width of raster image
yras=600 # your choice for pixel height of raster image
xc=xras/2 # center x pixel
yc=yras/2 # center y pixel  
# write PPM header:
oup=open("taskimage.ppm",'wb')
oup.write('P6\n')
oup.write('#radar image\n')
oup.write("%d %d\n" % (xras,yras)) 
oup.write("255\n")
outcolor="%c%c%c" % (100,100,100) # define a gray color for outside the range
# loop thru you pixels, finding the r and theta value of the piexel:
for j in range(yras):
	for i in range(xras):
		x=+(i+.5-xc) # x distance of pixel center from radar 
		y=-(j+.5-yc) # y distance of pixel center from radar
		r=sqrt(x*x+y*y) # distance from radar
		t=atan2(x,y)*180/3.1415 # azimuth, 0 deg is North, 90 deg is East
### approximately 7 lines needed here to complete the 1st part of task!
# your 7 lines here:

# end your code.
oup.close()		
print "wrote taskimage.ppm" 
## 2nd part of the task is defined on the web page