Program: REPROJ Purpose: Re-project a spatial map into another map with a different coordinate system. Category: COORDINATES, HEADER, MANIPULATION File: reproj.c Author: M. Vogelaar Keywords: INSET= Give set (subset(s)) with data to re-project: Input set (and subsets) which contain sky maps which should be reprojected. Maximum number of subsets is 2048. BOX= Area which should be reprojected: [entire structure] The input area that you select here is the part that will be reprojected. Data outside this box will not be transferred. DEFSET= Copy system from reference set, subsets: [manual input] If specified, the defaults for coordinate parameters OUTPOS=, CDELT=, CROTA=, EPOCH=, OUTBOX=, SKYSYS=, and PROSYS= will be taken from this set. ROTATEONLY= Rotation only (i.e.fixed sky-&proj.systems)? Y/[N] If you only want to rotate a map and do not change the sky and projection systems (ROTATEONLY=Y) then the keywords SKYSYS=, EPOCH=, CROTA= and PROSYS= are skipped. OUTPOS= Give proj. center for output: [copy from input map] Specify for the output map the position of the projection center (PC), i.e. a position somewhere in grid (0,0). The exact position (given by OUTPOS='PC') is derived from the fractional parts of the values of CRPIXn in the header in the set. This position can be specified either in grid coordinates or in world coordinates. The PC is the intersection of the line of sight with the celestial sphere and as such relates the projection system to the grid mesh. A change in PC thus implies a repro- jection of the projection system onto the grid and only equals a shift in pixels if both, the input and output systems are flat (axis name extension FLT). Note that the order of input is always x, y, i.e. longitude- latitude (e.g. for a RA-DEC map) or latitude-longitude (e.g. for a DEC-RA map) depending on the header values of 'CTYPE'. CDELT= New grid spacings (Dx,Dy) in ARCSEC: [copy from input map] Regrid the input to these new grid spacings. Note that the grid spacings are always in seconds of arc, NOT in the units of your header (CUNIT). CHANGE= Do you want to change sign of cdelt? [Y]/N If you entered a positive value for the grid spacing in longitude in an equatorial system, then CHANGE=Y will change its sign. ROTANGLE= Rotate map over .... degrees: [0.0] Specify an angle in degrees over which you want to rotate the input. The direction is counter clockwise for systems with a negative grid spacing in longitude. Note that this value will be added to the value of 'CROTA' (see description at CROTA=) as found in the header of the input set (see description!). If ROTATEONLY=N then also the keywords SKYSYS=, EPOCH= and PROSYS are asked. SKYSYS= Output sky system: [copy from input map] Change a map in a way that the x- and y axes of the image correspond to the principal axes of the new sky system. The following sky systems are implemented: Skysys CTYPE_l CTYPE_m Meaning =================================================== 1 RA DEC Equatorial (EPOCH 1950.0) 2 GLON GLAT Galactic 3 ELON ELAT Ecliptic (EPOCH 1950.0) 4 SLON SLAT Super galactic 5 RA DEC Equatorial (EPOCH 2000.0) EPOCH= Give epoch of new sky system: [copy from input map] For equatorial and ecliptic input sky systems it is possible to change the epoch. The default epoch is copied from the input set. If the input set has no epoch keyword in its header, then 2000.0 is assumed. If the sky system in the header is number 1 (see table above) and the epoch is 2000.0, then REPROJ changes this number automatically into 5. PROSYS= Output projection system: [copy from input map] Change the projection system. Default is the projection system of the input map. Projection systems can be identified in the axis names by the following postfixes: PROSYS postfix Meaning ================================================= 1 AIT Aitoff Equal Area projection 2 CYL Equivalent Cylindrical projection 3 FLT flat projection 4 TAN Gnomonic projection 5 SIN Orthographic projection 6 ARC Rectangular projection 7 GLS Transversal projection 8 NCP North Celestial Pole projection 9 STG Stereographic projection 10 MER Mercator projection CROTA= New value for map rotation (deg.): [copy from input map] Header item CROTA stores the map rotation. CROTA is the angle in degrees between the +y axis and the +m axis (latitude axis) at the position of the PC. The angle is counter-clockwise if the grid separation (CDELT) in longitude is < 0.0. If you want to rotate AND change the sky- and projection system as well, use the CROTA= keyword, else use ROTATEONLY=Y and ROTANGLE= (See description for additional information). OUTBOX= New box (in grids): [minimum size] The box of the input set is transformed using the given transformation parameters (OUTPOS=, CDELT=, CROTA= SKYSYS= and PROSYS=) to a box in the new system. This is the output box that just contains the entire input image after transformation and is therefore also the default box. DIMINISH= Decrease output dimensionality to 2? [Y]/N If you entered ONE 2-dim subset from a set with 3 or more axes, then the default output dimension is not copied from the input set, but is reduced to two. This avoids huge and almost empty sets to be created. Lost axes become so called hidden axis. OUTSET= Give output set, subset(s): This will be the destination of the reprojected set. A name is sufficient. SPEEDMAT= Size of 'position' interpolation box: [1,1] Accelerate the calculations by interpolating positions instead of transforming them using formulas for coordinate transformations. For interpolation boxes with a size bigger than 1 x 1 the result is obtained faster but less accurate. DATAMODE= Data acquisition (B)ilinear/(N)earest pix.: [B]/N Obtain output pixel values by (B)ilinear interpolation in a 2 x 2 matrix or get value of the (N)earest pixel. See description at: Interpolation of data (not positions). Description: This program re-projects a spatial input set into another set with a different coordinate system. A coordinate system is specified by the following parameters. -projection center, -grid spacing, -rotation angle, -sky system and epoch. -projection system New values for these parameters are either given by the user or read from the header of the set in DEFSET= Projection center: ================== The projection center (PC) is the intersection point of the line of sight with the celestial sphere. The PC in physical coordinates is attached to grid coordinate (0,0) but not necessarily in the center of that grid. Together with the sky & projection systems, the grid- spacing and the map rotation, it connects a physical coordinate system to a grid mesh. A PC can be specified in the standard GIPSY way: For spatial axes there are a number of prefixes: * ; for RA or DEC in resp. HMS and DMS. *1950 ; for RA or DEC in resp. HMS and DMS in EPOCH 1950.0 *xxxx.x ; for RA or DEC in resp. HMS and DMS in EPOCH xxxx.x G ; Galactic longitude or latitude in degrees E ; Ecliptic longitude or latitude in degrees S ; Supergalactic longitude or latitude in degrees The prefixes must be repeated for both directions. There must be a space between prefix and coordinate specification. Examples (e.g. an Equatorial map with longitude, latitude): OUTPOS=10 5 RA at grid 10 of the input map, DEC at grid 5 OUTPOS=* 10 12 8 * -67 8 9.6 RA = 10 hour, 12 min, 8 sec, DEC = -67 deg, 8 min, 9.6 sec., in a 2-d area and in the epoch as found in the header of the set: OUTPOS=*2000.0 3 14 38.02 *2000.0 41 13 54.8 Input of RA 3 h 14 m 38.02 s, DEC 41 d 13 m 54.84 s in epoch 2000.0: It is also possible to specify the grid center of the input map as the new projection center. This can be realized with 'AC' (axis center e.g. OUTPOS=AC). If the length of axis i is NAXISi and the reference pixel is CRPIXi, then the i-th coordinate is given by the expression NAXISi / 2 - CRPIXi. A projection system is given in the input system. If you transform to another sky system, then this projection center is transformed to the new system before it is inserted in the output header. Now you are sure that the projection center is the same physical location in the sky. Grid spacing: ============= REPROJ converts only spatial maps. For each axis in a map it must be possible to convert a spacing in header units to seconds of arc. The sign of the grid spacing in longitude determines the direction of rotation. For equatorial maps with a negative grid spacing in longitude (RA), rotation will be counter-clockwise. Rotation angle: =============== For rotations of maps where sky and projection systems are fixed, use ROTATEONLY=YES. Then the angle over which you want to rotate will be asked in ROTANGLE=. Otherwise the program will prompt you to give the new systems and the "header" rotation angle 'CROTA'. In both cases a rotation over x degrees implies a change in 'CROTA' with +/- x degrees. 'CROTA' is defined for the input PC and thus the rotation center is always this input center. For a rotation around a different center, program REPROJ has to run twice. The first time to change the PC, the second time to rotate the map. However, for maps which do not cover a too large fraction of the sky, and/or are not too close to a pole of the projection system, the result of a rotation is nearly independent of your choice of the rotation center. A rotation can then safely be performed around the default rotation center which is the original PC. Blanks: ======= If a pixel at certain position in the input set is a blank, then the corresponding pixel in the output set will be set to blank. If a pixel in the output set corresponds to a pixel outside the box or frame of the input set, this pixel is set to blank also. Interpolation of data (not positions): ====================================== For each grid position (integer x, y) in the OUTput set a grid position (floating x', y') in the INput set is calculated. The position (x', y') has distance dx, dy to the nearest pixel m1. If dx and dy were both 0.0 then the pixel value of m1 is returned. However, in most cases the values for dx and dy will not be equal to 0.0. Then the 3 closest neighbours are involved in an interpolation. If the pixel value at (x', y') is a blank, then a blank is returned to the output set. Else, an interpolation is used to calculate the pixel value. Given 4 values m1,m2,m3,m4 at positions: m4 m3 (0,1) o--------o (1,1) | | | | dy ^ | | dx | (0,0) o-->-----o (1,0) m1 m2 and two fractions: 0.0 <= dx <= 0.5 and 0.0 <= dy <= 0.5 then there are a number of interpolation schemes depending on the number of blanks pixels. If all pixels are non blank, a bi-linear interpolation is involved. If the start pixel m1 is blank then a blank is returned. If one of the other pixels is blank, the interpolation is in the plane of the (3) pixels that are not blank. For two non blank pixels there is a linear interpolation. The number of blanks in a map will therefore not increase (blot) or decrease (patch). If you change the values for the grid spacing, then only intensities will be conserved (not the flux!). Interpolation of positions (the 'speed matrix'): ================================================ For a 1000x1000 pixels map, a rotation over 45 degrees takes 63.3 cpu sec on a certain machine if all pixels were transformed. But if a 'speedmatrix' of 1000x1000 is used then the process takes 14.4 cpu sec. The less linear the projection is, the less accurate the re-projection. Re-projecting a map without a coordinate system: =============================================== It is not possible to shift or rotate an arbitrary map without a valid coordinate system. You can fit such a coordinate system if you know the physical position of some grids in your map by using ASTROM. A coordinate system of an existing set can be changed by program FIXHED. Appendix: Function transform() ============================== Here we outline the procedure of transforming one system to another. The core is function transform(). Four positions (each two coordinates) of the input box are transformed to the corresponding world coordinates. These are transformed back to grids for the new system. From this new set of grids one can compose a new box for the output set. Then for each grid in the output set the corresponding grid in the input set is calculated via intermediate world coordinates. These grids in the input set usually are not integers so an interpolation using neighbours is applied (or not if the user wants the nearest pixel only). This way we sample the output set on its equally spaced grids. An small but important correction is made for non-integer header values of CRPIX. These fractional parts are the cause of the fact that the projection center needs not to be at the center of grid (0,0). These offsets are processed in function transform(). Example: We want to extract one channel map from a GIPSY data cube and rotate that map 30 degrees (towards the East). The rotation is around the projection center of the input map. BOX= CDELT= DATAMODE= DEFSET= DIMINISH= y INSET= aurora velo 50 OUTBOX= OUTPOS= OUTSET= aurora_rot30 ROTANGLE= 30 ROTATEONLY= y SPEEDMAT= 1 1 Updates: Jul 10, 1990: VOG, Document created Apr 10, 1992: VOG, GDSPOS included for IN/OUTPOS= May 17, 1994: VOG, Changed integer position conversions Oct 4, 1994: VOG, Changed NINT in interpolation positions to INT and subtracted 1 for negative pixel positions. Sep 21, 1995: VOG, Rewritten in C Oct 31, 1995: VOG, Added 'getnewaxisname' function Jul 27, 2000: VOG, Changed local 'getipval' to external function interpol(). Aug 25, 2000: VOG, Added DIMINISH= keyword to decrease output dimensionality for one input subset from a > 2 dim cube. Apr 12, 2009: VOG,-Changed macro NINT to version with floor() -Grids corrected for offsets due to non integer values of CRPIX. -Changed default value of epoch from 1950 to 2000 -Added documentation to core function transform(). Jul 17, 2009: VOG, An extra transformation was added to transform CRVAL's (for the projection center) between two different sky systems. Oct 21, 2009: VOG, Removed bug in nearest integer data acquisition M[x][y] -> M[y][x] Mar 24, 2013: VOG, Output axes names lost hyphens in CTYPE. Repaired