Why WCS.jl package

After the introduction of package FITSIO.jl for handling FITS files in Julia, another package was required to manage the World Coordinate System(WCS) transformations of a FITS object. The WCS info. is present in the header of one of the HDU(s) of a FITS file. These transformations map the pixel locations in an image to their real-world units, such as their position on the sky sphere. These transformations can work both forward (from pixel to sky) and backward (from sky to pixel). This package wraps the WCSLIB C library.

Installation

WCS.jl is a registered package in Julia and can be installed from package manager as: pkg> add WCS

or it can be build from source as: julia> Pkg.clone("https://github.com/JuliaAstro/WCS.jl")

Usage

After installing the package, you can start using it with julia> using WCS

Creating a WCS object

WCS objects are defined by a type WCSTransform which contains the following fields based on the information present in the header:

    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    crval::Ptr{Cdouble}
    cunit::Ptr{Cvoid}
    ctype::Ptr{Cvoid}
    lonpole::Cdouble
    latpole::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    npv::Cint
    npvmax::Cint
    pv::Ptr{PVCard}
    nps::Cint
    npsmax::Cint
    ps::Ptr{PSCard}
    cd::Ptr{Cdouble}
    crota::Ptr{Cdouble}
    altlin::Cint
    velref::Cint
    alt::NTuple{4, UInt8}
    colnum::Cint
    colax::Ptr{Cint}
    cname::Ptr{Cvoid}
    crder::Ptr{Cdouble}
    csyer::Ptr{Cdouble}
    czphs::Ptr{Cdouble}
    cperi::Ptr{Cdouble}
    wcsname::NTuple{72, UInt8}
    timesys::NTuple{72, UInt8}
    trefpos::NTuple{72, UInt8}
    trefdir::NTuple{72, UInt8}
    plephem::NTuple{72, UInt8}
    timeunit::NTuple{72, UInt8}
    dateref::NTuple{72, UInt8}
    mjdref::NTuple{2, Cdouble}
    timeoffs::Cdouble
    dateobs::NTuple{72, UInt8}
    datebeg::NTuple{72, UInt8}
    dateavg::NTuple{72, UInt8}
    dateend::NTuple{72, UInt8}
    mjdobs::Cdouble
    mjdbeg::Cdouble
    mjdavg::Cdouble
    mjdend::Cdouble
    jepoch::Cdouble
    bepoch::Cdouble
    tstart::Cdouble
    tstop::Cdouble
    xposure::Cdouble
    telapse::Cdouble
    timsyer::Cdouble
    timrder::Cdouble
    timedel::Cdouble
    timepixr::Cdouble
    obsgeo::NTuple{6, Cdouble}
    obsorbit::NTuple{72, UInt8}
    radesys::NTuple{72, UInt8}
    equinox::Cdouble
    specsys::NTuple{72, UInt8}
    ssysobs::NTuple{72, UInt8}
    velosys::Cdouble
    zsource::Cdouble
    ssyssrc::NTuple{72, UInt8}
    velangl::Cdouble
    ntab::Cint
    nwtb::Cint
    tab::Ptr{Cvoid}  # Ptr{tabprm}
    wtb::Ptr{Cvoid}  # Ptr{wtbarr}
    lngtyp::NTuple{8, UInt8}
    lattyp::NTuple{8, UInt8}
    lng::Cint
    lat::Cint
    spec::Cint
    cubeface::Cint
    types::Ptr{Cint}
    lin::linprm
    cel::celprm
    spc::spcprm
    err::Ptr{WCSErr}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_crval::Ptr{Cdouble}
    m_cunit::Ptr{Cvoid}
    m_ctype::Ptr{Cvoid}
    m_pv::Ptr{PVCard}
    m_ps::Ptr{PSCard}
    m_cd::Ptr{Cdouble}
    m_crota::Ptr{Cdouble}
    m_colax::Ptr{Cint}
    m_cname::Ptr{Cvoid}
    m_crder::Ptr{Cdouble}
    m_csyer::Ptr{Cdouble}
    m_czphs::Ptr{Cdouble}
    m_cperi::Ptr{Cdouble}
    m_tab::Ptr{Cvoid}  # Ptr{tabprm}
    m_wtb::Ptr{Cvoid}  # Ptr{wtbarr}

Creating from scratch

To create a WCSTransform from scratch, we can use the constructor like:

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)

where the first argument is the number of axes naxis and keyword arguments can also be passed to set various attributes of the transform. Specifying keyword arguments is equivalent to setting them after construction and the above is equivalent to:

julia> wcs = WCSTransform(2)
julia> wcs.cdelt = [-0.066667, 0.066667]
julia> wcs.ctype = ["RA---AIR", "DEC--AIR"]
julia> wcs.crpix = [1000., 1000.]
julia> wcs.crval = [0., -90]
julia> wcs.pv = [(2, 1, 45.0)]

Both will return a WCSTransform object.

Creating from a Header

Currently, we can use from_header(header) to convert the header which is in string format to a WCSTransform object. We can also pass a keyword argument relax which 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.

There exisits a scope to make a WCSTransform object from a FITSHeader object and vice-versa.

Converting a WCSTransform back to Header

We can write a WCSTransform object back to a FITS header in form of a string by using to_header(wcs). The relax keyword can be used here too for controling how non-standard extensions to the WCS standard are handled.

Transformations

1) Pixel Coordinate to World Coordinate

Pixel Coordinate: Coordinates of each pixel in your N dimensional image. World Coordinate: Coordinates that serve to locate a measurement in some multidimensional parameter space. They include, for example, a measurable quantity such as the frequency or wavelength associated with each point in a spectrum or, more abstractly, the longitude and latitude in a conventional spherical coordinate system which define a direction in space.

Consider the above wcs object, the transformation can be done like:

julia> pixcoords = [0.0  24.0  45.0;  # x coordinates of a 3D image
                    0.0  38.0  98.0]  # y coordinates of a 3D image

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

worldcoords represents the World Coordinates of the corresponding pixel coordinates specified according to the wcs data.

2) World Coordinate to Pixel Coordinate

We can simply use:

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

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.

Interacting with WCS Transform

To get the attributes of a WCSTransform object, we can use get_property(wcs, :attribute) as:

julia> getproperty(wcs, :cdelt)
2-element Array{Float64,1}:
 -0.066667
  0.066667

Similarly, to set some attribute of a WCSTransform object, we can use setproperty!(wcs, :attribute, value) as: julia> setproperty!(wcs, :cdelt, [-0.0003, 0.0002])

Other features

check what underlying C library version is being used.
```
julia> WCS.wcslib_version()
v"6.2.0"
```

Rohit Kumar