File: wcsgrid.c

package info (click to toggle)
wcslib 7.4%2Bds-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 9,752 kB
  • sloc: ansic: 32,656; lex: 9,281; fortran: 6,634; sh: 3,369; sed: 497; pascal: 188; makefile: 15
file content (342 lines) | stat: -rw-r--r-- 9,496 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
/*============================================================================
  WCSLIB 7.4 - an implementation of the FITS WCS standard.
  Copyright (C) 1995-2021, Mark Calabretta

  This file is part of WCSLIB.

  WCSLIB is free software: you can redistribute it and/or modify it under the
  terms of the GNU Lesser General Public License as published by the Free
  Software Foundation, either version 3 of the License, or (at your option)
  any later version.

  WCSLIB 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 Lesser General Public License for
  more details.

  You should have received a copy of the GNU Lesser General Public License
  along with WCSLIB.  If not, see http://www.gnu.org/licenses.

  Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
  http://www.atnf.csiro.au/people/Mark.Calabretta
  $Id: wcsgrid.c,v 7.4 2021/01/31 02:24:52 mcalabre Exp $
*=============================================================================
*
* wcsgrid extracts the WCS keywords for an image from the specified FITS file
* and uses pgsbox() to plot a 2-D coordinate graticule for each representation
* found.  Refer to the usage notes below.
*
* TODO
*   Subimaging option.
*---------------------------------------------------------------------------*/

char usage[] =
"Usage: wcsgrid [-a<alt>] [-d<pgdev>] [-h<hdu>] [<fitsfile>]\n"
"\n"
"wcsgrid extracts the WCS keywords for an image from the specified FITS\n"
"file and uses pgsbox() to plot a 2-D coordinate graticule for each\n"
"alternate representation found.\n\n"
"The FITS file may be specified according to the syntax understood by\n"
"cfitsio, for example \"file.fits.gz+1\" refers to the first extension of\n"
"a gzip'd FITS file.  Use \"-\" or omit the file name for input from stdin.\n"
"\n"
"Options:\n"
"  -a<alt>      Plot a graticule only for the alternate representation\n"
"               specified (ignored if there is only one).\n"
"  -d<pgdev>    PGPLOT device type (default XWINDOW, use \"?\" for list).\n"
"  -h<hdu>      Move to HDU number (1-relative) which is expected to\n"
"               contain an image array.  (Useful for input from stdin.)\n";

#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include <fitsio.h>
#include <cpgplot.h>
#include <cpgsbox.h>

#include <wcshdr.h>
#include <wcsfix.h>
#include <wcsutil.h>
#include <wcs.h>
#include <getwcstab.h>


int main(int argc, char **argv)

{
  char alt = '\0', *header, idents[3][80], *infile, keyword[32], nlcprm[1],
       opt[2], pgdev[16];
  int  c0[] = {-1, -1, -1, -1, -1, -1, -1};
  int  alts[27], gcode[2], hdunum = 1, hdutype, i, ic, naxes, naxis[2],
       nkeyrec, nreject, nwcs, stat[NWCSFIX], status;
  float  blc[2], trc[2];
  double cache[257][4], grid1[3], grid2[3], nldprm[1];
  struct wcsprm *wcs;
  nlfunc_t pgwcsl_;
  fitsfile *fptr;


  // Parse options.
  strcpy(pgdev, "/XWINDOW");
  for (i = 1; i < argc && argv[i][0] == '-'; i++) {
    if (!argv[i][1]) break;

    switch (argv[i][1]) {
    case 'a':
      alt = toupper(argv[i][2]);
      break;

    case 'd':
      if (argv[i][2] == '?') {
        cpgldev();
        return 0;
      }

      if (argv[i][2] == '/') {
        strncpy(pgdev+1, argv[i]+3, 12);
      } else {
        strncpy(pgdev+1, argv[i]+2, 12);
      }
      wcsutil_null_fill(-16, pgdev);

      break;

    case 'h':
      hdunum = atoi(argv[i]+2);
      break;

    default:
      fprintf(stderr, "%s", usage);
      return 1;
    }
  }

  if (i < argc) {
    infile = argv[i++];

    if (i < argc) {
      fprintf(stderr, "%s", usage);
      return 1;
    }
  } else {
    infile = "-";
  }

  // Check accessibility of the input file.
  if (strcmp(infile, "-") && access(infile, R_OK) == -1) {
    printf("wcsgrid: Cannot access %s.\n", infile);
    return 1;
  }


  // Open the FITS file and move to the required HDU.
  status = 0;
  if (fits_open_file(&fptr, infile, READONLY, &status)) goto fitserr;
  if (fits_movabs_hdu(fptr, hdunum, &hdutype, &status)) goto fitserr;
  if (hdutype != IMAGE_HDU) {
    fprintf(stderr, "ERROR, HDU number %d does not contain an image array.\n",
      hdunum);
    return 1;
  }

  // Check that we have at least two image axes.
  if (fits_read_key(fptr, TINT, "NAXIS",  &naxes, NULL, &status)) {
    goto fitserr;
  }

  if (naxes < 2) {
    fprintf(stderr, "ERROR, HDU number %d does not contain a 2-D image.\n",
      hdunum);
    return 1;
  } else if (naxes > 2) {
    printf("HDU number %d contains a %d-D image array.\n", hdunum, naxes);
  }

  // Read in the FITS header, excluding COMMENT and HISTORY keyrecords.
  if (fits_hdr2str(fptr, 1, NULL, 0, &header, &nkeyrec, &status)) {
    goto fitserr;
  }


  // Interpret the WCS keywords.
  if ((status = wcspih(header, nkeyrec, WCSHDR_all, -3, &nreject, &nwcs,
                       &wcs))) {
    fprintf(stderr, "wcspih ERROR %d: %s.\n", status, wcshdr_errmsg[status]);
    return 1;
  }
  free(header);

  // Read -TAB arrays from the binary table extension (if necessary).
  if (fits_read_wcstab(fptr, wcs->nwtb, (wtbarr *)wcs->wtb, &status)) {
    goto fitserr;
  }

  // Translate non-standard WCS keyvalues.
  if ((status = wcsfix(7, 0, wcs, stat))) {
    status = 0;
    for (i = 0; i < NWCSFIX; i++) {
      if (stat[i] > 0) {
         fprintf(stderr, "wcsfix ERROR %d: %s.\n", stat[i],
                 wcsfix_errmsg[stat[i]]);

        // Ignore problems with CDi_ja and DATE-OBS.
        if (!(i == CDFIX || i == DATFIX)) status = 1;
      }
    }

    if (status) return 1;
  }

  // Sort out alternates.
  if (alt) {
    wcsidx(nwcs, &wcs, alts);

    if (alt == ' ') {
      if (alts[0] == -1) {
        fprintf(stderr, "WARNING, no primary coordinate representation, "
                        "doing all.\n");
        alt = '\0';
      }

    } else if (alt < 'A' || alt > 'Z') {
      fprintf(stderr, "WARNING, alternate specifier \"%c\" is invalid, "
                      "doing all.\n", alt);
      alt = '\0';

    } else {
      if (alts[alt - 'A' + 1] == -1) {
        fprintf(stderr, "WARNING, no alternate coordinate representation "
                        "\"%c\", doing all.\n", alt);
        alt = '\0';
      }
    }
  }

  // Get image dimensions from the header.
  sprintf(keyword, "NAXIS%d", wcs->lng + 1);
  fits_read_key(fptr, TINT, "NAXIS1", naxis,   NULL, &status);
  sprintf(keyword, "NAXIS%d", wcs->lat + 1);
  fits_read_key(fptr, TINT, "NAXIS2", naxis+1, NULL, &status);

  if ((naxis[0] < 2) || (naxis[1] < 2)) {
    fprintf(stderr, "ERROR, HDU number %d contains degenerate image axes.\n",
      hdunum);
    return 1;
  }

  fits_close_file(fptr, &status);


  // Plot setup.
  blc[0] = 0.5f;
  blc[1] = 0.5f;
  trc[0] = naxis[0] + 0.5f;
  trc[1] = naxis[1] + 0.5f;

  if (cpgbeg(0, pgdev, 1, 1) != 1) {
    fprintf(stderr, "ERROR, failed to open PGPLOT device %s.\n", pgdev);
    return 1;
  }
  cpgvstd();

  cpgwnad(blc[0], trc[0], blc[0], trc[1]);
  cpgask(1);
  cpgpage();

  // Compact lettering.
  cpgsch(0.8f);

  // Draw full grid lines.
  gcode[0] = 2;
  gcode[1] = 2;
  grid1[0] =    0.0;
  grid2[0] =    0.0;

  // These are for the projection boundary.
  grid1[1] = -180.0;
  grid1[2] =  180.0;
  grid2[1] =  -90.0;
  grid2[2] =   90.0;

  cpgsci(1);

  for (i = 0; i < nwcs; i++) {
    if (alt && (wcs+i)->alt[0] != alt) {
      continue;
    }

    if ((status = wcsset(wcs+i))) {
      fprintf(stderr, "wcsset ERROR %d: %s.\n", status, wcs_errmsg[status]);
      continue;
    }

    // Draw the frame.
    cpgbox("BC", 0.0f, 0, "BC", 0.0f, 0);

    // Axis labels; use CNAMEia in preference to CTYPEia.
    if ((wcs+i)->cname[0][0]) {
      strcpy(idents[0], (wcs+i)->cname[0]);
    } else {
      strcpy(idents[0], (wcs+i)->ctype[0]);
    }

    if ((wcs+i)->cname[1][0]) {
      strcpy(idents[1], (wcs+i)->cname[1]);
    } else {
      strcpy(idents[1], (wcs+i)->ctype[1]);
    }

    // Title; use WCSNAME.
    strcpy(idents[2], (wcs+i)->wcsname);
    if (strlen(idents[2])) {
      printf("\n%s\n", idents[2]);
    }

    // Formatting control for celestial coordinates.
    if (strncmp((wcs+i)->ctype[0], "RA", 2) == 0) {
      // Right ascension in HMS, declination in DMS.
      opt[0] = 'G';
      opt[1] = 'E';
    } else {
      // Other angles in decimal degrees.
      opt[0] = 'A';
      opt[1] = 'B';
    }

    // Draw the celestial grid.  The grid density is set for each world
    // coordinate by specifying LABDEN = 1224.
    ic = -1;
    cpgsbox(blc, trc, idents, opt, 0, 1224, c0, gcode, 0.0, 0, grid1, 0,
      grid2, 0, pgwcsl_, 1, WCSLEN, 1, nlcprm, (int *)(wcs+i), nldprm, 256,
      &ic, cache, &status);

    // Delimit the projection boundary.
    if ((wcs+i)->cel.prj.category != ZENITHAL) {
      // Reset to the native coordinate graticule.
      (wcs+i)->crval[0] = (wcs+i)->cel.prj.phi0;
      (wcs+i)->crval[1] = (wcs+i)->cel.prj.theta0;
      (wcs+i)->lonpole  = 999.0;
      (wcs+i)->latpole  = 999.0;
      status = wcsset(wcs+i);

      ic = -1;
      cpgsbox(blc, trc, idents, opt, -1, 0, c0, gcode, 0.0, 2, grid1, 2,
        grid2, 0, pgwcsl_, 1, WCSLEN, 1, nlcprm, (int *)(wcs+i), nldprm, 256,
        &ic, cache, &status);
    }

    cpgpage();
  }

  status = wcsvfree(&nwcs, &wcs);

  return 0;

fitserr:
  fits_report_error(stderr, status);
  fits_close_file(fptr, &status);
  return 1;
}