Program: FITSREPROJ Purpose: This task if a versatile tool to re-project data in a FITS file to another- or a modified world coordinate system. The output can be converted to a GIPSY set or a legacy GIPSY set (to be used for tasks like ROTCUR, ROTMOD etc.). Category: FITS, COORDINATES, HEADER, MANIPULATION File: fitsreproj.src Author: M.G.R. Vogelaar Keywords: LEGACY= Is final goal a legacy GIPSY set? ... Y/[N]: If the end product is a GIPSY set (FINAL=G), then a number of restrictions are applied (e.g. maximum number of pixels in output or the projection system). The set should be suitable to be used in legacy GIPSY tasks (ROTCUR, ROTMOD etc.). If you select 'F', there are no restrictions applied. You still can create a GIPSY set with hidden keyword MAKEGDS=Y. The header will be transferred unmodified. The result is a GIPSY set with a header that describes a world coordinate system that can be processed by a new generation of GIPSY applications based on the library WCSLIB. It can also be used in other GIPSY applications but unknown projection systems are all regarded as Gnomonic which is the default system if GIPSY doesn't recognize the projection. Be warned that an imported header can cause invalid spectral coordinate transformations. This task does not process other axis types than spatial. RFITS does some spectral axis parsing, FITSREAD not at all. FITSFILE= Enter name of FITS file ... [list files]: Name of a FITS file or a URL. The file can be a gzipped FITS file. E.g.: FITSFILE= http://www.atnf.csiro.au/people/mcalabre/data/WCS/1904-66_ZPN.fits.gz FITSFILE=/Mydata/m101.fits HDUNR= Enter number of Header Data Unit ... [0]: Prompt will only appear if the FITS file has multiple header data units (HDU's). An HDU overview is printed first. OPTION= Enter option from list above: ... [align to North]: 0. Re-project the data to the WCS of a second FITS file 1. Rotate data over a given angle (it has to create a 'classic' header first. This re-projection removes skew. 2. Rotate image so that it is aligned to the north. This re-projection removes skew. 3. Rotate to given angle. This re-projection removes skew. The value of FITS keyword CROTAn will be the given angle. 4. Give a number of FITS keywords and values which are inserted in the current header and re-project the data to the modified header. You can select sky definition and projection type from a list. 5. Make classic header and leave the data unchanged (i.e. do NOT re-project). You lose the WCS representation of skew in the destination header. Read the notes for more information. FITSDEST= Enter name of FITS file ... [list files]: Only for OPTION=0 A second FITS file with a header with spatial information. The current data will be re-projected so that it will be described by the WCS of this destination header. Only the axes limits can be changed. HDUNRDEST= Enter number of Header Data Unit ... [0]: Only for OPTION=0 Prompt will only appear if the destination FITS file has multiple header data units (HDU's). An HDU overview is printed first. SKYSYS= Sky system 0=eq, 1=ecl, 2=gal, 3=sup.gal .... [current]: Enter a sky system for the destination header. If you want to change the reference system, equinox or epoch, then you should enter 0 or 1. Otherwise the corresponding prompts are skipped because the program assumes that you want to copy the current sky system. REFSYS= Ref.sys 0=fk4, 1=fk4_no_e, 2=fk5, 3=icrs ... [current]: Enter a reference system for the destination header. Note that GIPSY doesn't support sky systems. EQUINOX= Enter equinox (e.g. J2000 or B1983.5) .... [current]: 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.', 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 EPOCH= Enter date of observation (e.g. MJD24034) .... [current]: Date of observation. Only prompted if FK4 or FK4-NO-E was selected. For the prefix, see also description at EQUINOX= For GIPSY sets we remove FITS keyword EPOCH (and if necessary add FITS keyword EQUINOX with a value). The reason is that GIPSY handles EPOCH and EQUINOX as if they were identical. Values for GIPSY keyword EPOCH= will be written into the destination header with FITS keyword MJD-OBS. PROJID= Projection of destination ... [copy from current]: Identifier for the output projection system. A menu is presented in a table. The first column is the identifier. If you required a legacy GIPSY set (LEGACY=Y), then the menu displays only those systems that are compatible to GIPSY's projection systems. PVi_j= Projection parameter PVi_j ... [abort parameter loop]: Some projection systems have extra parameters. You are prompted for these parameters with a reasonable default. WCSKEYS= Enter wcs keys in key:val format ... []: At this prompt you enter the FITS keywords and values for the destination header. Note that this is not trivial because it depends on the header which keywords one can use. Examples: -Change the scaling (pixel size) with: WCSKEYS=cdelt1:-0.01 cdelt2:0.01 -Change the projection center with: WCSKEYS=crval1:45.2 crval2:-33.432442 Note that the axis numbers could be different. You can find the axis numbers of the longitude and latitude axis in the LOG. One can also change CD or PC elements with this keyword. CLASSIC= Make header 'classic' (remove CD/PC)? ... [N]: This keyword belongs to option 4. If your header is not a 'classic' header and you want to set a new rotation angle with keyword CROTAn, or change the values of the scaling in CDELTn, then you need to use CLASSIC=Y. If you selected LEGACY=Y then the destination header will always be 'classic'. For example in: WCSKEYS=crota2:30.3 the new value for CROTA has only effect if the header is classic. ROTATION= Classic rotation over angle (deg): ... [0.0]:" This keyword belongs to option 1. It is the angle to rotate the image, not the final value for CROTA (except when its original value was 0.0) CROTAn= Rotate image so that CROTAn becomes angle (deg): ... [0]: This keyword belongs to option 3. It gives the final angle in degrees with respect to the North and stores it in keyword CROTAn where n is the axis number of the latitude axis. LIMITS_XX= Enter pixel limits on axis XX ... [1,max]: Set the new pixel coordinate limits for axis 'XX' (e.g. LIMITS_RA=, LIMITS_DEC=). The default values are 1 which is minimum value (first pixel on that axis) and the upper limit 'max' which is given by FITS header keyword NAXISn where n represents the axis number. One can enter a value smaller than 1 and bigger than 'max' to extend the size of the image. This can be useful if we want to preserve all the data after an image rotation. FITSOUT= Enter name of sliced FITS file ... [classic_xxx] The output if this task is a FITS file. The name is entered here. A default is the input name prepended by 'classic_'. If the input was a URL, then the default file name is 'classic_fromURL.fits'. IPORDER= nint=0, bilin.=1, quad sp=2, cubic sp=3, sp4=4, sp5=5 ... [3]: Interpolation mode. The destination map is sampled on a regular grid. These positions are transformed to an equivalent position in the input map. Usually this position in not an integer position so we have to decide how to find the best pixel value. There are a number of options: Assume coordinate (48.9, 2.4). Order 0: Nearest integer: [nint(x), nint(y)] -> [49, 2] Order 1: Bilinear interpolation in square with corners [int(x), int(y)], [int(x), int(y)+1] [int(x)+1, int(y)], [int(x)+1, int(y)+1] i.e. [48,2], [48,3], [49,2], [49,3] Order 2: Quadratic spline interpolation over the 9 points derived from nint(x)-1, nint(x), nint(x)+1, and same for y. Build with piecewise polynomials f(x) = a.X^2 + b.x +c. This is a polynomial of degree 2 and order 3 Order 3: Cubic spline, polynomials degree 3, order 4 Order 4: Spline. Polynomials degree 4, order 5 Order 5: Spline. Polynomials degree 5, order 6 Bilinear, is much faster than the default Cubic Spline, but one can expect some unwanted effects like Moire patterns for image rotations. IPCVAL= Values outside boundaries of input ... [blank]: Value used for points outside the boundaries of the input. Blanks are generated for example with rotations or with destination axis sizes bigger than the original. MAKEGDS= Do you want to convert it to a GIPSY set? ... Y/[N]: One can convert to a GIPSY set after the re-projection. If FINAL=G, this prompt will not appear and RFITS is called to convert the FITS file into a GIPSY set. I one selected FINAL=F, and MAKEGDS=Y then task FITSREAD will be called for the conversion. This is to avoid some unwanted parsing that otherwise is done in RFITS. GDSNAME= Name of GDS set: ... [original without extension]: The name of the GIPSY set. The default is copied from the value in FITSOUT=. Only the .fits extension is removed. Notes: Which option should I use? ************************** * I want to change the sky system of my FITS image from galactic to equatorial and do not want a GIPSY set -> OPTION=4, MAKEGDS=N * I have a very big FITS file with spatial maps and a spectral and Stokes axis which I want to rotate the image and give it another projection system and scaling -> LEGACY=N, OPTION=4, WCSKEYS=, MAKEGDS=N * I want to import my FITS file into a GIPSY set with a fair chance that it can be used with GIPSY legacy tasks like ROTCUR and I don't mind about image skew and now that the projection system is compatible ('AIT', 'FLT', 'CYL', 'TAN', 'SIN', 'ARC', 'GLS', 'NCP', 'STG', 'MER') -> LEGACY=Y, OPTION=5 This option will not address the validity of the projection system. * I want to import my FITS file into a GIPSY set with a fair chance that it can be used with GIPSY legacy tasks like ROTCUR and I want to remove image skew and want a projection system which can be processed by GIPSY -> LEGACY=N, OPTION=4, PROJID= * I want to rotate my image so that it is aligned with the North. Skew must be removed -> OPTION=2 * I want to rotate my FITS image data but conserve the image skew -> You have to calculate a new CD matrix yourself and insert it with OPTION=4 and WCSKEYS= * I want to decrease/increase the size of the output -> All options and LIMIT_XX= keywords. -If the destination WCS comes from an second FITS file then there is no check on validity of the projection system. (Legacy) GIPSY projection systems ********************************* GIPSY compatible projection systems are coded as: AIT, FLT, CYL, TAN, SIN, ARC, GLS, NCP, STG or MER From these, the following are not recognized by GIPSY: FLT, CYL Others (GLS and NCP) have an alias in WCSLIB Translation: GIPSY codes WCSLIB translation CYL (Equivalent cylindrical projection) CEA FLT (Flat Projection) CAR GLS (Global Sinusoidal projection) SFL NCP (North Celestial Pole) SIN WCSLIB codes: CEA (Cylindrical equal area) (lambda=1, default) CAR (Plate Carree) SFL (Sanson-Flamsteed projection) SIN (Sinusoidal with two PV elements) Quality of interpolation ************************ -The quality of a re-projection can be tested with WCSFLUX. Create two plots. One with the original data and one with the re-projected data and select a shape to compare the flux. Usually you see a (small) improvement of flux conservation with higher order Spline interpolations. 'Classic' headers ***************** This GIPSY task modifies a FITS file so that it can be imported in GIPSY. It creates a so called classic header and removes PC or CD elements from the header. It checks the size of the new data set. This size cannot exceed a maximum set by GIPSY. The task does not check for the compatibility of the projection type or the compatibility of the sky system. It doesn't check the compatibility of spectral axes. What is a classic FITS header? (See also http://fits.gsfc.nasa.gov/fits_standard.html) For the transformation between pixel coordinates and world coordinates, FITS supports three conventions. First some definitions: An intermediate pixel coordinate Qi is calculated from a pixel coordinates p with: Qi = Sigma_j=1^N Mij*(Pj-Rj) Rj are the pixel coordinate elements of a reference point (FITS header item CRPIXj), j is an index for the pixel axis and i for the world axis The matrix Mij must be non-singular and its dimension is NxN where N is the number of world coordinate axes (given by FITS header item NAXIS). The conversion of Qi to intermediate world coordinate Xi is a scale: Xi = Si*Qi Formalism 1 (PC keywords) ------------------------- Formalism 1 encodes Mij in so called PCi_j keywords and scale factor Si are the values of the CDELTi keywords from the FITS header. It is obvious that the value of CDELT should not be 0.0. Formalism 2 (CD keywords) ------------------------- If the matrix and scaling are combined we get for the intermediate WORLD COORDINATE Xi: Xi = Sigma_j=1^N (Si*Mij)(Pj-Rj) FITS keywords CDi_j encodes the product Si*Mij. The units of xi are given by FITS keyword CTYPEi. Formalism 3 (Classic) ---------------------- This is the oldest but now deprecated formalism. It uses CDELTi for the scaling and CROTAn for a rotation of the image plane. n is associated with the latitude axis so often one sees CROTA2 in the header if the latitude axis is the second axis in the dataset Following the FITS standard, a number of rules is set: 1) CDELT and CROTA may co-exist with the CDi_j keywords but must be ignored if an application supports the CD formalism. 2) CROTAn must not occur with PCi_j keywords 3) CRPIXj defaults to 0.0 4) CDELT defaults to 1.0 5) CROTA defaults to 0.0 6) PCi_j defaults to 1 if i==j and to 0 otherwise. The matrix must not be singular 7) CDi_j defaults to 0.0. The matrix must not be singular. 8) CDi_j and PCi_j must not appear together in a header. Alternate WCS axis descriptions -------------------------------- A World Coordinate System (WCS) can be described by an alternative set of keywords. For this keywords a character in the range [A..Z] is appended. In our software we made the assumption that the primary description contains all the necessary information to derive a 'classic' header and therefore we will ignore alternate header descriptions. Conversion to a formalism 3 ('classic') header ---------------------------------------------- Many FITS readers from the past are not upgraded to process FITS files with headers written using formalism 1 or 2. The purpose of this application is to convert a FITS file to a file that can be read and interpreted by old FITS readers. For GIPSY we require FITS headers to be written using formalism 3. If keywords are missing, they will be derived and a comment will be inserted about the keyword not being an original keyword. The method that converts the header, tries to process it first with WCSLIB (tools to interpret the world coordinate system as described in a FITS header). If this fails, then we are sure that the header is incorrect and we cannot proceed. One should use a FITS tool like 'fv' (the Interactive FITS File Editor from Nasa) to repair the header. The conversion process starts with exploring the spatial part of the header. It knows which axes are spatial and it reads the corresponding keywords (CDELTi, CROTAn, CDi_j, PCi_j and PC00i_00j (old format for PC elements). If there is no CD or PC matrix, then the conversion routine returns the unaltered original header. If it finds a PC matrix and no CD matrix then the header should contain CDELT keywords. With the values of these keywords we create a CD matrix: |cd11 cd12| = |cdelt1 0| * |pc11 pc12| |cd21 cd22| |0 cdelt2| |pc21 pc22| Note 1: We replaced notation i_j by ij so cd11 == CD1_1 Note 2: For the moment we restricted the problem to the 2 dim. spatial case because that is what we need to retrieve a value for CROTA, the rotation of the image.) Note 3: We assumed that the PC matrix did not represent transposed axes as in: | 0 1 0 | PC = | 0 0 1 | | 1 0 0 | If cd12 == 0.0 and cd12 == 0.0 then CROTA is obviously 0. There is no rotation and CDELT1 = cd11, CDELT2 = cd22 If one or both values of cd12, cd21 is not zero then we expect a value for CROTA unequal to zero and/or skew. We calculate the scaling parameters CDELT with: CDELT1 = sqrt(cd11*cd11+cd21*cd21) CDELT2 = sqrt(cd12*cd12+cd22*cd22) The determinant of the matrix is: det = cd11*cd22 - cd12*cd21 This value cannot be 0 because we required that the matrix is non-singular. Further we distinguish two situations: a determinant < 0 and a determinant > 0 (zero was already excluded). Then we derive two rotations. If these are equal, the image is not skewed. If they are not equal, we derive the rotation from the average of the two calculated angles. As a measure of skew, we calculated the difference between the two rotation angles. Here is a piece of the actual code: sign = 1.0 if det < 0.0: cdeltlon_cd = -cdeltlon_cd sign = -1.0 rot1_cd = atan2(-cd21, sign*cd11) rot2_cd = atan2(sign*cd12, cd22) rot_av = (rot1_cd+rot2_cd)/2.0 crota_cd = degrees(rot_av) skew = degrees(abs(rot1_cd-rot2_cd)) New values of CDELT and CROTA will be inserted in the new 'classic' header only if they were not part of the original header. The process continues with non spatial axes. For these axes we cannot derive a rotation. We only need to find a CDELTi if for axis i no value could be found in the header. If this value cannot be derived from the a CD matrix (usually with diagonal element CDi_i) then the default 1.0 is assumed. Note that there will be a warning about this in the comment string for the corresponding keyword in the new 'classic' FITS header. Finally there is some cleaning up. First all CD/PC elements are removed for all the axes in the data set. Second, some unwanted keywords are removed. The current list is: ["XTENSION", "EXTNAME", "EXTEND"] Data limits *********** In GIPSY we can attach meta data to a data structure of slices of a data structure. Even pixels can have their own header with keywords. To be able to identify such headers, a so called coordinate word is used. This coordinate word is currently a 32 bits integer. This implies that there is a limit to the size of a data cube in GIPSY. In practice this limit is 2^31-1 pixels. It does not matter what these pixels represent, i.e., the value of keyword BITPIX= in the FITS header is not important. So for the input data structure we prompt for axis limits. This enables a user to decrease the size of original data structure to meet the GIPSY requirements. Updates: Feb 26, 2011: VOG document created.