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