#!/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()