GNU Astronomy Utilities



12.4.6 Library demo - Warp to new grid

Gnuastro’s warp library (that you can access by including gnuastro/warp.h) allows you to resample an image from a grid to another entirely using the WCSLIB (while accounting for distortions if necessary; see Warp library (warp.h)). The Warp library uses a pixel-mixing or area-based resampling approach which is fully described in Resampling. The most generic uses cases for this library are already available in the Invoking Warp program. For a related demo (where the output grid and WCS are imported from another file), see Library demo - Warp to another image.

In the example below, we’ll assume you have the SDSS image downloaded in Downloading and validating input data. After downloading the image as described there, you will have r.fits in your current directory. We will therefore use r.fits as the input to the rest program here. The image is not aligned to the celestial coordinates, so we will align the pixel and WCS coordinates, but set the center of the pixel grid to be at (RA,Dec) of (202.4173735,47.3374525). We also give it a TAN projection with a pixel scale of 0.27 arcsecs, a defined center pixel. However, we’ll let the Warp library measure the proper output image size that will contain the aligned image.

To compile the demonstration program below, copy and paste the contents in a plain-text file (let’s assume you named it align-to-new.c) and use BuildProgram with this command: ‘astbuildprog align-to-new.c’. Please note that the demo program does not perform many sanity checks to avoid making it too complex and to highlight this particular feature in the library. For a robust method write programs with all the necessary sanity checks, see Gnuastro’s Warp source code, see Program source.

To encourage good coding practices, this script contains a copyright notice with a place holder for your name and your email (as you customize it for your own purpose). Always keep a one-line description and copyright notice like this in all your scripts, such “metadata” is very important to accompany every source file you write. Of course, when you write the source file from scratch and just learn how to use a single function from this manual, only your name/year should appear. The existing name of the original author of this example program is only for cases where you copy-paste this whole file.

/* Warp an image to a new grid.
 *
 * In the example below, We will use 'r.fits' as the input. The image is
 * not aligned to the celestial coordinates, so we will align the pixel
 * and WCS coordinates. We also give it a TAN projection. However, we’ll
 * let the Warp library measure the proper output image size that will
 * contain the aligned image.
 *
 * Copyright (C) 2024      Your Name <your@@email.address>
 * Copyright (C) 2022-2024 Pedram Ashofteh-Ardakani <pedramardakani@pm.me>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>

#include <gnuastro/wcs.h>       /* Contains gnuastro's fits.h */
#include <gnuastro/warp.h>      /* Contains gnuastro's data.h */
#include <gnuastro/array.h>     /* Contains gnuastro's type.h */

int
main(void)
{
  /* Input file's name and HDU. */
  char *filename="r.fits", *hdu="0";

  /* Output file name. */
  char *outname="align-to-new.fits";

  /* RA/Dec of the center of the central pixel of output. Please
   * change the center based on your input. */
  double center[]={202.4173735, 47.3374525};

  /* Coordinate and Projection algorithms of output. */
  char *ctype[2]={"RA---TAN", "DEC--TAN"};

  /* Output pixel scale (in units of degrees/pixel). */
  double cdelt[]={0.27/3600, 0.27/3600};

  /* For intermediate steps. */
  size_t two=2;

  /* Initialize the 'wa' struct with empty values and NULL pointers. */
  gal_warp_wcsalign_t wa=gal_warp_wcsalign_template();

  /* Set the width (and height!) of the output in pixels (as a 1D and
   * 2 element 'gal_data_t'). When it is NULL, the library will
   * calculate the appropriate width to fully fit the input image
   * after alignment. */
  wa.widthinpix=NULL;

  /* Set the number of threads to use. If the value is '0', the
   * library will estimate the maximum available threads at
   * run-time on the host operating system. */
  wa.numthreads=0;


  /* Read the input image and its WCS. */
  wa.input=gal_array_read_one_ch_to_type(filename, hdu, NULL,
					 GAL_TYPE_FLOAT64, -1, 0, NULL);
  wa.input->wcs=gal_wcs_read(filename, hdu, 0, 0, 0, &wa.input->nwcs,
                             NULL);


  /* Prepare the warp input structure. */
  wa.coveredfrac=1; wa.edgesampling=0;
  wa.ctype=gal_data_alloc(ctype, GAL_TYPE_STRING, 1, &two, NULL, 1,
                          -1, 0, NULL, NULL, NULL);
  wa.cdelt=gal_data_alloc(cdelt, GAL_TYPE_FLOAT64, 1, &two, NULL, 1,
                          -1, 0, NULL, NULL, NULL);
  wa.center=gal_data_alloc(center, GAL_TYPE_FLOAT64, 1, &two, NULL, 1,
                           -1, 0, NULL, NULL, NULL);


  /* Do the warp, then convert it to a 32-bit float. */
  gal_warp_wcsalign(&wa);
  wa.output=gal_data_copy_to_new_type_free(wa.output, GAL_TYPE_FLOAT32);


  /* WARNING: make sure there is no file with same name as 'out.fits'
   * or the result will be appended to its final HDU. */
  gal_fits_img_write(wa.output, outname, NULL, 0);


  /* Remove the pointers to arrays that we didn't allocate (and thus,
   * should not be freed by 'gal_data_free' below). */
  wa.cdelt->array=wa.center->array=wa.ctype->array=NULL;


  /* Clean up. */
  gal_data_free(wa.cdelt);   gal_data_free(wa.ctype);
  gal_data_free(wa.input);   gal_data_free(wa.output);
  gal_data_free(wa.center);  gal_data_free(wa.widthinpix);

  /* Give control back to the operating system. */
  return EXIT_SUCCESS;
}