#!/usr/bin/env python #---------------------------------------------------------------------- # FILE: wcsflux.py (GIPSY version) # PURPOSE: Calculate sum of image values and area defined by a polygon. # AUTHOR: M.G.R. Vogelaar, University of Groningen, The Netherlands # DATE: November 26, 2009 # UPDATE: November 26, 2009 # VERSION: 1.0 # # (C) University of Groningen # Kapteyn Astronomical Institute # Groningen, The Netherlands # E: gipsy@astro.rug.nl #---------------------------------------------------------------------- import gipsy import os.path # to check existence of a file import ioutils import datetime from kapteyn.maputils import colormaps as cmlist from kapteyn.positions import minmatch, unitfactor from kapteyn.maputils import FITSimage from kapteyn.mplutil import gipsy_connect from kapteyn import shapes from matplotlib.pyplot import figure, show, draw from matplotlib.figure import SubplotParams from numpy.ma import masked_invalid from math import sqrt from types import StringType from re import split as re_split __version__ = '1.1' #TODO: #-gebruiker toestaan om frames met de hand op te geven. #-Super title def parstr2dict(s): """ Convert string with parameters in key:value format to dictionary """ kwargs = {} if s != '': # Result of usertext is a list ss = re_split('[ ,\t]+', s) for sss in ss: if len(ss) > 0: try: k, v = sss.split(':') try: kwargs[k] = float(v) except: kwargs[k] = v except: pass return kwargs def plmarker(cb): """ #--------------------------------------------------------------------- This function is called with a frequency set by a timer function. It reads keyword MARKER= which must be specified while we are in Matplotlib's main loop. Positions are converted to pixel coordinates and these are plotted in all images on the canvas. The sets are stored in the global list 'allsets'. #--------------------------------------------------------------------- """ key = cb.key pos = gipsy.usertext(key, "", defval=None, default=2) if pos != None: for gipsyset in cb.allsets: grids, world, pix, units, errmes = gipsyset.wcspos(0, userpos=pos) if not len(grids): if errmes: gipsy.anyout(errmes) gipsy.cancel(key) gipsy.cancel(key) return else: for g,w in zip(grids, world): s = str(g) + "=" + str(w) + ' ' + str(units) gipsy.anyout(s) xy = pix.T x = xy[0]; y = xy[1] # Plot in all subsets for annim in gipsyset.allim: annim.frame.plot(x, y, '+', markersize=10, markeredgewidth=2, color='y') draw() gipsy.cancel(key) gipsy.cancel(key) # Repeat to remove last one def get_graticuleproperties(gipsyset, inpnr=None, deflev=1): #---------------------------------------------------------------------- """ Deal with properties of axes and graticule lines. """ #---------------------------------------------------------------------- if inpnr == None: pp = '=' else: pp = "%d="%inpnr spectral = longitude = latitude = False gipsyset.gratstart = [None,None] gipsyset.gratdelta = [None,None] gipsyset.gratfmt = [None,None] gipsyset.gratfun = [None,None] gipsyset.gratunits = [None,None] # Note for programmers: The 'gipsyset' object has a number # of attributes that are added in several stages, For instance # it has a projection object and from that object a sub-projection is derived. # Both have the attributes of the Projection class but in the axes order # as specified when you entered the image axes. # One can prove this with command: # print "Units", gipsyset.subproj.cunit, gipsyset.projection.cunit opts = ['X', 'Y'] for ax in [0, 1]: keystart = 'AXSTART' + opts[ax] + pp #gipsy.anyout("Axis number: %d"%gipsyset.projaxnum[ax]) mes = "Enter value of start label(s) on %s axis ... [calculated]:"%(opts[ax]) gipsyset.gratstart[ax] = gipsy.usertext(keystart, mes, defval=None, default=deflev) keydelta = 'AXDELTA' + opts[ax] + pp mes = "Enter value of step on %s axis ... [calculated]:"%(opts[ax]) gipsyset.gratdelta[ax] = gipsy.usertext(keydelta, mes, defval=None, default=deflev) # Ask for a spectral translation origunit = gipsyset.subproj.cunit[ax] axtype = gipsyset.subproj.ctype[ax].split('-')[0] if origunit != '': # Empty units (e.g. Stokes axis) cannot be converted. key = "AXUNITS" + opts[ax] + pp mes = "Enter units for %s (i.e. %s) axis ... [%s]:"%(opts[ax], axtype, origunit) gipsyset.gratunits[ax] = gipsy.usertext(key, mes, defval=None, default=deflev) return def imagespec(inpnr, taskname, deflev): """------------------------------------------------------------------------------- Purpose: This function is written to gather all relevant user input for one data set (either a FITS file or a GIPSY set). Input: inpnr- The number of the current data set. This number is appended to all relevant keywords. It should start with 1. Returns: A list with images (if more than one slice/subset was entered. Notes: Examples: -------------------------------------------------------------------------------""" pp = "%d="%inpnr # Create the string that needs to be appended to all keywords dim = 2 # Allowed dimension for images is 2 key = "DATASET" + pp # Append keyword with index number mes = 'Enter name of FITS file or GIPSY set:' hdukey = "HDUNUM" + pp altkey = "ALTHEAD" + pp axkey = "AXNAMES" + pp swapkey = "SWAPAXES" + pp if inpnr == 1: listopt = True else: listopt = False gipsyset = ioutils.wcsinp(key, mes, hdukey=hdukey, altkey=altkey, axkey=axkey, axsuffix=pp, swapkey=swapkey, requireddim=2, listopt=listopt) # Note that gipsyset has two projection objects called 'projection' and 'subproj' # Projection objects have an attribute 'gridmode' which is set to True by default # in a call to wcsinp() if gipsyset == None: # Nothing to do. Return 0 as the number of subsets return None, 0 gipsyset.printinfo() subdim = gipsyset.ndims(subset=True) # Set limits to the axes key = "BOX" + pp ublo, ubhi = gipsyset.wcsbox(key=key, mes="", default=1, showdev=0, option=0) # These are in order of (sub) projection object 'subproj' gipsyset.pxlim = gipsyset.axlimits[0] gipsyset.pylim = gipsyset.axlimits[1] # Ask for a spectral translation key = "SPECTRANS" + pp gipsyset.subspectraltrans, un = gipsyset.set_spectraltranslation(key, mes=None, inside=True, default=deflev) # At this stage we have either None or a spectral translation for one of the # subset axes! If there is a spectral axis in the subset, then ask the # units for this axis gipsyset.repspectraltrans = None gipsyset.repspectralunits = None gipsyset.orispectralunits = None if gipsyset.subspectraltrans is None: # There is no subset axis that is spectral. Perhaps one of the repeat # axes is spectral. gipsyset.repspectraltrans = None if len(gipsyset.subsets) > 1: # Use the same SPECTRANS= keyword mes = "Spectral translation for subset labeling (or ?) ... [native]:" # Now check the repeat axes with inside=False gipsyset.repspectraltrans, un = gipsyset.set_spectraltranslation(key, mes, inside=False, default=deflev) if not gipsyset.repspectraltrans is None: key = "SUNITS" + pp mes = "Units for spectral slice info ... [%s]:"%un gipsyset.repspectralunits = gipsy.userchar(key, mes, defval=None, default=deflev) # For spatial maps (images with two spatial axes, one longitude # and one latitude) one can change the celestial system for the graticules, # the WCS labels and the cursor information. gipsyset.u_skyout = None # Skyout set by user if gipsyset.spatialmap: key = "SKYOUT" + pp gipsyset.u_skyout = gipsyset.set_skyout(key=key, default=deflev) # Display a list with color map names and prompt the user # to set a color map #cmlist = cmapnames #cmlist.sort() key = "COLORMAP" + pp ok = False while not ok: st = gipsy.userchar(key, 'Enter name of color map or ? for list: ... [jet]:', defval='jet', default=deflev) if st == '?': gipsy.anyout("=========== Available color maps ===========") n = 0 maxcol = 6 while n < len(cmlist): col = 0 cmstr = "" while col < maxcol: if n < len(cmlist): if col > 0: cmstr += ", " + cmlist[n] else: cmstr += cmlist[n] n += 1 col += 1 gipsy.anyout(cmstr) ok = False gipsy.wkey(key) # First make empty gipsy.cancel(key) # Then cancel else: colormap = st if st != "": indx = minmatch(st, cmlist, case=1) if indx == None: if os.path.isfile(st): # Perhaps it is a lut on disk ok = True else: gipsy.reject(key, "Could not find a matching color map!" ) ok = False elif indx == -1: gipsy.reject(key, "Ambiguous input!") ok = False else: colormap = cmlist[indx] ok = True gipsy.wkey(key+colormap) else: gipsy.reject(key, "Nothing entered!") gipsy.anyout(taskname + ": Selected colormap: %s"%colormap) gipsyset.colormap = colormap # Ask user a color scale scales = {'1': 'linear', '2': 'log', '3': 'exp', '4': 'sqrt', '5': 'square'} # User can select a color scale key = "COLORSCALE" + pp mes = "1=linear 2=log 3=exp 4=sqrt 5=square ... [1]:" scale = gipsy.userint(key, mes, defval=1, default=deflev, nmax=1) if scale < 1: scale = 1 if scale > 5: scale = 5 gipsyset.colorscale = scales[str(scale)] # Must be one of the strings key = "CSLOPE" + pp mes = "Slope image values vs colours ... [1.0]:" gipsyset.cslope = gipsy.userdble(key, mes, defval=1.0, default=deflev, nmax=1) key = "COFFSET" + pp mes = "Offset image values vs colours ... [0.0]:" gipsyset.coffset = gipsy.userdble(key, mes, defval=0.0, default=deflev, nmax=1) # For the current subsets, find the min and max to set a # default for the color scaling. Use a masked array to avoid # counting blanks (which are covered by numpy's 'isinf()' method) mini = maxi = None for cword in gipsyset.subsets: z = gipsyset.subimage(cword) subdata = masked_invalid(z) dmin = subdata.min() dmax = subdata.max() if mini == None: mini = dmin else: if dmin < mini: mini = dmin if maxi == None: maxi = dmax else: if dmax > maxi: maxi = dmax # Get the min and max from the user, using default calculated previously key = "VMINMAX" + pp dmin = mini dmax = maxi mes = "Scale colors between image values ... [%.2g,%.2g]:"%(dmin, dmax) vminmax = gipsy.userdble(key, mes, defval=(dmin,dmax), default=deflev, nmax=2) if (len(vminmax) == 1): vminmax.append(dmax) gipsy.anyout(taskname+": Colors scale between %g and %g"%(vminmax[0], vminmax[1])) gipsyset.vminmax = vminmax # Process contours key = "WANTCONT" + pp mes = "Add contours? ... Y/[N]:" wantcont = gipsy.userlog(key, mes, default=deflev, defval=False, nmax=1) maxlevels = 100 levels = None colors = None if wantcont: key = "CLEVELS" + pp mes = "Enter contour levels between %g and %g ... [automatic]:"%(vminmax[0], vminmax[1]) levels = gipsy.userdble(key, mes, defval=None, default=deflev, nmax=maxlevels) gipsyset.levels = levels gipsy.anyout("b=blue, g=green, r=red, c=cyan, m=magenta, y=yellow, k=black, w=white") key = "CCOLORS" + pp mes = "Enter one or more colors for contours ... [r]:" colors = gipsy.userchar(key, mes, defval='r', default=deflev, nmax=maxlevels) gipsyset.colors = colors gipsyset.wantcont = wantcont key = "WANTCOLORBAR" + pp mes = "Add a colorbar? ... Y/[N]:" wantcolorbar = gipsy.userlog(key, mes, default=1, defval=False, nmax=1) gipsyset.wantcolorbar = wantcolorbar if wantcolorbar: key = "COLBARFRAME" + pp mes = "Enter x0,y0, width, height ... [calculated]:" colbarframe = gipsy.userreal(key, mes, default=deflev+4, defval=None, nmax=4) gipsyset.colbarframe = colbarframe gipsyset.cbframemode = '' if not colbarframe is None: key = "CBFRAMEMODE" + pp mes = "Frame coords: F=Norm. in frame, D=Data coords ... [Norm. in figure]: " s = gipsy.userchar(key, mes, defval='', default=deflev) gipsyset.cbframemode = s gipsyset.cborientation = 'V' key = "CBORIENT" + pp mes = "Horizontal color bar (H) or vertical (V)? ... H/[V]:" cborientation = gipsy.userchar(key, mes, default=deflev, defval='V', nmax=1) cborientation = cborientation.upper() if cborientation not in ['V', 'H']: cborientation = 'V' gipsyset.cborientation = cborientation if gipsyset.cborientation == 'V': gipsyset.cborientation = 'vertical' else: gipsyset.cborientation = 'horizontal' # One can label the color bar with the units of the data if wantcolorbar: if gipsyset.header.has_key('BUNIT'): cbunits = gipsyset.header['BUNIT'] else: cbunits = '' key = "CBLABEL" + pp mes = "Label for color bar ... [%s]:"%cbunits gipsyset.cbunits = gipsy.userchar(key, mes, default=deflev, defval=cbunits, nmax=1) # Does user wants to draw graticules and wcs labels? key = "GRATICULES" + pp mes = "Draw graticule lines and WCS labels? ... [Y]/N:" drawgraticules = gipsy.userlog(keyword=key, message=mes, defval=True, default=1) gipsyset.drawgraticules = drawgraticules if gipsyset.drawgraticules: get_graticuleproperties(gipsyset, inpnr, deflev=deflev) key = "GRIDLAB" + pp mes = "Do you want to plot grid labels? ... Y/[N]:" gridlabels = False gipsyset.gridlabels = gipsy.userlog(key, mes, default=deflev, defval=gridlabels) if gipsyset.gridlabels: key = "GRIDLABATT" + pp mes = "Plot attributes for plot object ... [color:g]:" defval = "color:g" s = gipsy.usertext(key, mes, default=1, defval=defval) gipsyset.gridlabelkwargs = parstr2dict(s) key = "BEAMPOS" + pp mes = "Enter position of beam ... [no beam]:" beampos = None gipsyset.beampos = gipsy.usertext(key, mes, default=deflev, defval=beampos) if not gipsyset.beampos is None: key = "FWHMUNITS" + pp mes = "Units of FWHM major&minor axes ... [arcmin]:" gipsyset.fwhmunits = gipsy.userchar(key, mes, default=1, defval='arcmin') key = "BEAMPARS" + pp beampars = [1.0, 1.0, 0.0] mes = "maj, min (%s), pa (deg) ... [1, 1, 0]:"%gipsyset.fwhmunits bpars = gipsy.userdble(key, mes, default=1, defval=beampars, nmax=3) for i, bp in enumerate(bpars): beampars[i] = bp gipsyset.beampars = beampars key = "BEAMATT" + pp mes = "Plot attributes for plot object ... [alpha:0.5]:" defval = "alpha:0.5" s = gipsy.usertext(key, mes, default=1, defval=defval) gipsyset.beamkwargs = parstr2dict(s) return gipsyset, len(gipsyset.subsets) def openscript(key, mes, taskname): script = None scriptname = gipsy.userchar(key, mes, defval=None, default=1) if not scriptname is None: script = [] commentline = "#"+'-'*70 FILE = open(scriptname, "w") script.append("#!/usr/bin/env python\n\n") from getpass import getuser username = getuser() from datetime import datetime date = datetime.now().strftime("%d-%m-%Y at %H:%M:%S") script.append(commentline+"\n") s = "# Script generated with GIPSY task %s by user: %s\n"%(taskname, username) script.append(s) s = "# Date: %s\n" % (date) script.append(s) script.append(commentline+"\n\n") script.append("from matplotlib import pyplot as plt\n") script.append("from kapteyn.maputils import FITSimage\n\n") script.append("taskname = '%s'\n" %taskname) else: FILE = None script = None scriptname = '' return FILE, script, scriptname def closescript(FILE, script): if FILE is None: return else: script.append("\nplt.show()\n") FILE.writelines(script) FILE.close() return def getfigure(xfig_cm=25, yfig_cm=25, canvascolor= '#efefef', script=None, deflev=1): # Ask user to enter a size for the figure canvas. # The default is the size of A4 is 210 x 297 (mm) cm2inch = 0.393700787 # 1cm is 0.393700787 inch #xfig_cm = 21.0 #yfig_cm = 29.7 key = "FIGSIZE=" mes = "Enter figure size (width,height) in cm ... [%.1f,%.1f]:" % (xfig_cm, yfig_cm) defval = (xfig_cm, yfig_cm) xycm = gipsy.userdble(keyword=key, message=mes, defval=defval, default=deflev, nmax=2) if (len(xycm) == 1): xycm.append(yfig_cm) xfig_inch = xycm[0] * cm2inch; yfig_inch = xycm[1] * cm2inch; # Create a figure key = "CANVASCOL=" mes = "Enter color of figure canvas ... [%s]:"%canvascolor canvascolor = gipsy.usertext(keyword=key, message=mes, defval=canvascolor, default=deflev) fig = figure(figsize=(xfig_inch, yfig_inch), facecolor=canvascolor) # Set a figure for Matplotlib if not script is None: s = "fig = plt.figure(figsize=(%g, %g), facecolor='%s')\n"%(xfig_inch, yfig_inch, canvascolor) script.append(s) # Are there any subplot parameters? key = "SUBPLOTPARS=" mes = "Enter parameters for subplots in param:value format ... [defaults]:" s = gipsy.usertext(keyword=key, message=mes, defval='', default=deflev) subplotpars = None subplotparsstr = '' if s != '': # Result of usertext is a list ss = re_split('[ ,\t]+', s) kwargs = {} for sss in ss: if len(ss) > 0: try: k, v = sss.split(':') kwargs[k] = float(v) # Create string for script. Take care of comma's if subplotparsstr == '': subplotparsstr += '%s=%s'%(k,v) else: subplotparsstr += ', %s=%s'%(k,v) except: pass if len(kwargs) > 0: fig.subplots_adjust(**kwargs) if not script is None: s = "fig.subplots_adjust(%s)\n"%subplotparsstr script.append(s) return fig def plotimages(fig, allsets, n, taskname, filename=None, script=None, deflev=1): # allsets is the list with images & subsets. n is the total number of imagespec # which can be higher than the number of elements in the allsets list, # because each subset is counted as an image. writetoscript = not script is None images = [] # A list with Annotatedimage objects rows = int(round(sqrt(n))) cols = n/rows if n%rows: cols +=1 key = "ROWSCOLS=" mes = "How many rows & columns for %d images? ... [%d,%d]:" % (n, rows,cols) rowscols = gipsy.userint(key, mes, defval=(rows,cols), default=deflev+4, nmax=2) ticks_fontsize = 11 - max(rowscols[0], rowscols[1]) if ticks_fontsize < 4: ticks_fontsize = 4 lab_fontsize = max(4, ticks_fontsize) if writetoscript: script.append("ticks_fontsize = %g\n"%ticks_fontsize) script.append("lab_fontsize = %g\n"%lab_fontsize) key = "SKIPLABELS=" mes = "Skip WCS labels not part of left and bottom axes? ... Y/[N]:" minlabs = gipsy.userlog(key, mes, defval=False, default=deflev) key = "TYPECLI=" mes = "Do you want mouse positions on command line? ... Y/[N]:" typecli = gipsy.userlog(key, mes, defval=False, default=deflev) nrows = rowscols[0] ncols = rowscols[1] # Does user wants some help text at bottom of figure? key = "HELPTEXT" mes = "Print help at bottom of figure ... Y/[N]:" wanthelptext = gipsy.userlog(key, mes, default=deflev, defval=False, nmax=1) r = c = 0 imnr = 1 framenr = 1 for gipsyset in allsets: allim = [] f = FITSimage(externalheader=gipsyset.projection.source, externaldata=gipsyset.image, externalname=gipsyset.dataname) if writetoscript: # Here we assume the image data came from a FITS file. # Currently there is no way to read directly from a GDS image # in the Kapteyn Package environment script.append("\nf = FITSimage('%s')\n"%f.filename) # Set a dummy slice (if any) to get the pixel aspect ratio f.set_imageaxes(gipsyset.projaxnum[0], gipsyset.projaxnum[1], slicepos=gipsyset.slicepositions[0]) gipsy.anyout(taskname + ": Pixel aspectratio %s: %f"%(gipsyset.dataname, f.get_pixelaspectratio())) key = "ASPECTRATIO%d="%imnr aspect = f.get_pixelaspectratio() mes = 'Enter pixel aspect ratio ... [%g]: '%aspect aspect = gipsy.userreal(key, mes, defval=aspect, default=deflev, nmax=1) cmap = gipsyset.colormap if writetoscript: script.append("cmap = '%s'\n" %gipsyset.colormap) for i, subset in enumerate(gipsyset.subsets): if writetoscript: script.append("\nclipmin = %g\n"%gipsyset.vminmax[0]) script.append("clipmax = %g\n"%gipsyset.vminmax[1]) script.append("pxlim = %s\n"%gipsyset.pxlim) script.append("pylim = %s\n"%gipsyset.pylim) if c == ncols: c = 0 r += 1 key = "TITLE%d_%d="%(imnr, i+1) if i == 0: deftitle = gipsyset.dataname.split('/')[-1] # Remove path mes = "Title for current plot ... [%s]:"%deftitle default = deflev else: deftitle = "" mes = "" default = 2 title = gipsy.usertext(key, mes, default=default, defval=deftitle) if title != '': key = "TITLEATT%d_%d="%(imnr, i+1) mes = "Plot attributes for plot object ... []:" defval = "" s = gipsy.usertext(key, mes, default=2, defval=defval) titlekwargs = parstr2dict(s) frame = fig.add_subplot(rows, cols, framenr); f.set_imageaxes(gipsyset.projaxnum[0], gipsyset.projaxnum[1], slicepos=gipsyset.slicepositions[i]) f.set_limits(pxlim=gipsyset.pxlim, pylim=gipsyset.pylim) if writetoscript: #script.append("axesrect = %s\n"%(str(axesrect))) #script.append("frame = fig.add_axes(axesrect)\n") script.append("frame = fig.add_subplot(%d, %d, %d)\n"%(rows, cols, framenr)) script.append("f.set_imageaxes(%d, %d , slicepos=%s)\n"%\ (gipsyset.projaxnum[0], gipsyset.projaxnum[1], gipsyset.slicepositions[i])) script.append("f.set_limits(pxlim=pxlim, pylim=pylim)\n") framenr += 1 if gipsyset.u_skyout != None: f.set_skyout(gipsyset.u_skyout) if writetoscript: script.append("f.set_skyout('%s')\n"%gipsyset.u_skyout) if gipsyset.subspectraltrans != None: f.set_spectrans(gipsyset.subspectraltrans) if writetoscript: script.append("f.set_spectrans('%s')\n"%gipsyset.subspectraltrans) annim = f.Annotatedimage(frame, clipmin=gipsyset.vminmax[0], clipmax=gipsyset.vminmax[1], cmap=cmap, basename=taskname, gridmode=True) annim.set_aspectratio(aspect) annim.Image() annim.cmap.set_scale(gipsyset.colorscale) annim.cmap.modify(gipsyset.cslope, gipsyset.coffset) # Done by constructor? annim.set_pixoffset(gipsyset.subproj.crpix[0], gipsyset.subproj.crpix[1]) if writetoscript: script.append("annim = f.Annotatedimage(frame, clipmin=clipmin, clipmax=clipmax, cmap=cmap, basename=taskname, gridmode=True)\n") script.append("annim.Image()\n") script.append("annim.cmap.set_scale('%s')\n"%gipsyset.colorscale) script.append("annim.cmap.modify(%g, %g)\n"%(gipsyset.cslope, gipsyset.coffset)) grat = None if gipsyset.drawgraticules: spectrans = gipsyset.subspectraltrans startx = gipsyset.gratstart[0] starty = gipsyset.gratstart[1] deltax = gipsyset.gratdelta[0] deltay = gipsyset.gratdelta[1] unitsx = gipsyset.gratunits[0] unitsy = gipsyset.gratunits[1] grat = annim.Graticule(spectrans=spectrans, startx=startx, starty=starty, deltax=deltax, deltay=deltay, unitsx=unitsx, unitsy=unitsy) if writetoscript: if spectrans is None: script.append("spectrans = None\n") else: script.append("spectrans = '%s'\n"%spectrans) if startx is None: script.append("startx = None\n") else: # In this context startx is a string! script.append('startx = "%s"\n'%startx) # Note the quotes. If startx strings contain grouping it can contain # single quotes. if starty is None: script.append("starty = None\n") else: # In this context starty is a string! script.append('starty = "%s"\n'%starty) if deltax is None: script.append("deltax = None\n") else: script.append('deltax = "%s"\n'%deltax) if deltay is None: script.append("deltay = None\n") else: script.append('deltay = "%s"\n'%deltay) if unitsx is None: script.append("unitsx = None\n") else: script.append('unitsx = "%s"\n'%unitsx) if unitsy is None: script.append("unitsy = None\n") else: script.append('unitsy = "%s"\n'%unitsy) script.append("grat = annim.Graticule(spectrans=spectrans,\n") script.append(" startx=startx, starty=starty,\n") script.append(" deltax=deltax, deltay=deltay,\n") script.append(" unitsx=unitsx, unitsy=unitsy)\n") if minlabs: if c == 0: grat.set_tickmode(plotaxis=("left"), mode="Native") grat.setp_axislabel(plotaxis="left", fontsize=lab_fontsize) if writetoscript: script.append('grat.setp_axislabel(plotaxis="left", fontsize=lab_fontsize)\n') script.append('grat.setp_axislabel(plotaxis="left", fontsize=lab_fontsize)\n') else: grat.setp_axislabel(plotaxis="left", visible=False) grat.set_tickmode(plotaxis=("left"), mode="NO") if writetoscript: script.append('grat.setp_axislabel(plotaxis="left", visible=False)\n') script.append('grat.set_tickmode(plotaxis=("left"), mode="NO")\n') if r == nrows - 1: grat.set_tickmode(plotaxis=("bottom"), mode="Native") grat.setp_axislabel(plotaxis="bottom", fontsize=lab_fontsize) if writetoscript: script.append('grat.set_tickmode(plotaxis=("bottom"), mode="Native")\n') script.append('grat.setp_axislabel(plotaxis="bottom", fontsize=lab_fontsize)\n') else: grat.setp_axislabel(plotaxis="bottom", visible=False) grat.set_tickmode(plotaxis=("bottom"), mode="NO") if writetoscript: script.append('grat.setp_axislabel(plotaxis="bottom", visible=False)\n') script.append('grat.set_tickmode(plotaxis=("bottom"), mode="NO")\n') else: grat.set_tickmode(plotaxis=("bottom", "left"), mode="Native") grat.setp_axislabel(plotaxis=("bottom", "left"), fontsize=lab_fontsize) if writetoscript: script.append('grat.set_tickmode(plotaxis=("bottom", "left"), mode="Native")\n') script.append('grat.setp_axislabel(plotaxis=("bottom", "left"), fontsize=lab_fontsize)\n') grat.setp_ticklabel(plotaxis=("left","bottom"), color='b', fontsize=ticks_fontsize) if writetoscript: script.append('grat.setp_ticklabel(plotaxis=("left","bottom"), color="b", fontsize=ticks_fontsize)\n') if gipsyset.wantcont: annim.Contours(levels=gipsyset.levels, colors=gipsyset.colors) if writetoscript: script.append("levels = %s\n"%(str(gipsyset.levels))) s = "colors = (" for col in gipsyset.colors: s += '"'+col+'",'; script.append("%s)\n"%(s)) script.append("annim.Contours(levels=levels, colors=colors)\n\n") if not gipsyset.beampos is None: beam = annim.Beam(major=gipsyset.beampars[0], minor=gipsyset.beampars[1], pa = gipsyset.beampars[2], pos=gipsyset.beampos, units=gipsyset.fwhmunits, **gipsyset.beamkwargs) # Fails on MPL systems v < 0.99.? if writetoscript: script.append('major = %g\n'%gipsyset.beampars[0]) script.append('minor = %g\n'%gipsyset.beampars[1]) script.append('beampa = %g\n'%gipsyset.beampars[2]) script.append('beampos = "%s"\n'%gipsyset.beampos) script.append('fwhmunits = "%s"\n'%gipsyset.fwhmunits) script.append('beamkwargs = %s\n'%str(gipsyset.beamkwargs)) script.append('beam = annim.Beam(major=major, minor=minor, pa=beampa, pos=beampos, units=fwhmunits, **beamkwargs)\n') if gipsyset.gridlabels: annim.Pixellabels(plotaxis=("top","right"), **gipsyset.gridlabelkwargs) if writetoscript: script.append('gridlabelkwargs = %s\n'%str(gipsyset.gridlabelkwargs)) script.append("annim.Pixellabels(plotaxis=('top','right'), **gridlabelkwargs)\n") if gipsyset.wantcolorbar: cb = None if writetoscript: script.append("cborientation = '%s'\n"%gipsyset.cborientation) if gipsyset.colbarframe is None: cb = annim.Colorbar(orientation=gipsyset.cborientation, fontsize=lab_fontsize-1) if writetoscript: script.append("cb = annim.Colorbar(orientation=cborientation, fontsize=lab_fontsize-1)\n") else: # The colorbar frame is given as a center position (xc, yc) # a width w and height h. The coordinates are normalized coordinates # in the current frame. They are converted to normalized coordinates # with respect to the figure, because this is needed to create # a Matplotlib Axes object. xc, yc, w, h = gipsyset.colbarframe x0 = xc-w/2.0; y0 = yc-h/2.0; x1 = xc+w/2.0; y1 = yc+h/2.0 mode = gipsyset.cbframemode single = False if mode.upper().startswith('F'): X0,Y0 = frame.transAxes.transform((x0,y0)) X1,Y1 = frame.transAxes.transform((x1,y1)) x0,y0 = fig.transFigure.inverted().transform((X0,Y0)) x1,y1 = fig.transFigure.inverted().transform((X1,Y1)) elif mode.upper().startswith('D'): X0,Y0 = frame.transData.transform((x0,y0)) X1,Y1 = frame.transData.transform((x1,y1)) x0,y0 = fig.transFigure.inverted().transform((X0,Y0)) x1,y1 = fig.transFigure.inverted().transform((X1,Y1)) else: single = True if not single or i == 0: f4 = (x0, y0, x1-x0, y1-y0) adjust = frame.get_adjustable() cfr = fig.add_axes(f4, adjustable=adjust, autoscale_on=False) cb = annim.Colorbar(orientation=gipsyset.cborientation, fontsize=lab_fontsize-1, frame=cfr) if writetoscript: script.append("adjust = frame.get_adjustable()\n") script.append("cfr = fig.add_axes(%s, adjustable=adjust, autoscale_on=False)\n"%str(f4)) script.append("cb = annim.Colorbar(orientation=cborientation, fontsize=lab_fontsize-1, frame=cfr)\n") if not cb is None and gipsyset.cbunits != '': cb.set_label(gipsyset.cbunits, fontsize=lab_fontsize-1) if writetoscript: script.append("cbunits = '%s'\n"%gipsyset.cbunits) script.append("cb.set_label(cbunits, fontsize=lab_fontsize-1)\n") if title != '': # First plot of a subset sequence # Titlekwargs belong to one subset. Therefore they are # not an attribute of 'gipsyset'. frame.set_title(title, **titlekwargs) if writetoscript: script.append('title = "%s"\n'%title) script.append('titlekwargs = %s\n'%str(titlekwargs)) script.append("frame.set_title(title, **titlekwargs)\n") if writetoscript: script.append("\n\n") # Add space to distinguish plots annim.plot() if writetoscript: script.append("annim.plot()\n") # If more than one subset was entered, then give world coordinates # for the slices. if len(gipsyset.subsets) > 1: # and not (gipsyset.subspectra is None): # If one of the axes outside the map is spectral, then ask # user to enter a spectral translation. vel, uni = f.slice2world(spectra=gipsyset.repspectraltrans, userunits=gipsyset.repspectralunits) # Skyout is also possible. sliceinfo = "" for v, u in zip(vel, uni): sliceinfo += "%g (%s) " % (v,u) if not grat is None: fr = grat.frame # This frame is only known after annim.plot() which sets frames else: fr = frame fr.text(0.98, 0.98, sliceinfo, horizontalalignment='right', verticalalignment='top', transform = fr.transAxes, fontsize=ticks_fontsize-1, color='w', bbox=dict(facecolor='red', alpha=0.7)) if writetoscript: if gipsyset.repspectraltrans is None: s = "repspectra = None\n" else: s = "repspectra = '%s'\n"%gipsyset.repspectraltrans script.append(s) if gipsyset.repspectralunits is None: s = "userunits = None\n" else: s = "userunits = '%s'\n"%gipsyset.repspectralunits script.append(s) script.append("vel, uni = f.slice2world(spectra=repspectra, userunits=userunits)\n") script.append('sliceinfo = ""\n') script.append('for v, u in zip(vel, uni):\n') script.append(' sliceinfo += "%g (%s) " % (v,u)\n') if not grat is None: script.append('fr = grat.frame\n') else: script.append('fr = frame\n') script.append("""fr.text(0.98, 0.98, sliceinfo, horizontalalignment='right', verticalalignment='top', transform = fr.transAxes, fontsize=ticks_fontsize-1, color='w', bbox=dict(facecolor='red', alpha=0.7))\n""") annim.interact_toolbarinfo() annim.interact_imagecolors() annim.interact_writepos(gipsy=True, typecli=typecli) if writetoscript: script.append("annim.interact_toolbarinfo()\n") script.append("annim.interact_imagecolors()\n") # Script runs outside GIPSY. But it must use grid coodinates # and not pixel coordinates. Therefore we set a pixel offset for # mouse positions. We set the projection object in grid mode # and tell the output of tagged positions to write to # a terminal instead of the GIPSY log and screen. script.append("annim.interact_writepos(gipsy=False, typecli=False)\n") c += 1 allim.append(annim) # All Annotatedimage objects per 'gipsyset' object images.append(annim) # List with all images of all 'gipsyset' objects cmap = annim.cmap # Use same colormap instance for all subsets if writetoscript: script.append("cmap = annim.cmap\n") imnr += 1 gipsyset.allim = allim # Later we need to be able to identify the subsets if wanthelptext: tdict = dict(color='g', fontsize=10, va='bottom', ha='left') helptxt = allsets[0].allim[0].get_colornavigation_info() fig.text(0.01, 0.01, helptxt, tdict) return images # Return list with all images of all 'gipsyset' objects #---------- Main program ----------------- gipsy.init() taskname = gipsy.myname() # Start prompting user for data sets askmore = True inp = 0 n = 0 # Counter for input of data sets. allsets = [] key = "DEVLEV=" mes = "Prompt all keywords (1) or only most important (2)? ... [2]:" deflev = gipsy.userint(key, mes, defval=2, default=1) if deflev != 2: deflev = 1 while askmore: inp += 1 gipsyset, plotted = imagespec(inp, taskname, deflev) askmore = plotted > 0 n += plotted if askmore: allsets.append(gipsyset) if n == 0: # No images were specified gipsy.finis() # Leave Hermes # Ask user of a name for a Python script that does the actual plotting. key = "SCRIPT=" mes = "Name of plot script ... [do not create script]:" fileobj, script, scriptname = openscript(key, mes, taskname) fig = getfigure(xfig_cm=40, script=script, deflev=deflev) images = plotimages(fig, allsets, n, taskname, script=script, deflev=deflev) closescript(fileobj, script) if not fileobj is None: gipsy.anyout(taskname + ": Wrote Python script [%s] to file"%scriptname) # Ask some keywords only for WCSFLUX if taskname.upper() == 'WCSFLUX': key = 'POLFROMFILE=' prompt = "Do you want to read polygon data from file? .... Y/[N] " fromfile = gipsy.userlog(keyword=key, message=prompt, defval=False, default=1) filename = None world = False if fromfile: prompt = "Enter name of file with polygon data:" key = "POLYFILENAME=" s = gipsy.usertext(keyword=key, message=prompt, defval='', default=1) if s == '': gipsy.anyout("No file name entered!") else: filename = s prompt = "Read Pixel- or World coords from shape file ... P/[W]:" key = "PIXORWORLD=" s = gipsy.usertext(keyword=key, message=prompt, defval='', default=1) world = False if s == '' or s in ['w', 'W']: world = True key = 'WCS=' prompt = "Markers propagate in Pixel- or World coordinates ... P/[W]:" s = gipsy.userchar(keyword=key, message=prompt, defval='W', default=1) if s.upper() == 'W': wcs = True else: wcs = False for ima in images: ima.world = wcs shapes = shapes.Shapecollection(images, fig, wcs=wcs, inputfilename=filename, inputwcs=world, gipsy=True) # Give user an option to set marker symbols at one or more positions gipsy_connect() gipsy.KeyCallback(plmarker, 'MARKER=', allsets=allsets) show() # At the end of the program, after closing the Matplotlib window, # we write the current settings of the image editor to the # corresponding keywords so that it will be possible to restore # the last plot using the GIPSY macro mechanism. for i, gset in enumerate(allsets): colormap0 = gset.allim[0].cmap gipsy.anyout(taskname + ": CSLOPE%d= %g COFFSET%d= %g"%(i+1, colormap0.slope, i+1, colormap0.shift)) gipsy.wkey("CSLOPE%d= %g"%(i+1, colormap0.slope)) gipsy.wkey("COFFSET%d= %g"%(i+1, colormap0.shift)) scales = {'LINEAR' : '1', 'LOG' : '2', 'EXP' : '3', 'SQRT' : '4', 'SQUARE': '5'} sc = colormap0.scale sc = scales[sc] gipsy.anyout(taskname + ": COLORSCALE%d= %s"%(i+1, sc)) gipsy.wkey("COLORSCALE%d= %s"%(i+1, sc)) # Module mplutil stores the source of the VariableColormap instance # In the context of this task, the source is (or should be) a string. Note that # for a series of subsets only the first stores the string, the others store # a VariableColormap instance. To be sure, check always whether the source is a string. cm = colormap0.source if type(cm) == StringType: gipsy.wkey("COLORMAP%d= %s"%(i+1, cm)) # The parameters that adjust the subplots key = "SUBPLOTPARS=" key += "left:%f right:%f top:%f bottom:%f hspace:%f wspace:%f"%(fig.subplotpars.left, fig.subplotpars.right, fig.subplotpars.top, fig.subplotpars.bottom, fig.subplotpars.wspace, fig.subplotpars.hspace) gipsy.wkey(key) cm2inch = 0.393700787 # 1cm is 0.393700787 inch xs, ys = fig.get_size_inches() gipsy.wkey("FIGSIZE=%f %f"%(xs/cm2inch, ys/cm2inch)) # End GIPSY task gipsy.finis()