The FITS standard defines the world coordinate system (WCS) as a mechanism to associate physical values to positions within a dataset. For example, it can be used to convert pixel coordinates in an image to celestial coordinates like the right ascension and declination. The functions in this section are mainly just wrappers over CFITSIO, WCSLIB and GSL library functions to help in common applications.

[**Tread safety**] Since WCSLIB version 5.18 (released in January 2018), most WCSLIB functions are thread safe^{255}.
Gnuastro has high-level functions to easily spin-off threads and speed up your programs.
For a fully working example see Library demo - multi-threaded operation.
However you still need to be cautious in the following scenarios below.

- Many users or operating systems may still use an older version.
- The
`wcsprm`

structure of WCSLIB is not thread-safe: you can’t use the same pointer on multiple threads. For example, if you use`gal_wcs_img_to_world`

simultaneously on multiple threads, you shouldn’t pass the same`wcsprm`

structure pointer. You can use`gal_wcs_copy`

to keep and use separate copies the main structure within each thread, and later free the copies with`gal_wcs_free`

.

The full set of functions and global constants that are defined by Gnuastro’s `gnuastro/wcs.h` are described below.

- Global integer:
**GAL_WCS_DISTORTION_TPD**¶ - Global integer:
**GAL_WCS_DISTORTION_SIP**¶ - Global integer:
**GAL_WCS_DISTORTION_TPV**¶ - Global integer:
**GAL_WCS_DISTORTION_DSS**¶ - Global integer:
**GAL_WCS_DISTORTION_WAT**¶ - Global integer:
**GAL_WCS_DISTORTION_INVALID**¶ -
Gnuastro identifiers of the various WCS distortion conventions, for more, see Calabretta et al. (2004, preprint)

^{256}. Among these, SIP is a prior distortion, the rest other are sequent distortions. TPD is a superset of all these, hence it has both prior and sequeal distortion coefficients. More information is given in the documentation of`dis.h`

, from the WCSLIB manual^{257}.

- Global integer:
**GAL_WCS_COORDSYS_EQB1950**¶ - Global integer:
**GAL_WCS_COORDSYS_EQJ2000**¶ - Global integer:
**GAL_WCS_COORDSYS_ECB1950**¶ - Global integer:
**GAL_WCS_COORDSYS_ECJ2000**¶ - Global integer:
**GAL_WCS_COORDSYS_GALACTIC**¶ - Global integer:
**GAL_WCS_COORDSYS_SUPERGALACTIC**¶ - Global integer:
**GAL_WCS_COORDSYS_INVALID**¶ -
Recognized WCS coordinate systems in Gnuastro.

`EQ`

and`EC`

stand for the EQuatorial and ECliptic coordinate systems. In the equatorial and ecliptic coordinates,`B1950`

stands for the Besselian 1950 epoch and`J2000`

stands for the Julian 2000 epoch.

- Global integer:
**GAL_WCS_LINEAR_MATRIX_PC**¶ - Global integer:
**GAL_WCS_LINEAR_MATRIX_CD**¶ - Global integer:
**GAL_WCS_LINEAR_MATRIX_INVALID**¶ Identifiers of the linear transformation matrix: either in the

`PCi_j`

or the`CDi_j`

formalism. For more, see the description of`--wcslinearmatrix`in Input/Output options.

- Macro:
**GAL_WCS_FLTERROR**¶ Limit of rounding for floating point errors.

- Function:

`struct wcsprm *`

**gal_wcs_create**`(double`

¶`*crpix`

, double`*crval`

, double`*cdelt`

, double`*pc`

, char`**cunit`

, char`**ctype`

, size_t`ndim`

, int`linearmatrix`

) Given all the most common standard components of the WCS standard, construct a

`struct wcsprm`

, initialize and set it for future processing. See the FITS WCS standard for more on these keywords. All the arrays must have`ndim`

elements with them except for`pc`

which should have`ndim*ndim`

elements (a square matrix). Also,`cunit`

and`ctype`

are arrays of strings. If`GAL_WCS_LINEAR_MATRIX_CD`

is passed to`linearmatrix`

then the output WCS structure will have a CD matrix (even though you have given a PC and CDELT matrix as input to this function). Otherwise, the output will have a PC and CDELT matrix (which is the recommended format by WCSLIB).#include <stdio.h> #include <stdlib.h> #include <gnuastro/wcs.h> int main(void) { int status; size_t ndim=2; struct wcsprm *wcs; double crpix[]={50, 50}; double pc[]={-1, 0, 0, 1}; double cdelt[]={0.4, 0.4}; double crval[]={178.23, 36.98}; char *cunit[]={"deg", "deg"}; char *ctype[]={"RA---TAN", "DEC--TAN"}; int linearmatrix = GAL_WCS_LINEAR_MATRIX_PC; /* Allocate and fill the 'wcsprm' structure. */ wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit, ctype, ndim, linearmatrix); printf("WCS structure created.\n"); /*... Add any operation with the WCS structure here ...*/ /* Free the WCS structure. */ gal_wcs_free(wcs); printf("WCS structure freed.\n"); /* Return successfully. */ return EXIT_SUCCESS; }

- Function:

`struct wcsprm *`

**gal_wcs_read_fitsptr**`(fitsfile`

¶`*fptr`

, int`linearmatrix`

, size_t`hstartwcs`

, size_t`hendwcs`

, int`*nwcs`

) Return the WCSLIB

`wcsprm`

structure that is read from the CFITSIO`fptr`

pointer to an opened FITS file. With older WCSLIB versions (in particular below version 5.18) this function may not be thread-safe.Also put the number of coordinate representations found into the space that

`nwcs`

points to. To read the WCS structure directly from a filename, see`gal_wcs_read`

below. After processing has finished, you should free the WCS structure that this function returns with`gal_wcs_free`

.The

`linearmatrix`

argument takes one of three values:`0`

,`GAL_WCS_LINEAR_MATRIX_PC`

and`GAL_WCS_LINEAR_MATRIX_CD`

. It will determine the format of the WCS when it is later written to file with`gal_wcs_write`

or`gal_wcs_write_in_fitsptr`

(which is called by`gal_fits_img_write`

) So if you do not want to write the WCS into a file later, just give it a value of`0`

. For more on the difference between these modes, see the description of`--wcslinearmatrix`in Input/Output options.If you do not want to search the full FITS header for WCS-related FITS keywords (for example, due to conflicting keywords), but only a specific range of the header keywords you can use the

`hstartwcs`

and`hendwcs`

arguments to specify the keyword number range (counting from zero). If`hendwcs`

is larger than`hstartwcs`

, then only keywords in the given range will be checked. Hence, to ignore this feature (and search the full FITS header), give both these arguments the same value.If the WCS information could not be read from the FITS file, this function will return a

`NULL`

pointer and put a zero in`nwcs`

. A WCSLIB error message will also be printed in`stderr`

if there was an error.This function is just a wrapper over WCSLIB’s

`wcspih`

function which is not thread-safe. Therefore, be sure to not call this function simultaneously (over multiple threads).

- Function:

`struct wcsprm *`

**gal_wcs_read**`(char`

¶`*filename`

, char`*hdu`

, int`linearmatrix`

, size_t`hstartwcs`

, size_t`hendwcs`

, int`*nwcs`

) [

**Not thread-safe**] Return the WCSLIB structure that is read from the HDU/extension`hdu`

of the file`filename`

. Also put the number of coordinate representations found into the space that`nwcs`

points to. Please see`gal_wcs_read_fitsptr`

for more.After processing has finished, you should free the WCS structure that this function returns with

`gal_wcs_free`

.

- Function:

`void`

**gal_wcs_free**`(struct wcsprm`

¶`*wcs`

) Free the contents

*and*the space that`wcs`

points to. WCSLIB’s`wcsfree`

function only frees the contents of the`wcsprm`

structure, not the actual pointer. However, Gnuastro’s`wcsprm`

creation and reading functions allocate the structure also. This higher-level function therefore simplifies the processing. A complete working example is given in the description of`gal_wcs_create`

.

- Function:

`char *`

**gal_wcs_dimension_name**`(struct wcsprm`

¶`*wcs`

, size_t`dimension`

) Return an allocated string array (that should be freed later) containing the first part of the

`CTYPEi`

FITS keyword (which contains the dimension name in the FITS standard). For example, if`CTYPE1`

is`RA---TAN`

, the string that function returns will be`RA`

. Recall that the second component of`CTYPEi`

contains the type of projection.

- Function:

`char *`

**gal_wcs_write_wcsstr**`(struct wcsprm`

¶`*wcs`

, int`*nkeyrec`

) Return an allocated string which contains the respective FITS keywords for the given WCS structure into it. The number of keywords is written in the space pointed by

`nkeyrec`

. Each FITS keyword is 80 characters wide (according to the FITS standard), and the next one is placed immediately after it, so the full string has`80*nkeyrec`

bytes. The output of this function can later be written into an opened FITS file using`gal_fits_key_write_wcsstr`

(see FITS header keywords).

- Function:

`void`

**gal_wcs_write**`(struct wcsprm`

¶`*wcs`

, char`*filename`

, char`*extname`

, gal_fits_list_key_t`*headers`

, char`*program_string`

) Write the given WCS structure into the second extension of an empty FITS header. The first/primary extension will be empty like the default format of all Gnuastro outputs. When

`extname!=NULL`

it will be used as the FITS extension name. Any set of extra headers can also be written through the`headers`

list and if`program_string!=NULL`

it will be used in a commented keyword title just above the written version information.

- Function:

`void`

**gal_wcs_write_in_fitsptr**`(fitsfile`

¶`*fptr`

, struct wcsprm`*wcs`

) Convert the input

`wcs`

structure (keeping the WCS programmatically) into FITS keywords and write them into the given FITS file pointer. This is a relatively low-level function which assumes the FITS file has already been opened with CFITSIO. If you just want to write the WCS into an empty file, you can use`gal_wcs_write`

(which internally calls this function after creating the FITS file and later closes it safely).

- Function:

`struct wcsprm *`

**gal_wcs_copy**`(struct wcsprm`

¶`*wcs`

) Return a fully allocated (independent) copy of

`wcs`

.

- Function:

`void`

**gal_wcs_remove_dimension**`(struct wcsprm`

¶`*wcs`

, size_t`fitsdim`

) Remove the given FITS dimension from the given

`wcs`

structure.

- Function:

`void`

**gal_wcs_on_tile**`(gal_data_t`

¶`*tile`

) Create a WCSLIB

`wcsprm`

structure for`tile`

using WCS parameters of the tile’s allocated block dataset, see Tessellation library (`tile.h`) for the definition of tiles. If`tile`

already has a WCS structure, this function will not do anything.In many cases, tiles are created for internal/low-level processing. Hence for performance reasons, when creating the tiles they do not have any WCS structure. When needed, this function can be used to add a WCS structure to each tile tile by copying the WCS structure of its block and correcting the reference point’s coordinates within the tile.

- Function:

`double *`

**gal_wcs_warp_matrix**`(struct wcsprm`

¶`*wcs`

) Return the Warping matrix of the given WCS structure as an array of double precision floating points. This will be the final matrix, irrespective of the type of storage in the WCS structure. Recall that the FITS standard has several methods to store the matrix. The output is an allocated square matrix with each side equal to the number of dimensions.

- Function:

`void`

**gal_wcs_clean_small_errors**`(struct wcsprm`

¶`*wcs`

) Errors can make small differences between the pixel-scale elements (

`CDELT`

) and can also lead to extremely small values in the`PC`

matrix. With this function, such errors will be “cleaned” as follows: 1) if the maximum difference between the`CDELT`

elements is smaller than the reference error, it will be set to the mean value. When the FITS keyword`CRDER`

(optional) is defined it will be used as a reference, if not the default value is`GAL_WCS_FLTERROR`

. 2) If any of the PC elements differ from 0, 1 or -1 by less than`GAL_WCS_FLTERROR`

, they will be rounded to the respective value.

- Function:

`void`

**gal_wcs_decompose_pc_cdelt**`(struct wcsprm`

¶`*wcs`

) Decompose the

`PCi_j`

and`CDELTi`

elements of`wcs`

. According to the FITS standard, in the`PCi_j`

WCS formalism, the rotation matrix elements \(m_{ij}\) are encoded in the`PCi_j`

keywords and the scale factors are encoded in the`CDELTi`

keywords. There is also another formalism (the`CDi_j`

formalism) which merges the two into one matrix.However, WCSLIB’s internal operations are apparently done in the

`PCi_j`

formalism. So its outputs are also all in that format by default. When the input is a`CDi_j`

, WCSLIB will still read the matrix directly into the`PCi_j`

matrix and the`CDELTi`

values are set to`1`

(one). This function is designed to correct such issues: after it is finished, the`CDELTi`

values in`wcs`

will correspond to the pixel scale, and the`PCi_j`

will correction show the rotation.

- Function:

`void`

**gal_wcs_to_cd**`(struct wcsprm`

¶`*wcs`

) Make sure that the WCS structure’s

`PCi_j`

and`CDi_j`

keywords have the same value and that the`CDELTi`

keywords have a value of 1.0. Also, set the`wcs->altlin=2`

(for the`CDi_j`

formalism). With these changes`gal_wcs_write_in_fitsptr`

(and thus`gal_wcs_write`

and`gal_fits_img_write`

and its derivatives) will have an output file in the format of`CDi_j`

.

- Function:

`int`

**gal_wcs_coordsys_from_string**`(char`

¶`*coordsys`

) Convert the given string to Gnuastro’s integer-based WCS coordinate system identifier (one of the

`GAL_WCS_COORDSYS_*`

, listed above). The expected strings can be seen in the description of the`--wcscoordsys`option of the Fits program, see Keyword inspection and manipulation.

- Function:

`int`

**gal_wcs_coordsys_identify**`(struct wcsprm`

¶`*wcs`

) Read the given WCS structure and return its coordinate system as one of Gnuastro’s WCS coordinate system identifiers (the macros

`GAL_WCS_COORDSYS_*`

, listed above).

- Function:

`struct wcsprm *`

**gal_wcs_coordsys_convert**`(struct wcsprm`

¶`*inwcs`

, int`coordsysid`

) Return a newly allocated WCS structure with the

`coordsysid`

coordinate system identifier. The Gnuastro WCS distortion identifiers are defined in the`GAL_WCS_COORDSYS_*`

macros mentioned above. Since the returned dataset is newly allocated, if you do not need the original dataset after this, use the WCSLIB library function`wcsfree`

to free the input, for example,`wcsfree(inwcs)`

.

- Function:

`int`

**gal_wcs_distortion_from_string**`(char`

¶`*distortion`

) Convert the given string (assumed to be a FITS-standard, string-based distortion identifier) to a Gnuastro’s integer-based distortion identifier (one of the

`GAL_WCS_DISTORTION_*`

macros defined above). The sting-based distortion identifiers have three characters and are all in capital letters.

- Function:

`int`

**gal_wcs_distortion_to_string**`(int`

¶`distortion`

) Convert the given Gnuastro integer-based distortion identifier (one of the

`GAL_WCS_DISTORTION_*`

macros defined above) to the string-based distortion identifier) of the FITS standard. The sting-based distortion identifiers have three characters and are all in capital letters.

- Function:

`int`

**gal_wcs_distortion_identify**`(struct wcsprm`

¶`*wcs`

) Returns the Gnuastro identifier for the distortion of the input WCS structure. The returned value is one of the

`GAL_WCS_DISTORTION_*`

macros defined above. When the input pointer to a structure is`NULL`

, or it does not contain a distortion, the returned value will be`GAL_WCS_DISTORTION_INVALID`

.

- Function:

`struct wcsprm *`

**gal_wcs_distortion_convert(struct**`wcsprm`

¶`*inwcs`

, int`outdisptype`

, size_t`*fitsize`

) Return a newly allocated WCS structure, where the distortion is implemented in a different standard, identified by the identifier

`outdisptype`

. The Gnuastro WCS distortion identifiers are defined in the`GAL_WCS_DISTORTION_*`

macros mentioned above.The available conversions in this function will grow. Currently it only supports converting TPV to SIP and vice versa, following the recipe of Shupe et al. (2012)

^{258}. Please get in touch with us if you need other types of conversions.For some conversions, direct analytical conversions do not exist. It is thus necessary to model and fit the two types. In such cases, it is also necessary to specify the

`fitsize`

array that is the size of the array along each C-ordered dimension, so you can simply pass the`dsize`

element of your`gal_data_t`

dataset, see Generic data container (`gal_data_t`

). Currently this is only necessary when converting TPV to SIP. For other conversions you may simply pass a`NULL`

pointer.For example, if you want to convert the TPV coefficients of your input

`image.fits`to SIP coefficients, you can use the following functions (which are also available as a command-line operation in Fits).int nwcs; gal_data_t *data=gal_fits_img_read("image.fits", "1", -1, 1); inwcs=gal_wcs_read("image.fits", "1", 0, 0, 0, &nwcs); data->wcs=gal_wcs_distortion_convert(inwcs, GAL_WCS_DISTORTION_TPV, NULL); wcsfree(inwcs); gal_fits_img_write(data, "tpv.fits", NULL, NULL);

- Function:

`double`

**gal_wcs_angular_distance_deg**`(double`

¶`r1`

, double`d1`

, double`r2`

, double`d2`

) Return the angular distance (in degrees) between a point located at (

`r1`

,`d1`

) to (`r2`

,`d2`

). All input coordinates are in degrees. The distance (along a great circle) on a sphere between two points is calculated with the equation below.$$\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)$$

However, since the pixel scales are usually very small numbers, this function will not use that direct formula. It will be use the Haversine formula which is better considering floating point errors:

$${\sin^2(d)\over 2}=\sin^2\left( {d_1-d_2\over 2} \right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)$$

- Function:

`void`

**gal_wcs_box_vertices_from_center**`(double`

¶`ra_center`

, double`dec_center`

, double`ra_delta`

, double`dec_delta`

, double`*out`

) Calculate the vertices of a rectangular box given the central RA and Dec and delta of each. The vertice coordinates are written in the space that

`out`

points to (assuming it has space for eight`double`

s).Given the spherical nature of the coordinate system, the vertice lengths can’t be calculated with a simple addition/subtraction. For the declination, a simple addition/subtraction is enough. Also, on the equator (where the RA is defined), a simple addition/subtraction along the RA is fine. However, at other declinations, the new RA after a shift needs special treatment, such that close to the poles, a shift of 1 degree can correspond to a new RA that is much more distant than the original RA. Assuming a point at Right Ascension (RA) and Declination of \(\alpha\) and \(\delta\), a shift of \(R\) degrees along the positive RA direction corresponds to a right ascension of \(\alpha+\frac{R}{\cos(\delta)}\). For more, see the description of

`box-vertices-on-sphere`

in Box shape operators.The 8 coordinates of the 4 vertices of the box are written in the order below. Where “bottom” corresponds to a lower declination and “top” to higher declination, “left” corresponds to a larger RA and “right” corresponds to lower RA.

out[0]: bottom-left RA out[1]: bottom-left Dec out[2]: bottom-right RA out[3]: bottom-right Dec out[4]: top-right RA out[5]: top-right Dec out[6]: top-left RA out[7]: top-left Dec

- Function:

`double *`

**gal_wcs_pixel_scale**`(struct wcsprm`

¶`*wcs`

) Return the pixel scale for each dimension of

`wcs`

in degrees. The output is an allocated array of double precision floating point type with one element for each dimension. If it is not successful, this function will return`NULL`

.

- Function:

`double`

**gal_wcs_pixel_area_arcsec2**`(struct wcsprm`

¶`*wcs`

) Return the pixel area of

`wcs`

in arc-second squared. This only works when the input dataset has at least two dimensions and the units of the first two dimensions (`CUNIT`

keywords) are`deg`

(for degrees). In other cases, this function will return a NaN.

- Function:

`int`

**gal_wcs_coverage**`(char`

¶`*filename`

, char`*hdu`

, size_t`*ondim`

, double`**ocenter`

, double`**owidth`

, double`**omin`

, double`**omax`

) Find the sky coverage of the image HDU (

`hdu`

) within`filename`. The number of dimensions is written into`ndim`

, and space for the various output arrays is internally allocated and filled with the respective values. Therefore you need to free them afterwards.Currently this function only supports images that are less than 180 degrees in width (which is usually the case!). This requirement has been necessary to account for images that cross the RA=0 hour circle on the sky. Please get in touch with us at mailto:bug-gnuastro@gnu.org if you have an image that is larger than 180 degrees so we try to find a solution based on need.

- Function:

`gal_data_t *`

**gal_wcs_world_to_img**`(gal_data_t`

¶`*coords`

, struct wcsprm`*wcs`

, int`inplace`

) Convert the linked list of world coordinates in

`coords`

to a linked list of image coordinates given the input WCS structure.`coords`

must be a linked list of data structures of float64 (‘double’) type, seeLinked lists (`list.h`) and List of`gal_data_t`

. The top (first popped/read) node of the linked list must be the first WCS coordinate (RA in an image usually) etc. Similarly, the top node of the output will be the first image coordinate (in the FITS standard). In case WCSLIB fails to convert any of the coordinates (for example, the RA of one coordinate is given as 400!), the respective element in the output will be written as NaN.If

`inplace`

is zero, then the output will be a newly allocated list and the input list will be untouched. However, if`inplace`

is non-zero, the output values will be written into the input’s already allocated array and the returned pointer will be the same pointer to`coords`

(in other words, you can ignore the returned value). Note that in the latter case, only the values will be changed, things like units or name (if present) will be untouched.

- Function:

`gal_data_t *`

**gal_wcs_img_to_world**`(gal_data_t`

¶`*coords`

, struct wcsprm`*wcs`

, int`inplace`

) Convert the linked list of image coordinates in

`coords`

to a linked list of world coordinates given the input WCS structure. See the description of`gal_wcs_world_to_img`

for more details.

https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/threads.html

https://www.atnf.csiro.au/people/mcalabre/WCS/dcs_20040422.pdf

https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/dis_8h.html

Proc. of SPIE Vol. 8451 84511M-1. https://doi.org/10.1117/12.925460, also available at http://web.ipac.caltech.edu/staff/shupe/reprints/SIP_to_PV_SPIE2012.pdf.

JavaScript license information

GNU Astronomy Utilities 0.20 manual, April 2023.