#!/usr/bin/python

import sys, vplot, getopt

opts, args = getopt.getopt(sys.argv[1:], 'hT:t:H:b:DS:r:s:', 
    ('help','title=', 'type=', 'header=', 'bins=', 'start=', 'radius=', 
     'size=', 'debug')
    )

# define default line numbers to find title, header, bins, and calm value
title = ''
title_line = 1
header_line = 2
bins_line = 3
calm_line = 0
bintype = 'counts'
nodir = 0
dirbins = []
binheaders = []    # speed bin headers
dirheaders = []    # direction headers (parsed but not used for anything?)
calmbin = 0
start_text = ''
max_radius = 0
image_size = '512x512'
debug = 0

colors = [
    (0,.5,0),(0,.7,.7),(.5,.6,.6),(.8,.6,.4),(1,0,0),(1,1,0),(0,1,0),(0,1,1)
    ]

for o, a in opts:

    if o == '-h' or  o == '--help':
        print """
NAME
    windrose-eps - process speed & direction bin data to graphical windrose 

SYNOPSIS
    $ windrose-eps [OPTIONS] infile outfile

DESCRIPTION
    Input file consists of a set of lines representing compass directions,
    and containing event counts for each of a set of windspeed "bins".
    The first line defines the speed bins, each value representing the
    low limit for its bin.  Following the direction and speed data is
    a line containing "calm" event counts, where "calm" is any speed
    below the least speed bin.

    Input and output filenames are required, but either may be given as
    "-" to represent stdin or stdout respectively.

    Graphical output is in eps (encapsulated postscript) format, which can
    be viewed by postscript viewers, included in other postscript pages,
    etc., or converted to other formats, e.g., using ImageMagick convert.

OPTIONS
    -h,  --help      : display this help/manpage info
    -TS, --title=S   : title for the windrose plot
    -TN, --title=N   : read title from line N (default -T0)
    -tS, --type=S    : bin data type is "counts", "frequency", or "percent"
    -HN, --header=N  : speed bins "header" is on line N (default -H1)
    -bN, --bins=N    : bin data begins on line N (default -b2)
    -cN, --calm=N    : calm events data is on line N
    -rN, --radius=N  : max radius as frequency or percent (include "%" sign)
    -D,  --nodir     : no direction identifier is included on data lines
    -SS, --start=S   : start string flagging bins header line
    -sS, --size=S    : image size in pixels, e.g., -s512x512 (default)
    --debug          : display debugging info

    note 1: in the descriptions above, S denotes a string, N a number
    note 2: string argument to -t or --type may be abbreviated
    note 3: line numbers start at 1; line 0 will not match any line
  
BUGS
    Input data format is arbitrary to match other data sets already in use.

AUTHOR
  Ken Irving <fnkci@uaf.edu>
  based on vplot.py and roseplot2.py by Brian Fiedler <brian_fiedler@yahoo.com>

"""
        sys.exit(1) # ... what is normal return value for -h ?

    if o == '-T' or o == '--title':
        if a.isdigit: # title can be string or line number containing title
            title_line = int(a)
        else:
            title = a

    if o == '-t' or o == '--type': 
        if 'counts'.find(a) == 0:
            type = 'counts'
        elif 'frequency'.find(a) == 0:
            type = 'frequency'
        elif 'percent'.find(a) == 0:
            type = 'percent'

    if o == '-H' or o == '--header': 
        header_line = int(a)

    if o == '-b' or o == '--bins': 
        bins_line = int(a)

    if o == '-c' or o == '--calm': 
        calm_line = int(a)

    if o == '-S' or o == '--start': 
        start_text = a

    if o == '-D' or o == '--nodir':
        no_dir_field = 1 

    if o == '-r' or o == '--radius': 
        if a.find("%") == -1:
            max_radius = float(a)
        else:
            max_radius = float(a.rstrip("%")) / 100.0

    if o == '-s' or o == '--size': 
        if a.find("x") == -1:
            sys.stderr.write('invalid -s or --size; see windrose-eps --help')
            sys.exit(1)
        image_size = a

    if o == '--debug':
        debug = 1 

image_x, image_y = image_size.split('x')
image_x = int(image_x)
image_y = int(image_y)

# using a subroutine for constants is a perlish notion...
def RADIUS_RATIO(): return 2.8 # really the inverse ratio...

if image_x <= image_y:
    rsize = image_x / RADIUS_RATIO()
else:
    rsize = image_y / RADIUS_RATIO()

if len(args) != 2: 
    sys.stderr.write('synopsis: windrose-eps [OPTIONS] infile outfile\n')
    sys.exit(1)

ifname = args[0]
if ifname == "-":
    ifname = "STDIN"
    infile = sys.stdin
else:
    infile = open(ifname,'r')

ofname = args[1]

if start_text:  # override other line specifiers...
    header_line = 0
    bins_line = 0
    calm_line = 0

# read the input data file:
numbins = 0 # initialize bin count to dummy value
sumbins = 0
# wf = []
maxsumdir = 0
n = 0

if debug:
    print "DEBUG:\n" \
        "\ttitle_line = ",  title_line, "\n", \
        "\theader_line = ", header_line, "\n", \
        "\tbins_line = ",   bins_line, "\n", \
        "\tcalm_line = ",   calm_line, "\n", \

find_start = 0

while 1:
    record = infile.readline()
    if not record: break
    n += 1

    if debug:
        print n, ':', record,

    if n == title_line:
        title += record
      # if debug: print "title:" + title,
        continue

    if start_text:
        testline = record.split()
        if len(testline) < 2: continue
        if testline[0] == start_text:
            header_line = n
            bins_line = n+1
            start_text = ''
        else:
            continue

    if n == header_line:
        binheaders = record.split()
      # print "DEBUG binheaders:", binheaders
        numbins = len(binheaders)
        if not nodir: dirheaders.append(binheaders.pop(0))
        continue

    if n >= bins_line:
        bins = record.split()
      # print "DEBUG:", bins
        
        if numbins == 0:  # first line of bins, no header
            numbins = len(bins)
            dirbins.append(bins)
        elif len(bins) == numbins:
            if not nodir: dirheaders.append(bins.pop(0)) 
            for i in range(len(bins)): bins[i] = float(bins[i])
            dirbins.append(bins)
            sumdir = reduce(lambda x, y: x+y, bins)
            if debug:
                print "DEBUG: sumdir:", sumdir
            if sumdir > maxsumdir: maxsumdir = sumdir
          # wf.append(sumdir)
            sumbins = sumbins + sumdir
        elif len(bins) == 1 and calm_line == 0:
            calmbin = float(bins[0])
            break
        elif len(bins) == 2 and calm_line == 0:
            calmbin = float(bins[1])
            break
        else:
            sys.stderr.write('mismatched number of bins')
        
sumbins += calmbin

  # # i,bb = eval(aline) #
  # b = reduce(lambda x, y: x+y, bb) # see "help(reduce)"
  # wf.append(b) # list of events in each angle bin 
  # wfs.append(bb) # list of lists of windspeed bins
  # sumb = sumb+b # total of all windspeed bins


if debug > 0:
    print \
        "DEBUG binheaders:\n", binheaders, \
      "\nDEBUG dirheaders:\n", dirheaders, \
      "\nDEBUG dirbins:\n", dirbins, \
      "\nDEBUG calmbin:\n", calmbin, \
      "\nDEBUG sumbins:\n", sumbins,

# sumb = sumb+calm+dirless # total of all events

spd_bin_totals = [] # sum counts for each speed, all directions
for i in range(len(binheaders)): spd_bin_totals.append(0)

for d in range(len(dirbins)):
    for s in range(len(dirbins[d])):
        spd_bin_totals[s] += dirbins[d][s]

if debug:
    print "speed ranges and totals:"
    print "total", sumbins, sumbins / sumbins * 100, "%"
    print "calm", calmbin, calmbin / sumbins * 100, "%"
    for i in range(len(spd_bin_totals)): 
        print i, binheaders[i], spd_bin_totals[i], spd_bin_totals[i] / sumbins * 100, "%"


## KLUDGE! alias my (new) data structures to BF's original ones...

wfs = dirbins
sumb = sumbins
calm = calmbin
dirless = 0
iw = binheaders

# n = len(wf) # should be 36 or some integral multiple of 360 degrees
n = len(wfs) # should be 36 or some integral multiple of 360 degrees

if max_radius > 0:
    mwf = sumb * max_radius
else:
    mwf = maxsumdir # max(wf)

# rsize = 170 # pt size for maximum "petal" length

sc = float(rsize)/mwf # scale factor, pts per number
da = 360./n
dah = .5*da
a1 = 90+dah # convert meteorology angle to postscript angle
r0 = 10 # an integer, so unit is pts

# now plot using vplot:
# a = vplot.eps_class()

a = vplot.eps_class( fname=ofname, bbx=image_x, bby=image_y )
a.scale(xmin = -1,xmax = 1,ymin = -1,ymax = 1)
a.linewidth(25L) # 25L means .25 pts, a very thin line

for i in range(0,n):
    a2 = a1
    a1 = a2 - da  # notice -da; meteorology is cw, postscript is ccw
    r2 = r0
    for n in range(len(wfs[i])):
        f = wfs[i][n]
        r1 = r2 # last outer radius becomes current inner radius of sector
        r2 = int(r1+sc*f)
        a.color(colors[n])
        a.sector(0., 0., r1, r2, a1, a2, 'F')

a.color(0,0,1)

if max_radius > 0:
    grad = range(1,int(max_radius*100+1))
    for i in grad: # make the blue frequency circles
        if len(grad) > 30:
            if i%5: continue
        elif len(grad) > 10:
            if i%2: continue
        rfreq = i * .01
        lab = "%3.0f%%" % (rfreq * 100)
        r = rfreq * sumb * sc + r0
        a.circle(0.,0.,int(r))
        a.text( int(a.ix(0.) - 24), int(a.jy(0.) - r - 10), 0, 12, lab)
else:
    for i in range(1,6): # make the blue frequency circles
        rfreq = i * .02
        lab = "%3.0f%%" % (rfreq * 100)
        r = rfreq * sumb * sc + r0
        a.circle(0.,0.,int(r))
        a.text( int(a.ix(0.) - 24), int(a.jy(0.) - r - 10), 0, 12, lab)

a.text(.05j, .90j, 0., 18, 'percent frequency per ' + str(da) + ' degrees')

a.color(0,0,0)
a.text(.25j, .95j, 0., 24, title)

# calmf = float(calm)/sumb * 100
# calmfs = "   calm = %f" % calmf 
# a.text(.1j, .05j, 0., 12, calmfs)
# 
# # dirlessf = float(dirless)/sumb * 100
# # dirlessfs = "dirless = %f" % dirlessf 
# # a.text(.1j, .02j, 0., 12, dirlessfs)
# 
# x = .5j
# y = .02j
# dx = .07j
# dy = .07j
# for n in range(len(iw)):
#     a.color(colors[n])
#     a.rect(x+n*dx, y, x+(n+1)*dx, y+dy, 'F')
#     a.color()
#     a.text(x+(n+1)*dx-.02j, y+dy+.01j, 0., 12, str(iw[n]))
# 
# a.close()

x = 0.05j # start speed key at 10% of x dimension
dx = 0.9j / (len(binheaders) + 1)
py = 0.1j
y = 0.03j
dy = 0.03j
py = 0.01j

a.color( (0,0,0) )
a.rect(x, y, x+dx, y+dy)
a.color()

calmf = float(calm)/sumb * 100
calmfs = "%.1f%%" % calmf 
a.text(x+py, y-dy, 0., 14, calmfs)
a.text(x+py, y+dy+py, 0., 14, "calm")

for i in range(0, len(binheaders) ):
    x += dx
    a.color(colors[i])
    a.rect(x, y, x+dx, y+dy, 'F')
    a.color()
    binf = float(spd_bin_totals[i])/sumb * 100
    if binf > 10: 
        binfs = "%.1f%%" % binf 
    else:
        binfs = "%.2f%%" % binf 
    a.text(x+py, y-dy, 0., 14, binfs)
    a.text(x+py, y+dy+py, 0., 14, binheaders[i] )

a.close()