#!/usr/bin/env python
# usage: plotrean.py filefromgrib.dat
# or just: plotrean.py
# This program makes PPM plots of NCEP/NCAR Reananalysis Data, extracted into text files
# To get an example file to play with, enter the following from the Linux command line:
# wgrib /REANALYSIS/1999/data/at00z12z/pres.sfc -d 1 -text -o pres1.dat
# This will extract the first record from 1999 in pres.sfc to the text file pres1.dat.
# (You are invited to snoop around in /REANALYSIS to see the data we have.)

from os import system #needed for the final statement
import sys, math 
#the program needs the name of the text file with the data:
try:
    infilename=sys.argv[1] #this works if you can specify the file name from the command line
except:
    infilename=raw_input('enter data file name==> ') #otherwise you are prompted for it
inf=open(infilename,'r')  #opens the file for reading
all=inf.readlines() #read all the lines into a array
header=all.pop(0).split() #remove the first line (a string) and split into substrings at the blanks
w=int(header[0]) #convert the first substring into an integer
h=int(header[1]) #convert the second substring into an integer
print "will assume w=",w,"   h=",h   
ouf=open('myplot.ppm','w') #open a PPM file in which to write the RGB triplets
#write header into PPM file:
ouf.write("P3\n") #you can use P3 for ascii, P6 for binary
ouf.write("%d %d\n" % (w, h)) 
ouf.write("255\n")
nc=0 #initialize counter
maxp=-1.e30 #a very negative number
minp=1.e30 #a very positive number
for s in all: #find maximum and minimum of data values
	p=float(s) #s is a string, and must be converted
	maxp=max(p,maxp)
	minp=min(p,minp)
#uncomment the next two lines to manually set the extremes that bound the color scale:
#maxp=105000.
#minp=95000.
for s in all:
	p=float(s) #somewhat inefficent, because we did this already...what the heck...
	v=(maxp-p)/(maxp-minp) 
#	v=math.sqrt(v) #uncomment this to accentuate smaller v (higher pressure)
	red=min(255,int(255*v))  #safety filter: don't let red exceed 255
	red=max(0,red) #safety filter: don't let red be less than 0
	blue=255-red 
	green=0
	ouf.write("%d %d %d\n" % (red,green,blue)) #use with P3
#	ouf.write("%c%c%c" % (red,green,blue)) #use with P6, PPM file will be smaller
	nc+=1
ouf.close()
print "file myplot.ppm was written"
print  "Here are the extremes on the colorscale:  max=",maxp,"  min=",minp
print h*w," triplets were expected, and ",nc," triplets were written"
if h*w==nc:
	print "Wonderful, that is what was supposed to happen!"
else:
	print "Oops.  Those numbers should have been equal."
system("display myplot.ppm") #optional, if want Linux to display it automatically
#
# Notes for the task for Fall 2005:
# The pres.sfc data is dominated by topography.  You might think that
# if we focus on the ocean surface, where there is no topography, that
# we should see "weather".  But oddly, this is not the case.  Instead
# we see a lot of numerical artifacts (Gibbs effect) which you can read
# about by typing:
# cat /REANALYSIS/1999/text/details/pres.sfc
# To see this nasty Gibbs effect better, modify this program to colorize
# only the points that are over the ocean.  Make the land points black.
# First, you need to get the land mask file:
# wgrib /REANALYSIS/1999/data/fixed/land.sfc -d 1 -text -o land.dat
#
# Hint (1): Suppose you have read the lines of land.dat into an array landmask.
# (And popped off the header).  Change the "for loops" to something like:
# for i in range(0,h*w):
#	p=float(all[i])
# 	l=int(landmask[i])
#   .....
# Hint (2): Beware of the high-altitude lake in Africa!  You may need to
# manually set maxp and minp to maximally accentuate the Gibbs effect.