Program: REPROJ
Purpose: Reproject 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 reproject:
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= Give input 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. centre for output: [copy from input map]
Specify for the output map the position of the
projection centre (PC), i.e. the position of grid (0,0).
This position can be specified either in grid coordinates
of the input map or in physical 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]
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 1950.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. The result is is obtained much faster
but less accurate. However for rotations and shifts this
is completely negligible. For rotations the default
values for the matrix are the maximum values (max.
acceleration).
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 reprojects a spatial input set into another
set with a different coordinate system. A coordinate
system is specified by the following parameters.
-projection centre,
-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 centre:
==================
The projection centre (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).
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 (Equatorial map, 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 centre of the input map
as the new projection centre. This can be realized with
'AC' (axis centre 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.
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 centre is always
this input centre.
For a rotation around a different centre, 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 centre.
A rotation can then safely be performed around the
default rotation centre 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
bilinear 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.
Reprojecting 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 rotation:
=================
i, j are grids in a map. Then the physical values are
calculated using:
x = (i-CRPIX_i) * CDELT_i
y = (j-CRPIX_j) * CDELT_j
(CRPIX is the reference pixel, if CRPIX=1 then the first
pixel is grid 0. CDELT is the grid spacing)
If rho is equal to CROTA (angle associated with the
latitude-like axis) and l, m are the direction cosines,
then:
l = x.cos(rho) - y.sin(rho)
m = y.cos(rho) + x.sin(rho)
where l is the longitude and m the latitude direction.
Rho is an angle from the +y axis to the +x axis, i.e.
counter clockwise if the grid spacing in longitude is
less than 0.
Example:
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
library interpol.
Aug 25, 2000: VOG, Added DIMINISH= keyword to decrease
output dimensionality for one input
subset from a > 2 dim cube.