Program: WCSFLUX Purpose: Plot set/subsets from a GIPSY set or slices from a FITS file and annotate those plots with world coordinate labels and or a graticule. Then add shapes in one image. The shapes are copied to other images using their world coordinates. Finally calculate the flux in those shapes. Category: COORDINATES, DISPLAY, PLOTTING, ANALYSIS File: wcsflux.src Author: M.G.R. Vogelaar Keywords: WCSFLUX specific keywords: ========================== POLFROMFILE= Do you want to read polygon data from file? .... Y/[N] Data used to construct a shape can be written to an Ascii file on disk. If you want to read it back, specify the file name with keyword POLYFILENAME= Shape data is written by pressing key 'w'. A file with shape data is written to disk with a file name that starts with 'shapes'. It also contains a date/time stamp. This name cannot be changed by a user. After reading a file with shape data, the plotting the data is started after pressing key 'r'. POLYFILENAME= Enter name of file with polygon data: Name of file on disk with data that defines a shape. Here is an example of a file with shape data. The data corresponds to an image with a world coordinate system and should be read back in a compatible system. The third a fourth column set a pixel coordinate. These start with 1 (FITS). Note that this is not the same as GIPSY grids (where grid 0 is equal to CRPIX). Column five and six are world coordinates given in the system of the image from which they are written to file. The first two columns are for internal administration. They indicate the associated shape (one of polygon, ellipse, circle rectangle or spline). The second number is the number of the object in its shape category. ! Saved at Sunday 28/11/2010 14:39:34 ! Data saved for image: ../Mapmaker/ngc6946.fits ! Slice position(s): [51] ! Format of this file: ! Shape # - object# - x pixel - y pixel - x world - y world ! --------------------------------------------------------- 0 0 84.867303 29.604512 40.707425 28.550258 0 0 138.401389 28.117454 40.687110 28.549761 0 0 138.401389 78.677425 40.687108 28.566615 0 0 138.401389 129.237395 40.687106 28.583468 0 0 80.406129 128.493866 40.709121 28.583221 0 0 84.123774 77.933896 40.707708 28.566367 0 0 84.123774 27.373925 40.707707 28.549514 PIXORWORLD= Read Pixel- or World coords from shape file ... P/[W]: Select for the file in POLYFILENAME= whether you want the positions in pixel coordinates (P) or in world coordinates (W). Note that pixel coordinates follow the FITS standard, i.e. they start with 1. These are not the same as grid coordinates in GIPSY because there we label the pixel that corresponds to the projection center (CRPIX) as grid 0. The advantage is that one can use the shape data files in both the GIPSY environment and in applications that follow the FITS standard. WCS= Markers propagate in Pixel- or World coords ... P/[W]: While plotting a marker at a position in the current image, the program looks for other images and tries to plot a corresponding position in those images. The correspondence can either be in pixel coordinates or in world coordinates. The procedure also works for maps with only one spatial axis (e.g. position-velocity diagrams). The graphical user interface (gui) shows a number of buttons. The text on the buttons explain what they do. At the bottom of the gui you will see a list with (keyboard) keys and their associated action. Keywords common in MAPMAKER and WCSFLUX ======================================= DEVLEV= Prompt all keywords (1) or only most important (2)? ... [2]: Many keywords are to adjust plot parameters and it is not useful to prompt all of them if you only want to have a quick plot on the screen for inspection. If you need different sky systems, spectral translations or an optimal result for a hardcopy, then set DEFLEV=1. Keywords in a loop: ------------------- Keywords prompted in a loop have a name followed by a number like KEYWORD1=, KEYWORD2= etc. The number corresponds to the data set which is also entered in a loop. DATASET1= Enter name of FITS file or GIPSY set ... [list files]: Specify the location of the data. This data can either be a GIPSY set or a FITS file. You can give the name of the file, then other keywords to specify axes will follow. Or your input contains the description of the slices (subsets). This syntax is the same for FITS files and GIPSY sets. Examples: 1) DATASET=AURORA 2) DATASET=AURORA FREQ 10:20 30 40 3) DATASET=m101.fits 4) DATASET=http://people.sc.fsu.edu/~burkardt/data/fit/datacube.fit FELO 2 5) DATASET=datacube.fit.gz FELO 0 ad 1) Just the name of the data. The program will prompt for additional information. Knowledge about the axes in the set is not necessary. ad 2) Using knowledge of the axes in a set, specify location and subsets (data slices). ad 3) Same as 1 but now for a FITS file ad 4) FITS files can be retrieved from an Internet location with an URL In this example one slice is defined. ad 5) FITS files can also be accessed in gzipped format. Note that the pixels for FITS files are not 1-based as is standard in FITS, but they follow GIPSY rules i.e. pixel 0 is found at the header value of 'crpix' of the corresponding axis. HDU1= Enter number of Header Data Unit ... [0]: If a FITS file is selected and the file contains more than one Header Data Units (hdu), the user is prompted to select the unit where the image data is stored. A table with options is written to screen. The input is a valid integer number. The default is 0. AXNAMES1= Enter two of RA,DEC,VELO ... [RA DEC]: If you did not specify subsets and the dimension of your data is greater than 2, you are prompted to enter two axis names that are part of the image you want to plot. In the prompt all the axes are included. The default are the first two axes in your set. The input are two strings which are minimal matched.against the known axes. Note that you can enter the axes in any order. Examples for set AURORA with RA, DEC, FREQ axes. AXNAMES1= # default is RA,DEC AXNAMES1=RA FREQ # XV map with FREQ along vertical axis AXNAMES1=FR R # XV map with FREQ along horizontal axis AXNAMES1=D R # Spatial DEC,RA map XXXX1= Enter pixel position(s) on XXXX ... [lo:hi]: XXXX represents a variable set of keywords with the same name as the repeat axes in a set. Repeat axes are all axes that are not part of the specified data slice (AXNAMES1=). The default values are all pixels along the axis resulting in all possible slices along the repeat axis. The syntax follows the rules for GIPSY input of integers. Example for set AURORA2 with: FREQ-OHEL from -4 to 5 DEC--NCP from -29 to 20 RA---NCP from -24 to 25 POL from 0 to 3 AXNAMES=DEC FR Enter grid position(s) on RA ..... [-24:25] RA=0 Enter grid position(s) on POL ..... [0:3] POL=1 This defines one slice with axes DEC and FREQ at RA=0 and POL=1 or: RA=-2:2 10 to get six slices at RA -2,-1,0,1,2 and 10 Note that these subset/slices could also be entered at the DATASET1= prompt as follows: DATASET1=AURORA2 ra -2:2 10 pol 1 SWAPAXES1= Do you want to swap axes RA and DEC? ... Y/[N]: If your input defined 2D slices then you you are able to swap the two axes in the subset. The GIPSY input syntax for set/subsets does not provide syntax to get for swapped axes. BOX1= Enter limits of map xlo,ylo, xhi,yhi ... [xlo,ylo, xhi,yhi]: Enter coordinates that extract an area smaller than or equal (the default) to the slice area. Pleas note: 1) that the syntax follows that of GIPSY and therefore the pixel positions are not 1-based as in FITS. The GIPSY syntax allows for a center coordinate and sizes. 2) that world coordinates are related to GIPSY's legacy coordinate system which will conflict with WCSlib. Therefore only grid coordinates and no world coordinates should be entered. 3) the order of the axes is always that of the axes in the data set. even if you forced swapping axes with AXNAMES= or SWAPAXES= Examples: Limits RAlo,VELOlo, RAhi,VELOhi ... [-50,-50,49,49]: BOX= BOX=-9 -9 10 10 BOX=0 0 D 20 20 Note that the second and third example result in the same area. SPECTRANS1= Spectral translation for axis xxx or ? for help ... [native]: If one of the image axes is of spectral type, you are prompted to enter a so called 'spectral translations' which sets a translation from the native coordinate to another (e.q. frequency to velocity). If you entered ? then you will see a list with available spectral translation All translations known by WCSlib are tested by the program and only the valid ones are listed. One can enter part of the string that represents the spectral translation. Example: SPECTRANS1=? =========== Available spectral translations =========== FREQ (Hz), ENER (J), WAVN (/m), VOPT-F2W (m/s), VRAD (m/s), VELO-F2V (m/s) WAVE-F2W (m), ZOPT-F2W, AWAV-F2A (m), BETA-F2V SPECTRANS1=b This selects the spectral translation BETA-F2V which has no units. Labeling and position information in the plot are from this moment on related to this translation. SUNITS1= Units for spectral slice info ... [native]: If more than one subset/slice is entered and one of the subset/slice axes is spectral, then you are asked to give an (alternative) unit for the numbers in the subset/slice info label. If, for example, the native axis is frequency and you entered SPECTRANS1=VOPT, then with SUNITS1=km/s you get the slice information optical velocities in km/s. COLORMAP1= Enter name of color map or ? for list ... [jet]: If you enter ? then a list with color maps is displayed in the Hermes' log file. Your selection is (part of) a string and is minimal matched with the entries in the list. If it is not in the list it could be a color map (lut) on disk. Color luts are created with key 'm' when an image is displayed. This way, one can create custom color maps. Note that the input is case sensitive. COLORSCALE1= 1=linear 2=log 3=exp 4=sqrt 5=square ... [1]: Set the color scaling to one of the list in the prompt. CSLOPE1= Slope image values vs colors ... [1.0]: Together with COFFSET1=, one defines how the range of image values correspond with a range of colors. Usually the slope and offset COFFSET1= Offset image values vs colors ... [0.0]: Together with CSLOPE1=, one defines how the range of image values correspond with a range of colors. VMINMAX1= Scale colors between image values ... [datamin, datamax]: Before the prompt, the minimum and maximum data values of you slices (subsets) within the given axes limits (BOX=) are calculated. These values are used as a default to distribute the colors of your color map. Data outside this range is clipped. WANTCONT1= Add contours? ... Y/[N]: To enhance the visibility of certain features in a map, one can add contour lines. If you confirm wanting contours you will be prompted for the wanted contour levels with CLEVELS= CLEVELS1= Enter contour levels between .. and .. ... [automatic]: Enter one or more levels in units of you image data for which contour lines will be plotted. The default forces the program to calculate levels. Usually these are not the most suitable levels because also noise is represented by one or more of these levels. The contours are labeled. CCOLORS1= Enter one or more colors for contours ... [r]: Distinguish contours by a color. The default is one color (red). In the log file a list is presented with possible colors and their abbreviation. WANTCOLORBAR1= Add a colorbar? ... Y/[N]: Create a bar with colors from the selected colormap and label these colors with the corresponding data values in the units of your data. COLBARFRAME1= Enter x0,y0, width, height ... [calculated]: If you want to specify a frame for a colorbar, then enter a center and a width and height. The coordinates in which this must be entered, depend on the value of CBFRAMEMODE1= CBFRAMEMODE1= Frame coords: F=Norm. in frame, D=Data coords ... [Norm. in figure]: The coordinates entered in COLBARFRAME1= are in normalized figure coordinates, i.e. the lower left coordinate of the figure is (0,0) and the upper right coordinate of the figure is (1,1). This is the default. If CBFRAMEMODE1=F then the coordinates are normalized coordinates but now with respect to the current subplot. If CBFRAMEMODE1=D, the coordinates are pixel coordinates with respect to the current subplot (Not yet tested and perhaps less useful). CBORIENT1= Horizontal color bar (H) or vertical (V)? ... H/[V]:" Orientation of the colorbar. CBLABEL1= Label for color bar ... [BUNITS from header]: To print a label near the colorbar with information about the data units, one should enter a string here. The default string is a unit as read from the header of Set or FITS file (FITS keyword BUNIT) SKYOUT1= Sky system 0=eq, 1=ecl, 2=gal, 3=sup.gal ... [native]? If both axes of your map are spatial axes (one longitude and one latitude axis), then you are prompted to enter a celestial system (sky system) for the coordinate transformations. The default value is the system native to the data as defined by the FITS or GIPSY header. The systems are given by number: 0 = Equatorial 1 = Ecliptic 2 = Galactic (II) 3 = Supergalactic Equatorial systems are further identified by an Equinox and a reference system. Those can be entered with the keywords EQUINOX= and REFSYS= You can only change the output sky system for spatial axes and not for maps with one spatial axis and one non spatial axis (e.q. position-velocity maps). RADESYS1= Ref.sys 0=fk4, 1=fk4_no_e, 2=fk5, 3=icrs, 4=j2000 ... [native]: Equatorial and ecliptic coordinates are relative to a reference system. The system is selected with an integer number. The fk4_no_e system is fk4 where the positions are not corrected for the elliptic terms of aberration. If nothing appropriate could be found in the header a reference system is assumed according the FITS standard rules. EQUINOX1= Enter equinox (eq. J2000, B1983.5) ... [native]: Equatorial and ecliptic coordinates are relative to a given equinox. The default value is the one found in the header. If the value is not listed in the header of your data set, an equinox is assumed according the FITS standard rules. The input is an integer or floating point number prefixed by one or more characters to set an epoch. B Besselian epoch. Example 'B 1950', 'b1950', 'B1983.5', '-B1100' J Julian epoch. Example: 'j2000.7', 'J 2000', '-j100.0' JD Julian date. This number of days (with decimals) that have elapsed since the initial epoch defined as noon Universal Time (UT) Monday, January 1, 4713 BC in the proleptic Julian calendar Example: 'JD2450123.7' MJD The Modified Julian Day (MJD) is the number of days that have elapsed since midnight at the beginning of Wednesday November 17, 1858. In terms of the Julian day: MJD = JD - 2400000.5 Example: 'mJD 24034', 'MJD50123.2' RJD The Reduced Julian Day (RJD): Julian date counted from nearly the same day as the MJD, but lacks the additional offset of 12 hours that MJD has. It therefore starts from the previous noon UT or TT, on Tuesday November 16, 1858. It is defined as: RJD = JD - 2400000 Example: 'rJD50123.2', 'Rjd 23433' F 1) DD/MM/YY Old FITS format Example: 'F29/11/57' 2) YYYY-MM-DD FITS format Example: 'F2000-01-01' 3) YYYY-MM-DDTHH:MM:SS FITS format with date and time. Example: 'F2002-04-04T09:42:42.1' GRATICULES1= Draw graticule lines and WCS labels? ... [Y]/N: A world coordinate system with parallels and meridians is plotted labeled with the corresponding constant values. The system is plotted in the projection defined in the header and in the celestial system defined with previous keywords SKYOUT=, RADESYS= and EQUINOX= If one wants to plot images without the graticule overlay, then enter GRATICULES=N GRIDLAB1= Do you want to plot grid labels? ... Y/[N]: Annotate the grids along the top and right plot axis. Note that grid (0,0) corresponds to pixel NINT(CRPIX_x), NINT(CRPIX_y). Function NINT() finds the nearest integer. CRPIX is the reference pixel. GRIDLABATT1= Plot attributes for plot object ... [color:g]:" Plot attributes in format key:value and value always without quotes. E.g.: GRIDLABATT1=color:#aabb00 fontsize=4 BEAMPOS1= Enter position of beam ... [no beam]: If a position is entered, the program assumes you want to plot a beam and it will ask for related parameters. The position is a string that follows the syntax for positions as described in wcspos.dc2 E.g. BEAM=23h10m30.2s -10d34m23s FWHMUNITS1= Units of FWHM major&minor axes ... [arcmin]: If a position is given in BEAMPOS=, one can enter the units in which the Full Width at Half Maximum are entered. The units must be compatible to degrees. BEAMPARS1= Maj, min (units), pa (deg) ... [1, 1, 0]: FWHM of major axis, minor axis and position angle of the beam. The units of the major and minor axis are given in keyword FWHMUNITS=. The angle is w.r.t. the North and is given in degrees. BEAMATT1= Plot attributes for plot object ... [alpha:0.5]: Plot attributes in format key:value and value always without quotes. e.g. BEAMKWARGS1= alpha:0.5 hatch:/ color:r ROWSCOLS= How many rows & columns for N images? ... [i,j]: Set a layout for a mosaic of plots. In the prompt you will see the total number of plots set by your DATASET= keywords and a suggestion is made for the number of rows and columns. SKIPLABELS= Skip WCS labels not part of left and bottom axes? ... Y/[N]: For mosaics of plots the labels can get in the way. For a better looking array of plots you can decide to get rid of (most) abundant labeling with SKIPLABELS=Y TYPECLI= Do you want mouse positions on command line ... Y/[N]: If set to Y, then clicking with the left mouse button in the displayed image will write the mouse position in grids on the Hermes command line. ASPECTRATIO1= Enter pixel aspect ratio ... [derived from pixel size]: This task tries to plot pixels in the correct size as given in world coordinates in the header (these sizes are given by FITS keywords CDELT). If one plots a map with uncorrelated axes, then this aspect ratio is not physical and a new, more realistic value is calculated. Often one wants to change this value to make a plot more or less square. TITLE1_m= Title for current plot ... [derived from file name] For each plot a title can be plotted. For a series of subsets the title for the first plot can be entered at a prompt, the other subset plot titles are asked with a hidden prompt. So if you want to change them please enter a title on the command line, e.g.: TITLE1_2=second subset 'm' is the plot number which starts with 1 for the first subset. TITLEATT1_m= Plot attributes for plot object ... []: Plot attributes in format key:value and value always without quotes. Handy attributes for a title are color, fontsize, y. The last one sets the position of the title in the y direction in normalized device coordinates, e.g. TITLEATT1_1=y:1.02 FIGSIZE= Enter figure size (width,height) in cm ... [25.0,25.0]: A figure size in centimeter. CANVASCOL= Enter color of figure canvas ... [#efefef]: The background of the figure can be changes with this keyword. SUBPLOTPARS= Enter parameters for subplots in param:value format ... [defaults]: The parameters are all in normalized figure coordinates: left Space in normalized coordinates to the left bottom Space in normalized coordinates to the bottom right Space in normalized coordinates to the right top Space in normalized coordinates to the top wspace The amount of width reserved for blank space between subplots hspace The amount of height reserved for white space between subplots If you change any of these parameters with the subplot configuration tool on the Matplotlib toolbar, then the new values are written back into keyword SUBPLOTPARS= This enables you to retrieve the same plot after using the GIPSY macro mechanism. SCRIPT= Name of plot script ... [do not create script]: If a name is entered, then a Python script with this name will be written on disk in the current directory. The file contains the commands to re-create the plot outside the context of GIPSY. The script can be used as a template. This template can be extended with features from the Kapteyn Package that are not available in the GIPSY context. MARKER= - Enter 'taskname MARKER=' on the Hermes command line and enter one or more positions. Markers follow the syntax as described in GIPSY module IOUTILS. Example: mapmaker marker=10 7.2 5.0 13.8 mapmaker marker={} 20h36m34.71s {} 59d57m21.4s Adjusted keywords for use in macro ================================== Some keyword values can be changed interactively while the plot is displayed. When the task is closed correctly then the corresponding keywords are updated too. Then you can retrieve the same plot with GIPSY's macro mechanism (!taskname on Hermes command line). The keywords are: CSLOPE= (after rescaling colors with mouse) COFFSET= (after rescaling colors with mouse) COLORSCALE= (after using a key to change the color scale) COLORMAP= (after page up/down for new colormap) SUBPLOTPARS= (after using subplot configuration tool) FIGSIZE= (after resizing the window) Notes: The images defined with the previous keywords are displayed on a Matplotlib figure canvas which is part of a new window showing a number of icons (for zooming, panning etc). (current URL: http://matplotlib.sourceforge.net/users/navigation_toolbar.html ) The aspect ratio of the image is taken so that the pixels are displayed with the correct sizes in world coordinates (if possible). If you resize the window, this aspect ratio remains constant. One of the toolbar buttons pops up a window that allows for the selection of the format of a hard copy. The standard figure size allows for a (encapsulated) PostScript output with everything within the bounding box. If one resizes the windows you get PostScript output with plots with a cutoff. PDF output scales correctly. So it is always possible to get the layout that you want by resizing the window. Then create a hardcopy in PDF and convert the PDF to PostScript with a utility (e.q. pdf2ps) Cursor: The cursor in a frame (where a plot is visible) generates an informative message in the toolbar of the figure canvas. This message shows the current grid position and (if possible) the corresponding world coordinate together with the image value at that position. Cursor + shift-left button: Position information is send to the log file. Cursor + right button: Change color map. Important: The markers propagate in the world coordinate system. Sometimes this gives unexpected results. Imagine you display a mosaic of position velocity diagrams extracted from a cube with RA, DEC and a FREQ axis. If the RA/DEC data is rotated (CROTA in the header is unequal to 0, e.g. 30 deg), then the spatial axes labels (plotted along the axes) are not at the same position in different slices. Therefore each slice/subset shows the marker at a different grid position because only the position in world coordinates remains unaltered. For these maps (with one spatial axis) one can mark features in an image using the right mouse button. These markers are copied on the same pixel position in each image that has the same projection system (otherwise we are comparing positions that are not related). Difference between pixel coordinates and grids: FITS defines it coordinates as follows. The first pixel is labeled 1. The last pixel is labeled NAXISn (where n is the axis number). The header includes for each axis a floating point position CRPIXn which tells us where on the axis you can put the label CRVALn, the world coordinate. Historically GIPSY users are used to another system. In GIPSY we define a coordinate system. Its origin (grid 0) is exactly located on CRPIXn. Therefore a box often contains negative numbers which is not the case if you would enter limits in a pixel system. Note that a grid system is not centered on a pixel if the value of CRPIXn is not an integer value. Different approach: This application is a demonstration of how a GIPSY application can be enhanced with standard libraries and Python modules to fulfill modern requirements for reading FITS files, plotting images and dealing with world coordinates. 1) It is a non gui (let's say classic GIPSY) application where keywords are prompted in a sequence. We demonstrated that it is easier for a user to enter the name of a data set first and then to enter the axes names of the subset (which we also call a slice) and to enter the grids on the repeat axes (axes that do not belong to the subset). 2) FITS files are usable as data sets. They are mapped onto a GIPSY set, transparent and invisible to the user. Also supported are FITS files addressed with an URL and gzipped files. 3) The use of WCSlib allows users to get precise information about positions in world coordinates in a map. There is support for spectral translation (e.g. a frequency axis in your image can be labeled in optical, radio or radial velocity etc.) For FITS files with non FITS extensions like FREQ-OHEL etc. the header is parsed to get a set of items that follow the FITS standard. There is also support for reference systems for equatorial and ecliptic celestial systems. We wrote a Python module, called celestial, which contains definitions and algorithms from the Explanatory Supplement to the Astronomical Almanac. Also the syntax for entering an epoch for the equinox is enhanced. Updates: Nov 27, 2010: VOG document created. Nov 11, 2010: VOG update to latest Kapteyn Package support Feb 20, 2011: VOG update documentation Sep 07, 2013: VOG adapted for Matplotlib 1.3