GNU Astronomy Utilities


Next: , Previous: , Up: Gnuastro library   [Contents][Index]


10.3.13 World Coordinate System (wcs.h)

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.

Macro: GAL_WCS_DISTORTION_TPD
Macro: GAL_WCS_DISTORTION_SIP
Macro: GAL_WCS_DISTORTION_TPV
Macro: GAL_WCS_DISTORTION_DSS
Macro: GAL_WCS_DISTORTION_WAT
Macro: GAL_WCS_DISTORTION_INVALID

Gnuastro identifiers of the various WCS distortion conventions, for more, see Calabretta et al. (2004, preprint)183. 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 manual184.

Macro: GAL_WCS_FLTERROR

Limit of rounding for floating point errors.

Function:
struct wcsprm *
gal_wcs_read_fitsptr (fitsfile *fptr, size_t hstartwcs, size_t hendwcs, int *nwcs)

[Not thread-safe] Return the WCSLIB wcsprm structure that is read from the CFITSIO fptr pointer to an opened FITS file. 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 can free the returned structure with WCSLIB’s wcsvfree keyword:

status = wcsvfree(&nwcs,&wcs);

If you don’t 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 couldn’t 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, 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.

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 won’t do anything.

In many cases, tiles are created for internal/low-level processing. Hence for performance reasons, when creating the tiles they don’t 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:
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 doesn’t 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)185. Please get in touch with us if you need other types of conversions.

For some conversions, direct analytical conversions don’t 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, &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 won’t 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:
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 its not successful, this function will return NULL.

Function:
double
gal_wcs_pixel_area_arcsec2 (struct wcsprm *wcs)

Return the pixel area of wcs in arcsecond squared. If the input WCS structure is not two dimensional and the units (CUNIT keywords) are not deg (for degrees), then this function will return a NaN.

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).

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.


Footnotes

(183)

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

(184)

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

(185)

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.


Next: , Previous: , Up: Gnuastro library   [Contents][Index]