#!/usr/bin/env python
import gipsy
import csv


#-----------------------------------------------------------------------
"""
Changelog.
09-09-2010:
-ATTRIBn= changed by ATTRIB_Xn where X is one of
the object abbreviations 'E', 'R' or 'N'.
-GRAT2= replaced by GRATOVER, added examples in dc1 document
-R_OPT= replaced by IOOPT=
-Attributes are specified as key:value
The order is not important.
-Added keywords for border graticule (so that one always has a
 decent border even if the projection was oblique
-Added keywords for title and x/y labels
-Changed ellipse etc. parameters in a file.
The syntax now is:
   string float float float
String is a position. e.g. "ga 20 ga 30" or "{eq,fk4,B1960} 101.23 {} -20.23"
or 20h10m23.6s -10d2m30s"
"""
#-----------------------------------------------------------------------

# next packages/modules are not part of the GIPSY distribution
# it is not sure that a user has all the necessary functionality available.
try:
   import numpy
except:
   raise Exception, "Cannot find package numpy!"
try:
   #from matplotlib import use
   #use('Qt4Agg')
   #use('GTKCairo')
   from matplotlib.pyplot import figure, setp, savefig, show
   from matplotlib.patches import Polygon
except:
   raise Exception, "Cannot find matplotlib!"
try:
   import pyfits
except:
   raise Exception, "Cannot find module pyfits for reading FITS files!"
try:
   from kapteyn import wcs, maputils, positions
   from kapteyn.mplutil import TimeCallback, gipsy_connect
   from kapteyn.tabarray import readColumns
   from kapteyn  import tabarray
   from kapteyn.positions import str2pos
except:
   raise Exception, "Cannot find the kapteyn package!"


# Global constant for maximum number of positions 
MAXPOS = 20000
MAXATTS = 10

epsilon = 0.0000000001
almostdec0 = 90.0 - epsilon   # Avoid plotting at the wrong side

linekwargs = ['alpha', 'antialiased', 'color', 'linestyle', 'ls',
              'linewidth', 'lw', 'marker', 'markeredgecolor', 'mec',
              'markeredgewidth', 'mew', 'markerfacecolor', 'mfc',
              'markersize', 'ms', 'visible', 'zorder']
polykwargs = ['alpha', 'antialiased', 'color', 'edgecolor', 'ec',
              'facecolor', 'fc', 'fill', 'linestyle', 'ls',
              'linewidth', 'lw', 'visible', 'zorder']
textkwargs = ['alpha', 'backgroundcolor', 'color', 'family',
              'horizontalalignment', 'ha', 'linespacing',
              'rotation', 'size', 'fontsize', 'style',
              'verticalalignment', 'va', 'visible', 'weight', 'zorder']
titlekwargs = textkwargs[:]; titlekwargs.append('x'); titlekwargs.append('y');
hardcopykwargs = ['papertype', 'orientation']
papertypelist = ['letter', 'legal', 'executive', 'ledger',
                 'a0', 'a1', 'a2','a3', 'a4', 'a5', 'a6', 'a7', 'a8', 'a9', 'a10',
                 'b0', 'b1', 'b2','b3', 'b4', 'b5', 'b6', 'b7', 'b8', 'b9', 'b10']
orientationlist = ['landscape', 'portrait']
floatkwargs = ['linewidth', 'lw', 'ms', 'markersize', 'alpha',
               'mew', 'markeredgewidth', 'rotation']
boolkwargs = ['visible', 'fill']
markersymbols = [ '+' , '*' , ',' , '.' , '1' , '2' , '3' , '4' ,
                  '<' , '>' , 'D' , 'H' , '^' , '_' , 'd' , 'h' ,
                  'o' , 'p' , 's' , 'v' , 'x' , ',' ]
linestyles =   [ '-' , '--' , '-.' , ':' , 'None' ]
colorlist = ['b', 'g', 'r', 'c', 'm', 'y',' k', 'w', '0', '1']
halignlist = ['center', 'right', 'left']
valignlist = ['center', 'top', 'bottom', 'baseline']


def cylrange():
   X = numpy.arange(0,400.0,30.0)
   # Replace last two (dummy) values by two values around 180 degrees   
   X[-1] = 180.0 - epsilon
   X[-2] = 180.0 + epsilon
   return X


def szp_border(mu):  
  xp = -mu * numpy.cos(theta*numpy.pi/180.0)* numpy.sin(phi*numpy.pi/180.0)
  yp =  mu * numpy.cos(theta*numpy.pi/180.0)* numpy.cos(phi*numpy.pi/180.0)
  zp =  mu * numpy.sin(theta*numpy.pi/180.0) + 1.0
  a = numpy.linspace(0.0,360.0,500)
  arad = a*numpy.pi/180.0
  rho = zp - 1.0
  sigma = xp*numpy.sin(arad) - yp*numpy.cos(arad)
  sq = numpy.sqrt(rho*rho+sigma*sigma)
  omega = numpy.arcsin(1/sq)
  psi = numpy.arctan2(sigma,rho)
  thetaxrad = psi - omega
  thetax = thetaxrad * 180.0/numpy.pi + 5
  return a, thetax
  

def sin_border(xi, eta):
  a = numpy.linspace(0,360,500)
  arad = a*numpy.pi/180.0
  thetaxrad = -numpy.arctan(xi*numpy.sin(arad)-eta*numpy.cos(arad))
  thetax = thetaxrad * 180.0/numpy.pi + 0.000001  # Little shift to avoid NaN's at border
  return a, thetax
  

def cube_border(proj):
   # Calculate perimeter of QUAD projection
   #dx = proj.crval[0]; dy = proj.crval[1]
   xlo, y = proj.topixel((-45.0-epsilon, 0.0))
   xhi, y = proj.topixel((315+epsilon, 0.0))
   x, ylo = proj.topixel((180, -45))
   x, yhi = proj.topixel((180, 45))
   x1, y = proj.topixel((45-epsilon, 0.0))
   x, ylolo = proj.topixel((0, -135+epsilon))
   x, yhihi = proj.topixel((0, 135-epsilon))
   perimeter = [(xlo,ylo), (x1,ylo), (x1,ylolo), (xhi,ylolo), (xhi,yhihi),
                (x1,yhihi), (x1,yhi), (xlo,yhi), (xlo,ylo)]
   return perimeter


def updatekwargs(key, mes, defval, obj, kwargs):
#-----------------------------------------------------------------
   """
   This function performs a check on the validity of keywords
   for Matplotlib.
   If 
   """
#-----------------------------------------------------------------
   ok = False
   while not ok:
      message = 'Try again'
      try:
         atts = gipsy.userchar(key, mes, default=1, defval=defval, nmax=MAXATTS)
         newkwargs = processkwargs(obj, atts, kwargs)
         ok = not newkwargs is None
      except AttributeError, (code, message):
         ok is False
      except ValueError, (code, message):
         ok is False
      if not ok:
         gipsy.reject(key, message)
   return newkwargs


def processkwargs(obj, atts, kwargs):
#-----------------------------------------------------------------
   """
   This function performs a check on the validity of keywords
   for Matplotlib.
   If 
   """
#-----------------------------------------------------------------
   for att in atts:
      ok = False
      if att == '?':
         key = '?'
      else:
         try:
            key, val = att.split(':',1)  # i.e. at most two elements
         except:
            raise ValueError(22, "Need key:value")
         
      if obj == 'line':
         if key == '?':
            gipsy.anyout("\nHELP --> Attributes for lines and markers:")
            gipsy.anyout(str(linekwargs))
         else:
            ok = key in linekwargs
      if obj == 'poly':
         if key == '?':
            gipsy.anyout("\nHELP --> Attributes for polygons:")
            gipsy.anyout(str(polykwargs))
         else:
            ok = key in polykwargs
      if obj == 'text':
         if key == '?':
            gipsy.anyout("\nHELP --> Attributes for text:")
            gipsy.anyout(str(textkwargs))
         else:
            ok = key in textkwargs
      if obj == 'hard':
         if key == '?':
            gipsy.anyout("\nHELP --> Attributes for hardcopy:")
            gipsy.anyout(str(hardcopykwargs))
         else:
            ok = key in hardcopykwargs
      if obj == 'title':
         if key == '?':
            gipsy.anyout("\nHELP --> Attributes for title:")
            gipsy.anyout(str(titlekwargs))
         else:
            ok = key in titlekwargs
            
      if not ok:
         if key != '?':
            mes = "Key in %s unknown to object %s" % (att, obj)
            raise AttributeError(23, mes)
         else:
            return None

      ok = False
      if key == 'alpha':
         if val == '?':
            gipsy.anyout("\nHELP --> Transparency:")
            gipsy.anyout("alpha must be a float in range 0 and 1\n")
         else:
            try:
               val = float(val)
            except:
               raise ValueError(25, "Argument [%s] must be float"%val)
            ok = (0.0 <= val <= 1.0)
            gipsy.anyout(str(ok))
      elif key in ['color', 'markeredgecolor', 'mec', 'markerfacecolor', 'mfc',
                 'edgecolor', 'ec', 'facecolor', 'fc', 'backgroundcolor']:
         if val == '?':
            gipsy.anyout("\nHELP --> Colors:")
            gipsy.anyout("b(lue), g(reen), r(ed), c(yan), m(agenta), y(ellow), k(black), w(hite)")
            gipsy.anyout("or a string with hex representation e.g. #AA0057")
            gipsy.anyout("or a number between 0 and 1 to set a gray scale, e.g. 0.75\n")
         else:
            ok = val in colorlist or val.startswith('#') or val.startswith('0.')
      elif key in ['ls', 'linestyle']:
         if val == '?':
            gipsy.anyout("\nHELP --> Line styles:")
            gipsy.anyout(str(linestyles))
            gipsy.anyout('')
         else:
            ok = val in linestyles
      elif key == 'marker':
         if val == '?':
            gipsy.anyout("\nHELP --> Marker symbols:")
            gipsy.anyout(str(markersymbols))
            gipsy.anyout('')
         else:
            ok = val in markersymbols
      elif key in ['horizontalalignment', 'ha']:
         if val == '?':
            gipsy.anyout("\nHELP --> Horizontal align:")
            gipsy.anyout(str(halignlist))
            gipsy.anyout('')
         else:
            ok = val in halignlist
      elif key in ['verticalalignment', 'va']:
         if val == '?':
            gipsy.anyout("\nHELP --> Vertical align:")
            gipsy.anyout(str(valignlist))
            gipsy.anyout('')
         else:
            ok = val in valignlist
      elif key == 'orientation':
         if val == '?':
            gipsy.anyout("\nHELP --> Paper orientation:")
            gipsy.anyout(str(orientationlist))
            gipsy.anyout('')
         else:
            ok = val in orientationlist
      elif key == 'papertype':
         if val == '?':
            gipsy.anyout("\nHELP --> Paper type:")
            gipsy.anyout(str(papertypelist))
            gipsy.anyout('')
         else:
            ok = val in papertypelist
      elif key in floatkwargs:
         if val == '?':
            gipsy.anyout("\nHELP --> Floating point parameters")
            gipsy.anyout("These parameters are represented by a floating point number")
            gipsy.anyout("e.g. lw:1.0 or lw:2.0 or lw:0.5\n")
         else:
            try:
               val = float(val)
            except:
               raise ValueError(25, "Argument [%s] must be float"%val)
            ok = True
      elif key in boolkwargs:
         if val == '?':
            gipsy.anyout("\nHELP --> Boolean parameters")
            gipsy.anyout("These parameters are represented by a boolean.")
            gipsy.anyout("The values are either True or False\n")
         else:
            try:
               val = {'true': True, 'false': False}.get(val.lower())
            except:
               raise ValueError(25, "Argument [%s] must be boolean"%val)
            ok = True
      else:    # All the others
         if val == '?':
            gipsy.anyout("No information available!")
         ok = True   # With a risk of crashing the program for wrong values
      if not ok or val == '?':
         if val != '?':
            mes = "Value [%s] in %s not allowed!" % (val, att)
            raise AttributeError(24, mes)
         else:
            return None
      
      kwargs.update({key:val})
   return kwargs


def settattribs(patch, **kwargs):
   if patch is None:
      return
   try:
      setp(patch, **kwargs)
   except AttributeError, message:
      errmes = str(message)
      gipsy.error(1, errmes)


projlist = []   # A global list with projections
class Projectioninfo(object):
#---------------------------------------------------------------------------
   """
   Set defaults for a projection. Append projection to global list
   'projlist'.
   """
#---------------------------------------------------------------------------
   def __init__(self, code, ptype, descr='', **pv):
      self.code = code
      self.descr = descr
      self.wylim = (-90, 90)
      self.wxlim = (0,360)
      self.startx = cylrange()
      self.starty = numpy.arange(-90,100,30.0)
      self.crval1 = 0.0
      if ptype == 'zenithal':
         self.crval2 = almostdec0
         self.startx = numpy.arange(0,360,30)
      else:
         self.crval2 = 0.0
      self.cdelt1 = -5.0
      self.cdelt2 = 5.0
      self.labels_lon = numpy.arange(0,360,30.0)
      self.labels_lat = [-60,-30, 30, 60]
      self.lon_constval = 0.0
      self.lat_constval = 0.0
      self.ptype = ptype
      if ptype in ['cylindrical', 'pseudocylindrical']:
         self.stretched = True
      else:
         self.stretched = False
      self.pv = pv
      projlist.append(self)
      


mu = 2.0; gamma = 20.0
p = Projectioninfo('AZP', 'zenithal', '(azimuthal) perspective projection (AZP)', 
                   PV2_1=mu, PV2_2=gamma)
lowval = (180.0/numpy.pi)*numpy.arcsin(-1.0/mu) + 0.00001
p.starty[0] = lowval

mu = 2.0; phi = 180.0; theta = 60
p = Projectioninfo('SZP', 'zenithal', 'Slant zenithal perspective (SZP)', PV2_1=mu, PV2_2=phi, PV2_3=theta)
p.mu = mu

p = Projectioninfo('TAN', 'zenithal', 'Gnomonic projection (TAN) diverges at theta=0');
p.wylim = (20,90)   # TAN diverges at 0
p.starty = numpy.arange(0,100,30.0); p.starty[0] = p.wylim[0]
p.lat_constval = 19

p = Projectioninfo('STG', 'zenithal', 'Stereographic projection (STG) diverges at theta=-90');
p.wylim = (-60,90)

xi =  -1/numpy.sqrt(6); eta = 1/numpy.sqrt(6)
p = Projectioninfo('SIN', 'zenithal', 'Slant orthograpic projection (SIN) with xi and eta',
                   PV2_1=xi, PV2_2=eta)
p.xi = xi,
p.eta = eta
p.lat_constval = 50.0

p = Projectioninfo('ARC', 'zenithal', 'Zenithal equidistant projection (ARC)');
p.starty[0] = -90.0 + 10*epsilon

p = Projectioninfo('ZPN', 'zenithal', 'Zenithal polynomial projection (ZPN) with PV2_n parameters',
                  PV2_1=1.0);
p.starty[0] = -90.0 + 10*epsilon

p = Projectioninfo('ZEA', 'zenithal', 'Zenith equal area projection (ZEA)')
p.starty[0] = -90.0 + 10*epsilon

p = Projectioninfo('AIR', 'zenithal', 'Airy projection (AIR)', PV2_1=45.0)

#p = Projectioninfo('CYP', 'cylindrical', r"""Gall's stereographic projection (CYP) with
#$\mu = 1$ and $\theta_x = 45^\circ$ (i.e. $\lambda=\frac{1}{2}\sqrt(2)$)""",
#                    PV2_1=1.0, PV2_2=numpy.sqrt(2.0)/2.0)
p = Projectioninfo('CYP', 'cylindrical', r"""Gall's stereographic projection (CYP) with $\mu = 1$ and $\theta_x = 45^\circ$ i.e. $\lambda=\frac{1}{2}\sqrt{2}$""",
                    PV2_1=1.0, PV2_2=numpy.sqrt(2.0)/2.0)

p = Projectioninfo('CEA', 'cylindrical', r"""Lambert's equal area projection (CEA) with $\lambda = 1$""",
                    PV2_1=1.0)

p = Projectioninfo('CAR', 'cylindrical', 'Plate Caree (CAR)');
p = Projectioninfo('MER', 'cylindrical', 'Mercator (MER)');
p.wylim = [-85,85]  # Mercator diverges at +/- 90 deg
p = Projectioninfo('SFL', 'cylindrical', 'Sanson-Flamsteed projection (SFL)')
p = Projectioninfo('PAR', 'cylindrical', 'Parabolic projection (PAR)')
p = Projectioninfo('MOL', 'cylindrical', 'Mollweide\'s projection (MOL).')
p = Projectioninfo('AIT', 'cylindrical', 'Hammer Aitoff');

theta_a = 45
t1 = 20.0; t2 = 70.0
eta = abs(t1-t2)/2.0
p = Projectioninfo('COP', 'conic', r"""Conic perspective projection (COP) with:
$\theta_a=45^\circ$, $\theta_1=20^\circ$ and $\theta_2=70^\circ$""",
                   PV2_1=theta_a, PV2_2=eta)
p.crval2 = theta_a
p.starty = numpy.arange(-30,90,30)

theta_a = -45
t1 = -20.0; t2 = -70.0
eta = abs(t1-t2)/2.0
p = Projectioninfo('COE', 'conic', r"""Conic equal area projection (COE) with:
$\theta_a=-45^\circ$, $\theta_1=-20^\circ$ and $\theta_2=-70^\circ$""",
                   PV2_1=theta_a, PV2_2=eta)
p.crval2 = theta_a
p.starty[-1] = almostdec0

theta_a = 45
t1 = 20.0; t2 = 70.0
eta = abs(t1-t2)/2.0
p = Projectioninfo('COD', 'conic', r"""Conic equidistant projection (COD) with:
$\theta_a=45^\circ$, $\theta_1=20^\circ$ and $\theta_2=70^\circ$""",
                   PV2_1=theta_a, PV2_2=eta)
p.crval2 = theta_a

theta_a = 45
t1 = 20.0; t2 = 70.0
eta = abs(t1-t2)/2.0
p = Projectioninfo('COO', 'conic', r"""Conic orthomorfic projection (COO) with:
$\theta_a=45^\circ$, $\theta_1=20^\circ$ and $\theta_2=70^\circ$""",
                   PV2_1=theta_a, PV2_2=eta)

p.crval2 = theta_a
p.wylim = (-30,90)  # Diverges at -90
p.starty = numpy.arange(-30,90,30)
p.startx = [0,30,60,90,120,150,180-epsilon,180+epsilon,210,240,270,300,330,360-epsilon]

theta1 = 45
p = Projectioninfo('BON', 'pseudoconic', r"""Bonne's equal area projection (BON) with $\theta_1=45^\circ$""",
                   PV2_1=theta1)

p = Projectioninfo('PCO', 'polyconic', 'Polyconic');
p = Projectioninfo('TSC', 'cube', "Tangential spherical cube projection (TSC)")
p = Projectioninfo('CSC', 'cube', "COBE quadrilateralized spherical cube projection (CSC)")
p = Projectioninfo('QSC', 'cube', "Quadrilateralized spherical cube projection (QSC)")


gipsy.init()                                     # Contact Hermes
# ------- Ask for projection --------

def ellipse_handler(cb):
   if cb.key == 'ELLIPSE=':
      annim.Epars = gipsy.userreal(cb.key, nmax=5)
   else:
      pars = annim.Epars
      attrs = gipsy.userchar(cb.key, nmax=4)
      gipsy.anyout(attrs[0])
      
      kwargs = {}
      if len(attrs) > 0:
         kwargs.update({'facecolor':attrs[0]})
      if len(attrs) > 1:
         kwargs.update({'edgecolor':attrs[1]})
      if len(attrs) > 2:
         kwargs.update({'linewidth':attrs[2]})
      if len(attrs) > 3:
         kwargs.update({'alpha':attrs[3]})
      gipsy.anyout(str(kwargs))

      skypol = cb.annim.Skypolygon(prescription="Ellipse", xc=pars[0], yc=pars[1],
                     major=pars[2], minor=pars[3], pa=pars[4], **kwargs)

      skypol.plot(cb.annim.frame)
      cb.annim.frame.figure.canvas.draw()


for i, p in enumerate(projlist):
   s = "%d: %s (%s)" %(i, p.descr, p.ptype)
   gipsy.anyout(s, dev=3)

ok = False
while not ok:
   key = "PROJECTION="
   mes = "Enter number [0,%d] of a valid projection ... [0]: " % (len(projlist)-1)
   projnr = gipsy.userint(key, mes, defval=0, default=1, nmax=1)
   ok = len(projlist) > projnr >= 0
   if not ok:
      gipsy.reject(key, 'You did not enter a number between 0 and %s'%(len(projlist)))
     
projuser = projlist[projnr]
gipsy.anyout("You selected projection: %s"%(projuser.descr))

#------- Ask for pixel size (CDELT) --------
key = "CDELT_XY="
cdeltx = projuser.cdelt1
cdelty = projuser.cdelt2
mes = "Enter pixel size in X and Y in deg. ... [%.2g,%.2g]: "%(cdeltx, cdelty)
cdeltx, cdelty = gipsy.userdble(key, mes, defval=(cdeltx, cdelty), default=1, nmax=2)
   
naxis1 = abs(400.0/cdeltx)
naxis2 = abs(200.0/cdelty)
if not projuser.stretched:
   naxis1 = max(naxis1,naxis2)
   naxis2 = naxis1

key = "NAXIS_XY="
mes = "Enter length of axes in pixels ... [%d %d]: "%(naxis1, naxis2)
naxis1, naxis2 = gipsy.userdble(key, mes, defval=(naxis1, naxis2), default=1, nmax=2)

key = "CRPIX_XY="
crpix1 = naxis1/2.0; crpix2 = naxis2/2.0;
if projuser.ptype == 'cube':
   crpix1 = naxis1*5.0/6.0
mes = "Pixel coordinates for center ... [%d %d]: "%(crpix1, crpix2)
crpix1, crpix2 = gipsy.userdble(key, mes, defval=(crpix1, crpix2), default=1, nmax=2)


#------- Ask for pixel size (CRVAL) --------
key = "CRVAL_XY="
crvalx = projuser.crval1
crvaly = projuser.crval2
mes = "Enter Projection center in X and Y ... [%.2g,%.2g]: "%(crvalx, crvaly)
crvalx, crvaly = gipsy.userdble(key, mes, defval=(crvalx, crvaly), default=1, nmax=2)

#------- Ask for sky system --------
key = "SKY="
defval = 1
mes = "Sky system 1=equat. 2=ecl 3=galactic 4=su.gal. ... [%s]: "%defval
sky = gipsy.userint(key, mes, defval=defval, default=1, nmax=1)
if sky < 1:
   sky = 1
if sky > 4:
   sky = 4
if sky == 1:
   ctypex = 'RA---'
   ctypey = 'DEC--'
if sky == 2:
   ctypex = 'ELON-'
   ctypey = 'ELAT-'
if sky == 3:
   ctypex = 'GLON-'
   ctypey = 'GLAT-'
if sky == 4:
   ctypex = 'SLON-'
   ctypey = 'SLAT-'


if crvalx != 0.0 and projuser.ptype == 'cylindrical':
   projuser.startx += crvalx


header = {'NAXIS'  : 2,
          'NAXIS1' : naxis1,
          'NAXIS2' : naxis2,
          'CTYPE1' : '%s%s'%(ctypex,projuser.code),
          'CRVAL1' : crvalx,
          'CRPIX1' : crpix1,
          'CUNIT1' : 'deg',
          'CDELT1' : cdeltx,
          'CTYPE2' : '%s%s'%(ctypey,projuser.code),
          'CRVAL2' : crvaly,
          'CRPIX2' : crpix2,
          'CUNIT2' : 'deg', 'CDELT2' : cdelty,
         }
# Add PV elemens if there are any.
if len(projuser.pv):
   header.update(projuser.pv)

# Figure size
if projuser.stretched:
   figsize = [12,6]
else:
   figsize = [12,9]
key = "FIGSIZE="
mes = "Figure width, height in inches ... [%g %g]: "% (figsize[0], figsize[1])
figsize = gipsy.userreal(key, mes, default=1, defval=figsize, nmax=2)
if len(figsize) == 1:
   figsize.append(figsize[0])
fig = figure(figsize=figsize)

# Matplotlib Axes object, which we call a frame
key = "FRAME="
x0 = y0 = 0.05
x1 = 0.95; y1 = 0.9
mes = "Frame x0,y0,x1,y1 (norm. coords) ... [%g %g %g %g]: "% (x0,y0,x1,y1)
framesize = gipsy.userreal(key, mes, default=1, defval=(x0,y0,x1,y1), nmax=4)
w = x1 - x0
h = y1 - y0
frame = fig.add_axes((x0,y0,w,h))

key = "MERIDIANS="
defval = projuser.startx[:]
mes = "Meridians at longitudes ... [calculated]: "
startx = gipsy.userreal(key, mes, default=1, defval=defval, nmax=100)

key = "PARALLELS="
defval = projuser.starty[:]
mes = "Parallels at latitudes ... [calculated]: "
starty = gipsy.userreal(key, mes, default=1, defval=defval, nmax=100)

f = maputils.FITSimage(externalheader=header)
annim = f.Annotatedimage(frame)
grat = annim.Graticule(wylim=projuser.wylim, wxlim=projuser.wxlim,
                       startx=startx , starty=starty)
key = "GR1_ATTS="
mes = "Plot atts base grat ... [color:0.75 lw:1]: "
defval = ''
kwargs = {'color':'0.75', 'lw':1}
kwargs = updatekwargs(key, mes, defval, 'line', kwargs)

# Change default label along X-axis
grat.setp_gratline(**kwargs)
key = "XLABEL="
mes = "Text to label X-axis ... [default]: "
defval = ''
lonlabel = gipsy.usertext(key, mes, default=1, defval=defval)

# Ask and set properties for label 
key = "LABX_ATTS="
mes = "Plot atts X-axis label ... [defaults]: "
defval = ''
kwargs = {}
lon_kwargs = updatekwargs(key, mes, defval, 'text',  kwargs)
if lonlabel == '':
   grat.setp_plotaxis("bottom", **lon_kwargs)
else:
   grat.setp_plotaxis("bottom", label=lonlabel, **lon_kwargs)

# Same for Y label
key = "YLABEL="
mes = "Text to label Y-axis ... [default]: "
defval = ''
latlabel = gipsy.usertext(key, mes, default=1, defval=defval)

key = "LABY_ATTS="
mes = "Plot atts Y-axis label ... [defaults]: "
defval = ''
kwargs = {}
lat_kwargs = updatekwargs(key, mes, defval, 'text',  kwargs)

if latlabel == '':
   grat.setp_plotaxis("left", **lat_kwargs)
else:
   grat.setp_plotaxis("left", label=latlabel, **lat_kwargs)

key = "ILABLONS="
mes = "Plot graticule labels at longitudes ... [calculated]: "
defval = startx
labels_lon = gipsy.userreal(key, mes, default=1, defval=defval, nmax=100)

# Ask plot properties for labels inside the plot.
key = "ILABLON_ATTS="
mes = "Plot atts lon. labels ... [color:r  fontsize:10]: "
defval = ''
kwargs = {'color':'r', 'fontsize':10}
lon_kwargs = updatekwargs(key, mes, defval, 'text',  kwargs)

key = "ILABLATS="
defval = []
for y in starty:
   if not y in [-90.0, 0.0, 90.0]:
      defval.append(y)
mes = "Plot graticule labels at latitudes ... [calculated]: "
labels_lat = gipsy.userreal(key, mes, default=1, defval=defval, nmax=100)

key = "ILABLAT_ATTS="
mes = "Plot atts lat. labels ... [color:b  fontsize:10]: "
defval = ''
kwargs = {'color':'b', 'fontsize':10}
atts = gipsy.userchar(key, mes, default=1, defval=defval, nmax=MAXATTS)
lat_kwargs = updatekwargs(key, mes, defval, 'text',  kwargs)

# Prevent plotting tick labels along the plot axes when the plot is
# not completeley contained in the frame and therefore crosses the 
# plotlines, which by default triggers plotting tick labels
grat.setp_ticklabel(visible=False)

# Format plot labels inside the plot
key = "LABFORMAT="
mes = "Enter format for longitude labels ... [Dms]: "
defval = "Dms"
fmt = gipsy.userchar(key, mes, default=1, defval=defval, nmax=1)
il1 = grat.Insidelabels(wcsaxis=0,
                  world=labels_lon, constval=projuser.lat_constval, fmt=fmt)
il1.setp_label(**lon_kwargs)
# The latitudes are plotted at the first longitude entered
# in MERIDIANS= (i.e. in startx)
il2 = grat.Insidelabels(wcsaxis=1, 
                  world=labels_lat, constval=startx[0], fmt='Dms')
il2.setp_label(**lat_kwargs)

# Overlay more graticules
grat2 = None
key = "GRATOVER="
mes = "Enter sky system for overlay 2nd graticule ... [skip]: "
defval = ''
sky2 = gipsy.userchar(key, mes, default=1, defval=defval)
if sky2 != '':
   # User wants a second graticule
   grat2 = annim.Graticule(skyout=sky2, wylim=projuser.wylim, wxlim=projuser.wxlim,
                           startx=projuser.startx , starty=projuser.starty)
   grat2.setp_axislabel(plotaxis=("left","bottom"), visible=False)
   key = "GR2_ATTS="
   mes = "Plot atts 2nd grat ... [color:'#aa2222' lw:0.8 alpha:0.7]: "
   defval = ''
   kwargs = {'color':'#aa2222', 'lw':0.8, 'alpha':0.7}
   kwargs = updatekwargs(key, mes, defval, 'line', kwargs)
   grat2.setp_gratline(**kwargs)
   
   # Line at latitude 0, e.g. for galactic plane:
   key = "LAT0_ATTS="
   mes = "Plot atts for lat. 0  ... [color:r lw:1 alpha:1]: "
   defval = ''
   kwargs = {'color':'r', 'lw':1.0, 'alpha':1.0}
   kwargs = updatekwargs(key, mes, defval, 'line', kwargs)
   grat2.setp_gratline(wcsaxis=1, position=0.0, **kwargs)

# We need a non oblique border at all times:
header_border = header.copy()
header_border['CRVAL1'] = 0.0

gratborder = None
lineborder = None
cubeperimeter = None
if projuser.ptype in ['cylindrical', 'polyconic', 'pseudoconic']:
   header_border['CRVAL2'] = 0.0
   gratborder = annim.Graticule(header=header_border, axnum=(1,2), wylim=projuser.wylim, wxlim=projuser.wxlim,
                             startx=(180-epsilon, 180+epsilon), skipy=True)
   gratborder.setp_axislabel(plotaxis=("left","bottom"), visible=False)
   lineborder = None
elif projuser.code == 'SZP':
   phi, thetax = szp_border(projuser.mu)
   lineborder = grat.addgratline(phi, thetax, pixels=False)
elif projuser.code == 'SIN':
   phi, thetax = sin_border(projuser.xi, projuser.eta)
   lineborder = grat.addgratline(phi, thetax, pixels=False)
elif projuser.ptype == 'cube':
   header_border['CRVAL2'] = 0.0
   projcube = wcs.Projection(header_border)
   lineborder = cube_border(projcube)  # This is not a 'gratline' but vertices for a polygon

  
# Plot attributes for border graticule
key = "BOR_ATTS="
mes = "Plot atts. of base grat ... [color:k lw:1]: "
defval = ''
kwargs = {'color':'k', 'lw':1}
kwargs = updatekwargs(key, mes, defval, 'line', kwargs)
if not gratborder is None:
   gratborder.setp_lineswcs0(**kwargs)
   gratborder.setp_lineswcs1(**kwargs)
   # Prevent plotting tick labels along the plot axes when the plot is
   # not completeley contained in the frame and therefore crosses the 
   # plotlines, which by default triggers plotting tick labels
   gratborder.setp_ticklabel(visible=False)
if not lineborder is None:
   if projuser.ptype == 'cube':
      Xp, Yp = zip(*lineborder)
      frame.plot(Xp, Yp, **kwargs)
   else:
      grat.setp_linespecial(lineborder, **kwargs)




# File name and parameters for hardcopy on disk.
key = "HARDCOPY="
mes = "Name of hardcopy of plot on disk ... [skip]: "
defval = ''
hardcopy = gipsy.userchar(key, mes, default=1, defval=defval, nmax=1)
if hardcopy != '':
   key = "HC_ATTS="
   mes = "Atts hardcopy ... [papertype:a4 orientation:landscape]: "
   defval = ''
   kwargs = {'papertype':'a4', 'orientation':'landscape'}
   hc_kwargs = updatekwargs(key, mes, defval, 'hard', kwargs)

# Set title for Matplotlib
key = "TITLE="
defval = projuser.descr
mes = "Title above plot ... [%s]: " % defval
title  = gipsy.usertext(key, mes, default=1, defval=defval)
#title = r"""%s"""%title
key = "TITLE_ATTS="
mes = "Atts title ... [color:g y:1.02]: "
defval = ''
kwargs = {'color':'g', 'y':1.02}
kwargs = updatekwargs(key, mes, defval, 'title', kwargs)
frame.set_title(title, **kwargs)

# Start loop asking keywords to specify the required shapes and attributes
basekey = "SHAPE"
parkey  = "PARAMS"
attrkey = "ATTRIBS"
legname = "LEGEND"
i = 0
cont = True
# Start loop asking keywords to specify the required shapes and attributes
# for objects to be plotted onjto an all sky figure.
legendpatches = []
legendnames = []
while cont:
   i += 1
   # Select object (shape, marker or image)
   key = "SHAPE"+"%d="%i
   mes = "E(ell),R(ect),N(-poly),M(arker),I(mage) ... [start plot]: "
   defval = ''
   ts = gipsy.usercharu(key, mes, default=1, defval=defval)
   cont = ts != defval
   if cont:
      sp = None         # The last sky polygon
      if ts in ['E','R', 'N', 'P']:
         # Get the attributes
         key = attrkey+"_%s%d=" % (ts, i)
         mes = "Enter plot attributes ...... [fc:r ec:g lw:1 alpha:0.6]: "
         defval = ''
         kwargs = {'fc':'r', 'ec':'g', 'lw':1.0, 'alpha':0.6}
         kwargs = updatekwargs(key, mes, defval, 'poly', kwargs)
         # Ask origin of the shape parameters, file or manual input
         key = "IOOPT"+"%d="%i
         defval = 'M'
         mes = "Parameters from F(ile) or M(anual) input ... [%s]: "%defval
         org = gipsy.usercharu(key, mes, default=1, defval=defval, nmax=1)
         if org == 'F':
            key = "FILENAME"+"%d="%i
            mes = "Enter name of file:"
            fn = gipsy.userchar(key, mes, nmax=1)
            fp = open(fn, "rb")
            reader = csv.reader(fp, delimiter=' ')
            for row in reader:
               xcyc,major,minor,pa = row
               major = float(major); minor= float(minor); pa = float(pa)
               nangles = float(minor)
               world, pixels, units, errmes = positions.str2pos(xcyc, annim.projection)
               if errmes != '':
                  gipsy.error(4, errmes)
               xc = world.T[0]; yc = world.T[1]
               sp = annim.Skypolygon(prescription=ts, xc=xc, yc=yc, major=major,
                                     minor=minor, nangles=nangles, pa=pa, **kwargs)
               settattribs(sp.p1, **kwargs)
               settattribs(sp.p2, **kwargs)
            fp.close()  # Close the parameters file
         elif org == 'M':
            # First manual set of parameters is the position of the origin
            key = "CPOS"+"%d="%i
            mes = "Central position xc, yc: "
            badinput = True
            while badinput:
               pos = gipsy.usertext(key, mes)
               world, pixels, units, errmes = positions.str2pos(pos, annim.projection)
               if errmes != '':
                  gipsy.reject(key, errmes)
               else:
                  badinput = False
            xc = world.T[0]; yc = world.T[1]
            # And ask the other parameters. Adjust message for each kind of shape
            if ts == 'E':
               mes = "Enter maj, min, pa: "
            elif ts == 'R':
               mes = "Enter height, width, pa: "
            elif ts == 'N':
               mes = "Enter radius, nangles, pa: "
            key = parkey+"%d="%i
            pars = gipsy.userreal(key, mes, default=4, nmax=3)
            major, minor, pa = pars
            nangles = minor
            sp = annim.Skypolygon(prescription=ts, xc=xc, yc=yc, major=major,
                                  minor=minor, nangles=nangles, pa=pa)
            settattribs(sp.p1, **kwargs)
            settattribs(sp.p2, **kwargs)

      # Marker
      if ts == 'M':
         key = "M_ORG"+"%d="%i
         defval = 'P'
         mes = "Markers as L(ong,lat) or as P(ositions) ... [%s] "%defval
         org = gipsy.usercharu(key, mes, default=1, defval=defval, nmax=1)
         if org == 'L':
            key = "LONS"+"%d="%i
            mes = "Enter longitude(s) in deg.: "
            lons = gipsy.userreal(key, mes, nmax=MAXPOS);
            key = "LATS"+"%d="%i
            mes = "Enter %d latitudes (deg.): " % len(lons)
            lats = gipsy.userreal(key, mes, nmax=len(lons));
         elif org == 'P':
            key = "POSITIONS"+"%d="%i
            mes = "Enter positions: "
            badinput = True
            while badinput:
               pos = gipsy.usertext(key, mes, nmax=4096)
               world, pixels, units, errmes = positions.str2pos(pos, annim.projection)
               if errmes != '':
                  gipsy.reject(key, errmes)
               else:
                  badinput = False
            lons = world.T[0]; lats = world.T[1]
         key = "M_OPT"+"%d="%i
         defval = 'M'
         mes = "Plot as M(arkers), I(rregular polygon) or B(oth) ... [%s]: "%defval
         opt = gipsy.usercharu(key, mes, default=1, defval=defval, nmax=1)
         if opt in ['M', 'B']:
            key = attrkey+"_M%d="%i
            mes = "Enter plot attributes ... [color:r marker:o ms:2 mew:1 alpha:0.6]: "
            kwargs = {'color':'r', 'marker':'o', 'ms':1.0, 'mew':1.0, 'alpha':0.6}
            defval = ''
            kwargs = updatekwargs(key, mes, defval, 'line', kwargs)
            sp = annim.Marker(x=lons, y=lats, mode='world', **kwargs)
            #legendpatches.append(sp)
            #legendnames.append("Bla bla")
         if opt in ['I', 'B']:
            key = attrkey+"_I%d="%i
            mes = "Enter plot attributes ...... [fc:r ec:g lw:1 alpha:0.6]: "
            defval = ''
            kwargs = {'facecolor':'r', 'edgecolor':'g', 'linewidth':1.0, 'alpha':0.6}
            kwargs = updatekwargs(key, mes, defval, 'poly', kwargs)
            sp = annim.Skypolygon(prescription=None, lons=lons, lats=lats)
            settattribs(sp.p1, **kwargs)
            settattribs(sp.p2, **kwargs)
      
      if not sp is None:
         key = legname+"%d="%i
         mes = "Entry for legend ... [skip]: "
         s = gipsy.usertext(key, mes, default=1, defval='')
         if s != '':
            legendpatches.append(sp)
            legendnames.append(r'%s'%s)

annim.plot()
# Do text labels after plotting maputils objects. Then gratborder's frame
# is known so that labels can be plotted on top of the grid lines.

lastframe = grat.frame
if not grat2 is None:
  lastframe = grat.frame
if not gratborder is None:
  lastframe = gratborder.frame

cont = True
i = 0
# Start loop asking keywords to specify the text labels
while cont:
   i += 1
   # Select object (shape, marker or image)
   key = "ANNTXT"+"%d="%i
   mes = "Enter text to plot in figure ... [stop loop]: "
   defval = ''
   labtxt = gipsy.usertext(key, mes, default=1, defval=defval)
   cont = labtxt != defval
   if cont:
      # First manual set of parameters is the position of the origin
      key = "ANNPOS"+"%d="%i
      mes = "Central position xc, yc: "
      badinput = True
      while badinput:
         pos = gipsy.usertext(key, mes)
         world, pixels, units, errmes = positions.str2pos(pos, annim.projection)
         if errmes != '':
            gipsy.reject(key, errmes)
         else:
            badinput = False
      xc = pixels.T[0]; yc = pixels.T[1]
      key = "ANNATTS"+"%d="%i
      mes = "Enter annotation attributes ... [color:r fontsize:10]: "
      kwargs = {'color':'r', 'fontsize':10}
      defval = ''
      kwargs = updatekwargs(key, mes, defval, 'text', kwargs)
      lastframe.text(xc, yc, labtxt, **kwargs)
      


# Do legend after plotting maputils objects. Then the patches are known

if len(legendpatches) > 0:
   legendobjects = []
   key = "LEGENDLOC="
   mes = "Enter a legend location (number 0-10) ... [1]: "
   defval = 1
   loc = gipsy.userint(key, mes, default=1, defval=defval, nmax=1)
   for l,n in zip(legendpatches, legendnames):
      legendobjects.append(l.patch)
      try:
         print len(l.patch)
      except:
         pass
   
   leg = lastframe.legend(legendobjects, legendnames, shadow=True, loc=loc)
   # matplotlib.text.Text instances
   for t in leg.get_texts():
      t.set_fontsize(8)    # the legend text fontsize

# Plot alternative borders
if cubeperimeter != None:
   #p = Polygon(cubeperimeter, facecolor='#d6eaef', lw=2)
   #frame.add_patch(p)
   Xp, Yp = zip(*cubeperimeter)
   frame.plot(Xp, Yp, color='r')


if hardcopy != '':
   savefig(hardcopy, **hc_kwargs)
annim.interact_toolbarinfo()
annim.interact_writepos(gipsy=True, wcsfmt="%f",zfmt=None, pixfmt=None, hmsdms=False)

gipsy_connect()

#gipsy.KeyCallback(ellipse_handler, 'E_ATTR=', annim=annim)
#gipsy.KeyCallback(ellipse_handler, 'ELLIPSE=', annim=annim)
#gipsy.xeq('!skyshapes.py', "DUMMY=")

show()
gipsy.finis()