GNU Astronomy Utilities



6.4.4.1 Align pixels with WCS considering distortions

When none of the linear warps191 are requested, Warp will align the input’s pixel axes with it’s WCS axes. In the process, any possibly existing distortion is also removed (such as TPV and SIP). Usually, the WCS axes are the Right Ascension and Declination in equatorial coordinates. The output image’s pixel grid is highly customizable through the options in this section. To learn about Warp’s strategy to build the new pixel grid, see Resampling. For strong distortions (that produce strong curvatures), you can fine-tune the area-based resampling with --edgesampling, as described below.

On the other hand, sometimes you need to Warp an input image to the exact same grid of an already available reference FITS image with an existing WCS. If that image is already aligned, finding its center, number of pixels and pixel scale can be annoying (and just increase the complexity of your script). On the other hand, if that image is not aligned (for example, has a certain rotation in the sky, and has a different distortion), there are too many WCS parameters to set (some are not yet available explicitly in the options here)! For such scenarios, Warp has the --gridfile option. When --gridfile is called, the options below that are used to define the output’s WCS will be ignored (these options: --center, --widthinpix, --cdelt, --ctype). In this case, the output’s WCS and pixel grid will exactly match the image given to --gridfile (including any rotation, pixel scale, or distortion or projection).

Set --cdelt explicitly when you plan to stack many warped images: To align some images and later stack them, it is necessary to be sure the pixel sizes of all the images are the same exactly. Most of the time the measured (during astrometry) pixel scale of the separate exposures, will be different in the second or third digit number after the decimal point. It is a normal/statistical error in measuring the astrometry. On a large image, these slight differences can cause different output sizes (of one or two pixels on a very large image).

You can fix this by explicitly setting the pixel scale of each warped exposure with Warp’s --cdelt option that is described below. For good strategies of setting the pixel scale, see Moiré pattern in stacking and its correction.

Another problem that may arise when aligning images to new pixel grids is the aliasing or visible Moiré patterns on the output image. This artifact should be removed if you are stacking several exposures, especially with a pointing pattern. If not see Moiré pattern in stacking and its correction for ways to mitigate the visible patterns. See the description of --gridfile below for more.

Known issue: Warp’s WCS-based aligning works best with WCSLIB version 7.12 (released in September 2022) and above. If you have an older version of WCSLIB, you might get a wcss2p error otherwise.

-c FLT,FLT
--center=FLT,FLT

WCS coordinates of the center of the central pixel of the output image. Since a central pixel is only defined with an odd number of pixels along both dimensions, the output will always have an odd number of pixels. When --center or --gridfile aren’t given, the output will have the same central WCS coordinate as the input.

Usually, the WCS coordinates are Right Ascension and Declination (when the first three characters of CTYPE1 and CTYPE2 are respectively RA- and DEC). For more on the CTYPEi keyword values, see --ctype below.

-w INT[,INT]
--width=INT[,INT]

Width and height of the output image in units of WCS (usually degrees). If you want the values to be read as pixels, also call the --widthinpix option with --width. If a single value is given, Warp will use the same value for the second dimension (creating a square output). When --width or --gridfile aren’t given, Warp will calculate the necessary size of the output pixel grid to fully contain the input image.

Usually the WCS coordinates are in units of degrees (defined by the CUNITi keywords of the FITS standard). But entering a certain number of arcseconds or arcminutes for the width can be annoying (you will usually need to go to the calculator!). To simplify such situations, this option also accepts division. For example --width=1/60,2/60 will make an aligned warp that is 1 arcmin along Right Ascension and 2 arcminutes along the Declination.

With the --widthinpix option the values will be interpreted as numbers of pixels. In this scenario, this option should be given odd integer(s) that are greater than 1. This ensures that the output image can have a central pixel. Recall that through the --center option, you specify the WCS coordinate of the center of the central pixel. The central coordinate of an image with an even number of pixels will be on the edge of two pixels, so a “central” pixel is not well defined. If any of the given values are even, Warp will automatically add a single pixel (to make it an odd integer) and print a warning message.

--widthinpix

When called, the values given to the --width option will be interpreted as the number of pixels along each dimension(s). See the description of --width for more.

-x FLT[,FLT]
--cdelt=FLT[,FLT]

Coordinate deltas or increments (CDELTi in the FITS standard), or the pixel scale in both dimensions. If a single value is given, it will be used for both axes. In this way, the output’s pixels will be squares on the sky at the reference point (as is usually expected!). When --cdelt or --gridfile aren’t given, Warp will read the input’s pixel scale and choose the larger of CDELT1 or CDELT2 so the output pixels are square.

Usually (when dealing with RA and Dec, and the CUNITis have a value of deg), the units of the given values are degrees/pixel. Warp allows you to easily convert from arcsec to degrees by simply appending a /3600 to the value. For example, for an output image of pixel scale 0.27 arcsec/pixel, you can use --cdelt=0.27/3600.

--ctype=STR,STR

The coordinate types of the output (CTYPE1 and CTYPE2 keywords in the FITS standard), separated by a comma. By default the value to this option is ‘RA---TAN,DEC--TAN’. However, if --gridfile is given, this option is ignored.

If you don’t call --ctype or --gridfile, the output WCS coordinates will be Right Ascension and Declination, while the output’s projection will be Gnomonic, also known as Tangential (TAN). This combination is the most common in extra-galactic imaging surveys. For other coordinates and projections in your output use other values, as described below.

According to the FITS standard version 4.0192: CTYPEi is the “type for the Intermediate-coordinate Axis \(i\). Any coordinate type that is not covered by this Standard or an officially recognized FITS convention shall be taken to be linear. All non-linear coordinate system names must be expressed in ‘4–3’ form: the first four characters specify the coordinate type, the fifth character is a hyphen (-), and the remaining three characters specify an algorithm code for computing the world coordinate value. Coordinate types with names of fewer than four characters are padded on the right with hyphens, and algorithm codes with fewer than three characters are padded on the right with SPACE. Algorithm codes should be three characters” (see list of algorithm codes below).

You can use any of the projection algorithms (last three characters of each coordinate’s type) provided by your host WCSLIB (a mandatory dependency of Gnuastro; see WCSLIB). For a very elaborate and complete description of projection algorithms in the FITS WCS standard, see Calabretta and Greisen 2002. Wikipedia also has a nice article on Map projections. As an example, WCSLIB 7.12 (released in September 2022) has the following projection algorithms:

AZP

Zenithal/azimuthal perspective

SZP

Slant zenithal perspective

TAN

Gnomonic (tangential)

STG

Stereographic

SIN

Orthographic/synthesis

ARC

Zenithal/azimuthal equidistant

ZPN

Zenithal/azimuthal polynomial

ZEA

Zenithal/azimuthal equal area

AIR

Airy

CYP

Cylindrical perspective

CEA

Cylindrical equal area

CAR

Plate carree

MER

Mercator

SFL

Sanson-Flamsteed

PAR

Parabolic

MOL

Mollweide

AIT

Hammer-Aitoff

COP

Conic perspective

COE

Conic equal area

COD

Conic equidistant

COO

Conic orthomorphic

BON

Bonne

PCO

Polyconic

TSC

Tangential spherical cube

CSC

COBE spherical cube

QSC

Quadrilateralized spherical cube

HPX

HEALPix

XPH

HEALPix polar, aka "butterfly"

-G
--gridfile

FITS filename containing the final pixel grid and WCS for the output image. The HDU/extension containing should be specified with --gridhdu or its short option -H. The HDU should contain a WCS, otherwise, Warp will abort with a crash. When this option is used, Warp will read the respective WCS and the size of the image to resample the input. Since this WCS of this HDU contains everything needed to construct the WCS the options above will be ignored when --gridfile is called: --cdelt, --center, and --widthinpix.

In the example below, let’s use this option to put the image of M51 in one survey (J-PLUS) into the pixel grid of another survey (SDSS) containing M51. The J-PLUS field of view is very large (almost \(1.5\times1.5\) deg\(^2\), in \(9500\times9500\) pixels), while the field of view of SDSS in each filter is small (almost \(0.3\times0.25\) deg\(^2\) in \(2048\times1489\) pixels). With the first two commands, we’ll first download the two images, then we’ll extract the portion of the J-PLUS image that overlaps with the SDSS image and align it exactly to SDSS’s pixel grid. Note that these are the two images that were used in two of Gnuastro’s tutorials: Building the extended PSF and Detecting large extended targets.

## Download the J-PLUS DR2 image of M51 in the r filter.
$ jplusbase="http://archive.cefca.es/catalogues/vo/siap"
$ wget $jplusbase/jplus-dr2/get_fits?id=67510 \
       -O jplus.fits.fz

## Download the SDSS image in r filter and decompress it
## (Bzip2 is not a standard FITS compression algorithm).
$ sdssbase=https://dr12.sdss.org/sas/dr12/boss/photoObj/frames
$ wget $sdssbase/301/3716/6/frame-r-003716-6-0117.fits.bz2 \
       -O sdss.fits.bz2
$ bunzip2 sdss.fits.bz2

## Warp and crop the J-PLUS image so the output exactly
## matches the SDSS pixel gid.
$ astwarp jplus.fits.fz --gridfile=sdss.fits --gridhdu=0 \
          --output=jplus-on-sdss.fits

## View the two images side-by-side:
$ astscript-fits-view sdss.fits jplus-on-sdss.fits

As the example above shows, this option can therefore be very useful when comparing images from multiple surveys. But there are other very interesting use cases also. For example, when you are making a mock dataset and need to add distortion to the image so it matches the distortion of your camera. Through --gridhdu, you can easily insert that distortion over the mock image and put the mock image in the pixel grid of an exposure.

-H
--gridhdu

The HDU/extension of the reference WCS file specified with option --wcsfile or its short version -H (see the description of --wcsfile for more).

--edgesampling=INT

Number of extra samplings along the edge of a pixel. By default the value is 0 (the output pixel’s polygon over the input will be a quadrilateral (a polygon with four edges/vertices).

Warp uses pixel mixing to derive the output pixel values. For a complete introduction, see Resampling, and in particular its later part on distortions. To account for this possible curvature due to distortion, you can use this option. For example, --edgesampling=1 will add one extra vertice in the middle of each edge of the output pixel, producing an 8-vertice polygon. Similarly, --edgesampling=5 will put 5 extra vertices along each edge, thus sampling the shape (and possible curvature) of the output pixel over an input pixel with \(4+5\times4=24\) vertice polygon. Since the polygon clipping will happen for every output pixel, a higher value to this option can significantly reduce the running speed and increase the RAM usage of Warp; so use it with caution: in most cases the default --edgesampling=0 is sufficient.

To visually inspect the curvature effect on pixel area of the input image, see option --pixelareaonwcs in Pixel information images.

--checkmaxfrac

Check each output pixel’s maximum coverage on the input data and append as the ‘MAX-FRAC’ HDU/extension to the output aligned image. This option provides an easy visual inspection for possible recurring patterns or fringes caused by aligning to a new pixel grid. For more detail about the origin of these patterns and how to mitigate them see Moiré pattern in stacking and its correction.

Note that the ‘MAX-FRAC’ HDU/extension is not showing the patterns themselves; It represents the largest area coverage on the input data for that particular pixel. The values can be in the range between 0 to 1, where 1 means the pixel is covering at least one complete pixel of the input data. On the other hand, 0 means that the pixel is not covering any pixels of the input at all.


Footnotes

(191)

For linear warps, see Linear warps to be called explicitly.

(192)

FITS standard version 4.0: https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf