WCS.jl

License

Astronomical World Coordinate System library for Julia. This package wraps the WCSLIB C library. This is a tool made for users who are already familiar with WCS transformations. If you are not, please reference the following manuscripts

Installation

From the REPL, press ] to enter Pkg mode

pkg> add WCS

Usage

Import the library:

There are many ways to utilize WCS transformations. Let's make one for a 2-dimensional array (like an image) from scratch:

julia> wcs = WCSTransform(2;
           cdelt = [-0.066667, 0.066667],
           ctype = ["RA---AIR", "DEC--AIR"],
           crpix = [-234.75, 8.3393],
           crval = [0., -90],
           pv    = [(2, 1, 45.0)],
       )
WCSTransform(naxis=2, cdelt=[-0.066667, 0.066667], crval=[0.0, -90.0], crpix=[-234.75, 8.3393])

We can also create one directly from a header loaded from FITSIO.jl if it contains the appropriate keywords:

julia> using FITSIO: read_header

julia> using Downloads: download

julia> header = (read_header ∘ download)("https://fits.gsfc.nasa.gov/nrao_data/samples/wcs/m82rad.fits");

julia> WCS.from_header(header; ignore_rejected = true)
1-element Vector{WCSTransform}:
 WCSTransform(naxis=4, cdelt=[-0.0002247777738, 0.0002273333375, 5.0e7, 1.0], crval=[147.927100003, 69.9079999924, 1.49e9, 0.0], crpix=[0.0, 0.0, 0.0, 1.0])

header file

julia> typeof(header)
FITSIO.FITSHeader

julia> header
SIMPLE  =                    T
BITPIX  =                   16
NAXIS   =                    4
NAXIS1  =                  558
NAXIS2  =                  262
NAXIS3  =                    1
NAXIS4  =                    1
EXTEND  =                    T / Tables following main image
BLOCKED =                    T / Tape may be blocked
OBJECT  = 'U05322  '
TELESCOP= 'VLA     '
INSTRUME= 'VLA     '
OBSERVER= 'AC243   '
DATE-OBS= '06/02/89'
DATE-MAP= '13/09/94'
BSCALE  =     1.95069823037e-6 / REAL = TAPE * BSCALE + BZERO
BZERO   =      0.0548100471497
BUNIT   = 'JY/BEAM '
EPOCH   =               1950.0 / EPOCH OF RA DEC
BLANK   =               -32768 / Blanked pixel tape value
VELREF  =                    3 / >256 RADIO, 1 LSR 2 HEL 3 OBS
ALTRVAL =             100000.0 / ALTENATE FREQ/VEL REF VALUE
ALTRPIX =                  1.0 / ALTENATE FREQ/VEL REF PIXEL
YSHIFT  =      1.421085472e-14 / NET SHIFT OF PHASE CENTER INY
DATAMAX =         0.1187149212 / MAX PIXEL VALUE
DATAMIN =      -0.009094826877 / MIN PIXEL VALUE
CTYPE1  = 'RA---ARC'
CRVAL1  =        147.927100003
CDELT1  =     -0.0002247777738
CRPIX1  =                301.5
CROTA1  =                  0.0
CTYPE2  = 'DEC--ARC'
CRVAL2  =        69.9079999924
CDELT2  =      0.0002273333375
CRPIX2  =                 95.5
CROTA2  =          24.26000023
CTYPE3  = 'FREQ    '
CRVAL3  =               1.49e9
CDELT3  =                5.0e7
CRPIX3  =                  1.0
CROTA3  =                  0.0
CTYPE4  = 'STOKES  '
CRVAL4  =                  1.0
CDELT4  =                  1.0
CRPIX4  =                  1.0
CROTA4  =                  0.0
HISTORY --------------------------------------------------------------------
HISTORY /Begin "HISTORY" information found in fits tape header by IMLOD
HISTORY BLOCKED =                    T /Tape may be blocked
HISTORY --------------------------------------------------------------------
HISTORY /Begin "HISTORY" information found in fits tape header by IMLOD
HISTORY BLOCKED =                    T /Tape may be blocked
HISTORY AIPS   IMNAME='NGC3034     ' IMCLASS='SUBIM ' IMSEQ=   1     /
HISTORY AIPS   USERNO=  135            /
HISTORY AIPS   CLEAN BMAJ=  5.0000E-04 BMIN=  5.0000E-04 BPA=   0.00
HISTORY AIPS   CLEAN NITER=     2000 PRODUCT=1   / NORMAL
HISTORY /END FITS tape header "HISTORY" information
HISTORY --------------------------------------------------------------------
HISTORY IMLOD OUTNAME ='
HISTORY IMLOD OUTSEQ =    0     INTAPE = 1     OUTDISK= 2
HISTORY IMLOD INFILE = 'FITS2:NGC3034.A                                 '
HISTORY IMLOD RELEASE = '15JUL94'
HISTORY SUBIM RELEASE ='15JUL94 '  /********* Start 25-MAR-1994   09:37:01
HISTORY SUBIM INNAME='M82 in radio'   INCLASS='SUBIM '
HISTORY SUBIM INSEQ=    1   INDISK=  2
HISTORY SUBIM INTYPE ='MA'      USERID=  103
HISTORY SUBIM OUTNAME='M82 in radio'   OUTCLASS='SUBIM '
HISTORY SUBIM OUTSEQ=    2   OUTDISK=  2
HISTORY SUBIM BLC    =    1,    1,    1,    1,    1,    1,    1
HISTORY SUBIM TRC    =  256,  256,    1,    1,    1,    1,    1
HISTORY SUBIM XINC   =    2      YINC =    2
HISTORY SUBIM OPCODE = 'AVE '
HISTORY HGEOM RELEASE ='15JUL94 '  /********* Start 25-MAR-1994   09:37:24
HISTORY HGEOM INNAME='M82 in radio'   INCLASS='SUBIM '
HISTORY HGEOM INSEQ=    2   INDISK=  2
HISTORY HGEOM IN2NAME='M82 in blue '   IN2CLASS='SUBIM '
HISTORY HGEOM IN2SEQ=    1   IN2DISK=  2
HISTORY HGEOM OUTNAME='M82 in radio'   OUTCLASS='      '
HISTORY HGEOM OUTSEQ=    1   OUTDISK=  2
HISTORY HGEOM BLC =    1,    1,   1,   1,   1,   1,   1
HISTORY HGEOM TRC =  128,  128,   1,   1,   1,   1,   1
HISTORY HGEOM IMSIZE =  558,  263          / Output image size
HISTORY HGEOM  / Interpolation order used was BiCubic
HISTORY HGEOM  / Indeterminate pixels filled with magic values
HISTORY HGEOM  /     131164 Pixels blanked due to memory limits or geometry
HISTORY HGEOM  /          0 Pixels blanked due to input blanked pixels
HISTORY AIPS   IMNAME='M82 in radio' IMCLASS='HGEOM ' IMSEQ=   1     /
HISTORY AIPS   USERNO=  103            /
HISTORY AIPS   CLEAN BMAJ=  5.0000E-04 BMIN=  5.0000E-04 BPA=   0.00
HISTORY AIPS   CLEAN NITER=     2000 PRODUCT=1   / NORMAL
HISTORY /END FITS tape header "HISTORY" information
HISTORY --------------------------------------------------------------------
HISTORY IMLOD OUTNAME ='
HISTORY IMLOD OUTSEQ =    0     INTAPE = 1     OUTDISK= 4
HISTORY IMLOD INFILE = 'FITS2:M82RADIO.FITS                             '
HISTORY IMLOD RELEASE = '15JAN95'
HISTORY SUBIM RELEASE ='15JAN95 '  /********* Start 13-SEP-1994   13:00:09
HISTORY SUBIM INNAME='M82 in radio'   INCLASS='HGEOM '
HISTORY SUBIM INSEQ=    1   INDISK=  4
HISTORY SUBIM INTYPE ='MA'      USERID=  103
HISTORY SUBIM OUTNAME='M82         '   OUTCLASS='Radio '
HISTORY SUBIM OUTSEQ=    3   OUTDISK=  4
HISTORY SUBIM BLC    =    1,    1,    1,    1,    1,    1,    1
HISTORY SUBIM TRC    =  558,  262,    1,    1,    1,    1,    1
HISTORY SUBIM XINC   =    1      YINC =    1
HISTORY SUBIM OPCODE = '    '
ORIGIN  = 'AIPSGorilla (NRAOCV IPX) 15JAN95'
DATE    = '13/09/94'           / File written on dd/mm/yy
HISTORY AIPS   IMNAME='M82         ' IMCLASS='Radio ' IMSEQ=   3     /
HISTORY AIPS   USERNO=  103            /
HISTORY AIPS   CLEAN BMAJ=  5.0000E-04 BMIN=  5.0000E-04 BPA=   0.00
HISTORY AIPS   CLEAN NITER=     2000 PRODUCT=1   / NORMAL

Now we can do conversions between pixel and world coordinates.

Note

WCS transforms use pixel (0, 0) as the top-left corner (not 1-indexed!!)

julia> pixcoords = [0.0  24.0  45.0;  # x coordinates
                    0.0  38.0  98.0]  # y coordinates
2×3 Matrix{Float64}:
 0.0  24.0  45.0
 0.0  38.0  98.0

julia> worldcoords = pix_to_world(wcs, pixcoords)
2×3 Matrix{Float64}:
 267.965   276.539   287.771
 -73.7366  -71.9741  -69.6781

julia> pixcoords = world_to_pix(wcs, worldcoords)
2×3 Matrix{Float64}:
  1.16529e-12  24.0  45.0
 -7.10543e-14  38.0  98.0

API/Reference

WCS.WCSTransform — Type

WCSTransform(naxis; kwds...)

Construct a WCS transformation with the given number of axes naxis. Keyword arguments can be passed to set various attributes of the transform. Specifying keyword arguments is equivalent to setting them after construction:

julia> wcs = WCSTransform(2; crpix=[1000., 1000.])
WCSTransform(naxis=2, cdelt=[1.0, 1.0], crval=[0.0, 0.0], crpix=[1000.0, 1000.0])

is equivalent to:

julia> wcs = WCSTransform(2)
WCSTransform(naxis=2, cdelt=[1.0, 1.0], crval=[0.0, 0.0], crpix=[0.0, 0.0])

julia> wcs.crpix = [1000., 1000.];

julia> wcs
WCSTransform(naxis=2, cdelt=[1.0, 1.0], crval=[0.0, 0.0], crpix=[1000.0, 1000.0])

Properties

Below is the entire list of public properties for a WCSTransform

Keyword	Type	Description
naxis	`Int32`	Number of dimensions
crval	`Vector{Float64}` (length `naxis`)	coordinate value at reference point
crpix	`Vector{Float64}` (length `naxis`)	array location of the reference point in pixels
cdelt	`Vector{Float64}` (length `naxis`)	coordinate increment at reference point
crder	`Vector{Float64}` (length `naxis`)	random error in coordinate
csyer	`Vector{Float64}` (length `naxis`)	systematic error in coordinate
ctype	`Vector{String}` (length `naxis`)	axis type (8 characters)
crota	`Vector{Float64}` (length `naxis`)	rotation from stated coordinate type
cunit	`Vector{String}` (length `naxis`)	units of axes
cunit	`Vector{String}` (length `naxis`)	names of axes
pc	`Matrix{Float64}` (size `naxis` × `naxis`)	linear transformation matrix
cd	`Matrix{Float64}` (size `naxis` × `naxis`)	linear transformation matrix (with scale)
equinox	`Float64`	the equinox associated with dynamical equatorial or ecliptic coordinate systems
latpole	`Float64`	the native latitude of the celestial pole
lonpole	`Float64`	the native longitude of the celestial pole
mjdavg	`Float64`	Modified Julian Date corresponding to `DATE-AVG`
mjdobs	`Float64`	Modified Julian Date corresponding to `DATE-OBS`
restfrq	`Float64`	rest frequency (Hz)
restwav	`Float64`	rest wavelength (m)
velangl	`Float64`	velocity angle
velosys	`Float64`	relative radial velocity
zsource	`Float64`	the redshift of the source
colnum	`Int32`	column of FITS binary table associated with this WCS
dateavg	`String`	representative mid-point of the date of observation
dateobs	`String`	start of the date of observation
radesys	`String`	the equatorial or ecliptic coordinate system type
specsys	`String`	spectral reference frame (standard of rest)
ssysobs	`String`	spectral reference frame
ssyssrc	`String`	spectral reference frame for redshift
wcsname	`String`	name of this coordinate representation
obsgeo	`Vector{Float64}` (length 3 or 6)	location of the observer in a standard terrestrial reference frame
alt	`String`	character code for alternate coordinate descriptions

source

WCS.from_header — Method

from_header(header[; relax=WCS.HDR_ALL, ctrl=0, ignore_rejected=false, table=false])

Parse the FITS image header in the String header, returning a Vector{WCSTransform} giving all the transforms defined in the header. The relax determines the treatment of non-standard keywords. The default is to accept all known non-standard keywords. Use relax=WCS.HDR_NONE to ignore all non-standard keywords. Use, e.g., relax=(WCS.HDR_RADECSYS & WCS.HDR_CROTAia) to only accept selected non-standard keywords.

source

WCS.obsfix — Method

obsfix(ctrl::Integer, wcs::WCSTransform)

Complete the obsgeo field wcs of observatory coordinates. That is, if only the (x,y,z) Cartesian coordinate triplet or the (l,b,h) geodetic coordinate triplet are set, then it derives the other triplet from it. If both triplets are set, then it checks for consistency at the level of 1 metre.

Parameters

ctrl: flag that controls behaviour if one triplet is defined and the other is only partially defined:
- 0: Reset only the undefined elements of an incomplete coordinate triplet.
- 1: Reset all elements of an incomplete triplet.
- 2: Don't make any changes, check for consistency only. Returns an error if either of the two triplets is incomplete.
wcs: Coordinate transformation parameters. Its obsgeo field may be changed.

Returns

-1: No change required (not an error).
0: Success.
1: Null wcsprm pointer passed.
5: Invalid parameter value.

source

WCS.pix_to_world! — Method

pix_to_world!(wcs, pixcoords, worldcoords[; stat, imcoords, phi, theta])

Convert the array of pixel coordinates pixcoords to world coordinates according to the WCSTransform wcs, storing the result in the worldcoords and stat arrays. pixcoords should be a 2-d array where pixcoords[:, i] is the i-th set of coordinates, or a 1-d array representing a single set of coordinates. worldcoords must be the same size and type as pixcoords.

If given, the arrays stat, imcoords, phi, theta will be used to store intermediate results. Their sizes and types must all match pixcoords, except for stat which should be the same size but of type Cint (typically Int32).

source

WCS.pix_to_world — Method

pix_to_world(wcs, pixcoords)

Convert the array of pixel coordinates pixcoords to world coordinates according to the WCSTransform wcs. pixcoords should be a 2-d array where pixcoords[:, i] is the i-th set of coordinates, or a 1-d array representing a single set of coordinates.

The return value is the same shape as pixcoords.

source

WCS.to_header — Method

to_header(wcs[; relax=WCS.HDR_NONE])

Encode the WCSTransform wcs as a FITS header string. The relax keyword controls how non-standard extensions to the WCS standard are handled.

source

WCS.world_to_pix! — Method

world_to_pix!(wcs, worldcoords, pixcoords[; stat, phi, theta, imcoords])

Convert the array of pixel coordinates worldcoords to pixel coordinates according to the WCSTransform wcs, storing the result in the pixcoords array. worldcoords should be a 2-d array where worldcoords[:, i] is the i-th set of coordinates, or a 1-d array representing a single set of coordinates. pixcoords must be the same size and type as worldcoords.

If given, the arrays stat, imcoords, phi, theta will be used to store intermediate results. Their sizes and types must all match worldcoords, except for stat which should be the same size but of type Cint (typically Int32).

source

WCS.world_to_pix — Method

world_to_pix(wcs, worldcoords)

Convert the array of world coordinates worldcoords to pixel coordinates according to the WCSTransform wcs. worldcoords is a 2-d array where worldcoords[:, i] is the i-th set of coordinates, or a 1-d array representing a single set of coordinates.

The return value is the same size as worldcoords.

source