[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 1f8b3ee: All programs: --wcslinearmatrix to se
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 1f8b3ee: All programs: --wcslinearmatrix to select PCi_j and CDi_j in output |
|
Date: |
Mon, 22 Mar 2021 22:24:28 -0400 (EDT) |
branch: master
commit 1f8b3eed8e32e2377c64b1df67e2942647c360f9
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
All programs: --wcslinearmatrix to select PCi_j and CDi_j in output
Until now, Gnuastro's programs would only write their output WCS linear
matrix keywords in the 'PCi_j' formalism (which was the default of
WCSLIB). For the TPV distortion, it would use the 'CDi_j' formalism
(because of a problem in WCSLIB for using 'PCi_j' with the TPV, TNX and ZPX
distortions) until it is fully implemented in WCSLIB and print a warning
about it. However, in a communation with Mark Calabretta (author of
WCSLIB), he informed me that implementing the correction would be very
complicated (in the P.S).
Since the TPV distortion will not work with the 'PCi_j' within WCSLIB any
more, with this commit, the warning has been removed. Furthermore, with
this commit, since the low-level infrastructure already existed to write
the WCS as a 'CDi_j', I generalized it as a common option to all programs:
--wcslinearmatrix. The default mode is still the 'PCi_j', but in some cases
it may be necessary for the users to write the output matrix as 'CDi_j'.
Since this is a pretty low-level operation, the best to implement it was to
use this in the part where we read the WCS structure and internally tweak
the 'wcsprm' so the 'PCi_j' and 'CDi_j' matrices are identical ('CDELTi'
would be 1.0). Then when writing the WCS structure and convert the 'PCi_j'
keywords into 'CDi_j' and delete the 'CDELTi' keywords. Therefore it was
necessary to give the 'wcslinearmatrix' option to the 'gal_wcs_read',
'gal_wcs_read_fitsptr'. Also, to 'gal_wcs_to_cd' function is now a public
function in the library (to let programs also use this option directly: in
particular MakeProfiles).
P.S. Note from WCSLIB's author: "I had also said in a previous email in
that thread that I was working on a fix. In fact, having spent quite a bit
of time on it, I realised that I was creating a monster (without a
compelling reason to do so), and ended up backing out of it completely.
Complications arose in handling TPV in the case where there were additional
axes. The bookkeeping required to keep track of the special handling of
PCi_j + CDELTi for TPV (and TNX and ZPX) bled over into other parts of
WCSLIB that didn't ought to have to deal with such issues, and it made the
code look and feel ugly. (I have often had to complicate the WCSLIB code
to handle past non-standard usage, and it just seems to defeat all of the
effort that went into developing the WCS standard!) I ended up adding a
warning against splitting CDi_j into PCi_j + CDELTi for TPV, (and TNX, and
ZPX)."
---
NEWS | 9 ++++
bin/arithmetic/operands.c | 5 +-
bin/arithmetic/ui.c | 3 +-
bin/convertt/ui.c | 3 +-
bin/convolve/ui.c | 3 +-
bin/crop/ui.c | 9 ++--
bin/fits/fits.c | 6 ++-
bin/fits/keywords.c | 3 +-
bin/fits/ui.c | 1 +
bin/gnuastro.conf | 1 +
bin/mkcatalog/ui.c | 3 +-
bin/mknoise/ui.c | 3 +-
bin/mkprof/ui.c | 7 ++-
bin/noisechisel/ui.c | 3 +-
bin/segment/ui.c | 6 ++-
bin/statistics/statistics.c | 3 +-
bin/statistics/ui.c | 3 +-
bin/table/arithmetic.c | 3 +-
bin/warp/ui.c | 3 +-
doc/gnuastro.texi | 43 +++++++++++++---
lib/gnuastro-internal/commonopts.h | 14 +++++
lib/gnuastro-internal/options.h | 8 ++-
lib/gnuastro/wcs.h | 22 ++++++--
lib/options.c | 49 ++++++++++++++++++
lib/wcs.c | 101 ++++++++++++++++++++-----------------
25 files changed, 232 insertions(+), 82 deletions(-)
diff --git a/NEWS b/NEWS
index f291113..f8a3110 100644
--- a/NEWS
+++ b/NEWS
@@ -22,6 +22,13 @@ See the end of the file for license conditions.
with the profile's value in one column and any requested measure in
the other columns (any MakeCatalog measurement is possible).
+ All programs:
+ --wcslinearmatrix: new option in all programs that let's you select the
+ output WCS linear matrix format. It takes one of two values: 'pc' (for
+ the 'PCi_j' formalism) and 'cd' (for 'CDi_j'). Until now, the outputs
+ were always stored in the 'PCi_j' formalism (which is still the
+ recommended format).
+
Book:
- New "Image surface brightness limit" section added to the third
tutorial (on "Detecting large extended targets"). It describes the
@@ -141,6 +148,8 @@ See the end of the file for license conditions.
desired keyword's value is no longer mandatory. When not given, the
smallest numeric datatype that can keep the value will be found and
used.
+ - gal_wcs_read: allows specifying the linear matrix of the WCS.
+ - gal_wcs_read_fitsptr: allows specifying the linear matrix of the WCS.
** Bugs fixed
bug #60082: Arithmetic library crash for integer operators like modulo
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 1c06500..cb10d3c 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -135,8 +135,9 @@ operands_add(struct arithmeticparams *p, char *filename,
gal_data_t *data)
: NULL);
/* Read the WCS. */
- p->refdata.wcs=gal_wcs_read(filename, newnode->hdu, 0, 0,
- &p->refdata.nwcs);
+ p->refdata.wcs=gal_wcs_read(filename, newnode->hdu,
+ p->cp.wcslinearmatrix,
+ 0, 0, &p->refdata.nwcs);
/* Remove extra (length of 1) dimensions (if we had an
image HDU). */
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 93dcbba..5f13635 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -377,7 +377,8 @@ ui_preparations(struct arithmeticparams *p)
dsize=gal_fits_img_info_dim(p->wcsfile, p->wcshdu, &ndim);
/* Read the WCS. */
- p->refdata.wcs=gal_wcs_read(p->wcsfile, p->wcshdu, 0, 0,
+ p->refdata.wcs=gal_wcs_read(p->wcsfile, p->wcshdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->refdata.nwcs);
if(p->refdata.wcs)
{
diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index d6e306b..48a3d4e 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -528,7 +528,8 @@ ui_make_channels_ll(struct converttparams *p)
/* Read in the array and its WCS information. */
data=gal_fits_img_read(name->v, hdu, p->cp.minmapsize,
p->cp.quietmmap);
- data->wcs=gal_wcs_read(name->v, hdu, 0, 0, &data->nwcs);
+ data->wcs=gal_wcs_read(name->v, hdu, p->cp.wcslinearmatrix,
+ 0, 0, &data->nwcs);
data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize,
data->wcs);
gal_list_data_add(&p->chll, data);
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index e466975..2c2582f 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -438,7 +438,8 @@ ui_read_input(struct convolveparams *p)
INPUT_USE_TYPE,
p->cp.minmapsize,
p->cp.quietmmap);
- p->input->wcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0,
+ p->input->wcs=gal_wcs_read(p->filename, p->cp.hdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
p->input->dsize,
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index a1bc980..8776f15 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -855,8 +855,9 @@ ui_preparations_to_img_mode(struct cropparams *p)
size_t i;
int nwcs;
double *darr, *pixscale;
- struct wcsprm *wcs=gal_wcs_read(p->inputs->v, p->cp.hdu, p->hstartwcs,
- p->hendwcs, &nwcs);
+ struct wcsprm *wcs=gal_wcs_read(p->inputs->v, p->cp.hdu,
+ p->cp.wcslinearmatrix,
+ p->hstartwcs, p->hendwcs, &nwcs);
/* Make sure a WCS actually exists. */
if(wcs==NULL)
@@ -1018,8 +1019,8 @@ ui_preparations(struct cropparams *p)
tmpfits=gal_fits_hdu_open_format(img->name, p->cp.hdu, 0);
gal_fits_img_info(tmpfits, &p->type, &img->ndim, &img->dsize,
NULL, NULL);
- img->wcs=gal_wcs_read_fitsptr(tmpfits, p->hstartwcs, p->hendwcs,
- &img->nwcs);
+ img->wcs=gal_wcs_read_fitsptr(tmpfits, p->cp.wcslinearmatrix,
+ p->hstartwcs, p->hendwcs, &img->nwcs);
if(img->wcs)
{
gal_wcs_decompose_pc_cdelt(img->wcs);
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 19a7a1b..718d50c 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -323,7 +323,8 @@ fits_pixelscale(struct fitsparams *p)
double multip, *pixelscale;
/* Read the desired WCS. */
- wcs=gal_wcs_read(p->input->v, p->cp.hdu, 0, 0, &nwcs);
+ wcs=gal_wcs_read(p->input->v, p->cp.hdu, p->cp.wcslinearmatrix,
+ 0, 0, &nwcs);
/* If a WCS doesn't exist, let the user know and return. */
if(wcs)
@@ -474,7 +475,8 @@ fits_skycoverage(struct fitsparams *p)
}
/* For the range type of coverage. */
- wcs=gal_wcs_read(p->input->v, p->cp.hdu, 0, 0, &nwcs);
+ wcs=gal_wcs_read(p->input->v, p->cp.hdu, p->cp.wcslinearmatrix,
+ 0, 0, &nwcs);
printf("\nSky coverage by range along dimensions:\n");
for(i=0;i<ndim;++i)
printf(" %-8s %-15.10g%-15.10g\n", gal_wcs_dimension_name(wcs, i),
diff --git a/bin/fits/keywords.c b/bin/fits/keywords.c
index 9a325f6..9b48fd7 100644
--- a/bin/fits/keywords.c
+++ b/bin/fits/keywords.c
@@ -463,7 +463,8 @@ keywords_distortion_wcs(struct fitsparams *p)
}
/* Read the input's WCS and make sure one exists. */
- inwcs=gal_wcs_read(p->input->v, p->cp.hdu, 0, 0, &nwcs);
+ inwcs=gal_wcs_read(p->input->v, p->cp.hdu, p->cp.wcslinearmatrix,
+ 0, 0, &nwcs);
if(inwcs==NULL)
error(EXIT_FAILURE, 0, "%s (hdu %s): doesn't have any WCS structure "
"for converting its distortion",
diff --git a/bin/fits/ui.c b/bin/fits/ui.c
index c9cbc63..fb741fb 100644
--- a/bin/fits/ui.c
+++ b/bin/fits/ui.c
@@ -120,6 +120,7 @@ ui_initialize_options(struct fitsparams *p,
case GAL_OPTIONS_KEY_SEARCHIN:
case GAL_OPTIONS_KEY_IGNORECASE:
case GAL_OPTIONS_KEY_TYPE:
+ case GAL_OPTIONS_KEY_WCSLINEARMATRIX:
case GAL_OPTIONS_KEY_DONTDELETE:
case GAL_OPTIONS_KEY_LOG:
case GAL_OPTIONS_KEY_NUMTHREADS:
diff --git a/bin/gnuastro.conf b/bin/gnuastro.conf
index 6c180c0..a93d7cb 100644
--- a/bin/gnuastro.conf
+++ b/bin/gnuastro.conf
@@ -36,6 +36,7 @@
# Output:
tableformat fits-binary
+ wcslinearmatrix pc
# Operating mode
quietmmap 0
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index d84c9c3..aa88388 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -538,7 +538,8 @@ ui_wcs_info(struct mkcatalogparams *p)
size_t i;
/* Read the WCS meta-data. */
- p->objects->wcs=gal_wcs_read(p->objectsfile, p->cp.hdu, 0, 0,
+ p->objects->wcs=gal_wcs_read(p->objectsfile, p->cp.hdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->objects->nwcs);
/* Read the basic WCS information. */
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index f39bbc9..4c99355 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -311,7 +311,8 @@ ui_preparations(struct mknoiseparams *p)
p->input=gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu, NULL,
GAL_TYPE_FLOAT64, p->cp.minmapsize,
p->cp.quietmmap);
- p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs);
+ p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, p->cp.wcslinearmatrix,
+ 0, 0, &p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim, p->input->dsize,
p->input->wcs);
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 5e1ad36..9b29674 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -1359,6 +1359,10 @@ ui_prepare_wcs(struct mkprofparams *p)
if(status)
error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
wcs_errmsg[status]);
+
+ /* Convert it to CD if the user wanted it. */
+ if(p->cp.wcslinearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
+ gal_wcs_to_cd(wcs);
}
@@ -1384,7 +1388,8 @@ ui_prepare_canvas(struct mkprofparams *p)
the background image and the number of its dimensions. So
'ndim==0' and what 'dsize' points to is irrelevant. */
tdsize=gal_fits_img_info_dim(p->backname, p->backhdu, &tndim);
- p->wcs=gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->nwcs);
+ p->wcs=gal_wcs_read(p->backname, p->backhdu, p->cp.wcslinearmatrix,
+ 0, 0, &p->nwcs);
tndim=gal_dimension_remove_extra(tndim, tdsize, p->wcs);
free(tdsize);
if(p->nomerged==0)
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index df71670..403f5c9 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -567,7 +567,8 @@ ui_preparations_read_input(struct noisechiselparams *p)
NULL, GAL_TYPE_FLOAT32,
p->cp.minmapsize,
p->cp.quietmmap);
- p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu, 0, 0,
+ p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
p->input->dsize,
diff --git a/bin/segment/ui.c b/bin/segment/ui.c
index d65b103..93e3ed8 100644
--- a/bin/segment/ui.c
+++ b/bin/segment/ui.c
@@ -411,8 +411,10 @@ ui_prepare_inputs(struct segmentparams *p)
/* Read the input as a single precision floating point dataset. */
p->input = gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
NULL, GAL_TYPE_FLOAT32,
- p->cp.minmapsize, p->cp.quietmmap);
- p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu, 0, 0,
+ p->cp.minmapsize,
+ p->cp.quietmmap);
+ p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
p->input->dsize,
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 48986eb..afacdc0 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -832,7 +832,8 @@ histogram_2d(struct statisticsparams *p)
cunit[0] = p->input->unit; cunit[1] = p->input->next->unit;
ctype[0] = histogram_2d_set_ctype(p->input->name, "X");
ctype[1] = histogram_2d_set_ctype(p->input->next->name, "Y");
- img->wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit, ctype, 2);
+ img->wcs=gal_wcs_create(crpix, crval, cdelt, pc, cunit, ctype, 2,
+ p->cp.wcslinearmatrix);
/* Write the output. */
output=statistics_output_name(p, suf, &isfits);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 0b2ad5c..7c2b7fc 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -971,7 +971,8 @@ ui_preparations(struct statisticsparams *p)
p->inputformat=INPUT_FORMAT_IMAGE;
p->input=gal_array_read_one_ch(p->inputname, cp->hdu, NULL,
cp->minmapsize, p->cp.quietmmap);
- p->input->wcs=gal_wcs_read(p->inputname, cp->hdu, 0, 0,
+ p->input->wcs=gal_wcs_read(p->inputname, cp->hdu,
+ p->cp.wcslinearmatrix, 0, 0,
&p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
p->input->dsize,
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index dfa576c..b964535 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -159,7 +159,8 @@ arithmetic_init_wcs(struct tableparams *p, char *operator)
"for the '%s' operator", operator);
/* Read the WCS. */
- p->wcs=gal_wcs_read(p->wcsfile, p->wcshdu, 0, 0, &p->nwcs);
+ p->wcs=gal_wcs_read(p->wcsfile, p->wcshdu, p->cp.wcslinearmatrix,
+ 0, 0, &p->nwcs);
if(p->wcs==NULL)
error(EXIT_FAILURE, 0, "%s (hdu: %s): no WCS could be read by "
"WCSLIB", p->wcsfile, p->wcshdu);
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index a91fb14..72659c1 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -352,7 +352,8 @@ ui_check_options_and_arguments(struct warpparams *p)
p->cp.quietmmap);
/* Read the WCS and remove one-element wide dimension(s). */
- p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, p->hstartwcs,
+ p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu,
+ p->cp.wcslinearmatrix, p->hstartwcs,
p->hendwcs, &p->input->nwcs);
p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
p->input->dsize,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4161e29..ec5d28f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -7115,6 +7115,17 @@ But there are two types of FITS tables: FITS ASCII, and
FITS binary.
Thus, with this option, the program is able to identify which type you want.
The currently recognized values to this option are:
+@item --wcslinearmatrix=STR
+Select the linear transformation matrix of the output's WCS.
+This option only takes two values: @code{pc} (for the @code{PCi_j} formalism)
and @code{cd} (for @code{CDi_j}).
+For more on the different formalisms, please see Section 8.1 of the FITS
standard@footnote{@url{https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf}},
version 4.0.
+
+@cindex @code{CDELT}
+In short, in the @code{PCi_j} formalism, we only keep the linear rotation
matrix in these keywords and put the scaling factor (or the pixel scale in
astronomical imaging) in the @code{CDELTi} keywords.
+In the @code{CDi_j} formalism, we blend the scaling into the rotation into a
single matrix and keep that matrix in these FITS keywords.
+By default, Gnuastro uses the @code{PCi_j} formalism, because it greatly helps
in human readability of the raw keywords and is also the default mode of WCSLIB.
+However, in some circumstances it may be necessary to have the keywords in the
CD format; for example when you need to feed the outputs into other software
that don't follow the full FITS standard and only recognize the @code{CDi_j}
formalism.
+
@table @command
@item txt
A plain text table with white-space characters between the columns (see
@@ -25595,22 +25606,32 @@ TPD is a superset of all these, hence it has both
prior and sequeal distortion c
More information is given in the documentation of @code{dis.h}, from the
WCSLIB
manual@footnote{@url{https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/dis_8h.html}}.
@end deffn
+@deffn Macro GAL_WCS_LINEAR_MATRIX_PC
+@deffnx Macro GAL_WCS_LINEAR_MATRIX_CD
+Identifiers of the linear transformation matrix: either in the @code{PCi_j} or
the @code{CDi_j} formalism.
+For more, see the description of @option{--wcslinearmatrix} in @ref{Input
output options}.
+@end deffn
+
+
@deffn Macro GAL_WCS_FLTERROR
Limit of rounding for floating point errors.
@end deffn
-@deftypefun {struct wcsprm *} gal_wcs_read_fitsptr (fitsfile @code{*fptr},
size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
-[@strong{Not thread-safe}] Return the WCSLIB @code{wcsprm} structure that
-is read from the CFITSIO @code{fptr} pointer to an opened FITS file. Also
-put the number of coordinate representations found into the space that
-@code{nwcs} points to. To read the WCS structure directly from a filename,
-see @code{gal_wcs_read} below. After processing has finished, you can free
-the returned structure with WCSLIB's @code{wcsvfree} keyword:
+@deftypefun {struct wcsprm *} gal_wcs_read_fitsptr (fitsfile @code{*fptr}, int
@code{linearmatrix}, size_t @code{hstartwcs}, size_t @code{hendwcs}, int
@code{*nwcs})
+[@strong{Not thread-safe}] Return the WCSLIB @code{wcsprm} structure that is
read from the CFITSIO @code{fptr} pointer to an opened FITS file.
+Also put the number of coordinate representations found into the space that
@code{nwcs} points to.
+To read the WCS structure directly from a filename, see @code{gal_wcs_read}
below.
+After processing has finished, you can free the returned structure with
WCSLIB's @code{wcsvfree} keyword:
@example
status = wcsvfree(&nwcs,&wcs);
@end example
+The @code{linearmatrix} argument takes one of three values: @code{0},
@code{GAL_WCS_LINEAR_MATRIX_PC} and @code{GAL_WCS_LINEAR_MATRIX_CD}.
+It will determine the format of the WCS when it is later written to file with
@code{gal_wcs_write} or @code{gal_wcs_write_in_fitsptr} (which is called by
@code{gal_fits_img_write})
+So if you don't want to write the WCS into a file later, just give it a value
of @code{0}.
+For more on the difference between these modes, see the description of
@option{--wcslinearmatrix} in @ref{Input output options}.
+
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 @code{hstartwcs} and @code{hendwcs} arguments
to specify the keyword number range (counting from zero).
If @code{hendwcs} is larger than @code{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.
@@ -25622,7 +25643,7 @@ This function is just a wrapper over WCSLIB's
@code{wcspih} function which is no
Therefore, be sure to not call this function simultaneously (over multiple
threads).
@end deftypefun
-@deftypefun {struct wcsprm *} gal_wcs_read (char @code{*filename}, char
@code{*hdu}, size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
+@deftypefun {struct wcsprm *} gal_wcs_read (char @code{*filename}, char
@code{*hdu}, int @code{linearmatrix}, size_t @code{hstartwcs}, size_t
@code{hendwcs}, int @code{*nwcs})
[@strong{Not thread-safe}] Return the WCSLIB structure that is read from the
HDU/extension @code{hdu} of the file @code{filename}.
Also put the number of coordinate representations found into the space that
@code{nwcs} points to.
Please see @code{gal_wcs_read_fitsptr} for more.
@@ -25708,6 +25729,12 @@ correspond to the pixel scale, and the @code{PCi_j}
will correction show
the rotation.
@end deftypefun
+@deftypefun void gal_wcs_to_cd (struct wcsprm @code{*wcs})
+Make sure that the WCS structure's @code{PCi_j} and @code{CDi_j} keywords have
the same value and that the @code{CDELTi} keywords have a value of 1.0.
+Also, set the @code{wcs->altlin=2} (for the @code{CDi_j} formalism).
+With these changes @code{gal_wcs_write_in_fitsptr} (and thus
@code{gal_wcs_write} and @code{gal_fits_img_write} and its derivates) will have
an output file in the format of @code{CDi_j}.
+@end deftypefun
+
@deftypefun int gal_wcs_distortion_from_string (char @code{*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 @code{GAL_WCS_DISTORTION_*} macros defined above).
The sting-based distortion identifiers have three characters and are all in
capital letters.
diff --git a/lib/gnuastro-internal/commonopts.h
b/lib/gnuastro-internal/commonopts.h
index d55f98b..bec62fa 100644
--- a/lib/gnuastro-internal/commonopts.h
+++ b/lib/gnuastro-internal/commonopts.h
@@ -276,6 +276,20 @@ struct argp_option gal_commonopts_options[] =
gal_options_read_tableformat
},
{
+ "wcslinearmatrix",
+ GAL_OPTIONS_KEY_WCSLINEARMATRIX,
+ "STR",
+ 0,
+ "WCS linear matrix of output ('pc' or 'cd').",
+ GAL_OPTIONS_GROUP_OUTPUT,
+ &cp->wcslinearmatrix,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_wcslinearmatrix
+ },
+ {
"dontdelete",
GAL_OPTIONS_KEY_DONTDELETE,
0,
diff --git a/lib/gnuastro-internal/options.h b/lib/gnuastro-internal/options.h
index fdb7478..343ebf8 100644
--- a/lib/gnuastro-internal/options.h
+++ b/lib/gnuastro-internal/options.h
@@ -125,6 +125,7 @@ enum options_common_keys
GAL_OPTIONS_KEY_INTERPONLYBLANK,
GAL_OPTIONS_KEY_INTERPMETRIC,
GAL_OPTIONS_KEY_INTERPNUMNGB,
+ GAL_OPTIONS_KEY_WCSLINEARMATRIX,
};
@@ -194,9 +195,10 @@ struct gal_options_common_params
/* Output. */
char *output; /* Directory containg output. */
uint8_t type; /* Data type of output. */
+ uint8_t tableformat; /* Internal code for output table format. */
+ uint8_t wcslinearmatrix; /* WCS matrix to use (PC or CD). */
uint8_t dontdelete; /* ==1: Don't delete existing file. */
uint8_t keepinputdir; /* Keep input directory for auto output. */
- uint8_t tableformat; /* Internal code for output table format. */
/* Operating modes. */
uint8_t quiet; /* Only print errors. */
@@ -276,6 +278,10 @@ gal_options_read_searchin(struct argp_option *option, char
*arg,
char *filename, size_t lineno, void *junk);
void *
+gal_options_read_wcslinearmatrix(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk);
+
+void *
gal_options_read_tableformat(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 8ad29ba..beef0bc 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -70,6 +70,15 @@ enum gal_wcs_distortions
GAL_WCS_DISTORTION_WAT, /* The WAT polynomial distortion. */
};
+/* Macros to identify the type of distortion for conversions. */
+enum gal_wcs_linear_matrix
+{
+ GAL_WCS_MATRIX_INVALID, /* Invalid (=0 by C standard). */
+
+ GAL_WCS_LINEAR_MATRIX_PC,
+ GAL_WCS_LINEAR_MATRIX_CD,
+};
+
@@ -78,16 +87,17 @@ enum gal_wcs_distortions
*********** Read WCS ***********
*************************************************************/
struct wcsprm *
-gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t hendwcs,
- int *nwcs);
+gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, size_t hstartwcs,
+ size_t hendwcs, int *nwcs);
struct wcsprm *
-gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
+gal_wcs_read(char *filename, char *hdu, int linearmatrix, size_t hstartwcs,
size_t hendwcs, int *nwcs);
struct wcsprm *
gal_wcs_create(double *crpix, double *crval, double *cdelt,
- double *pc, char **cunit, char **ctype, size_t ndim);
+ double *pc, char **cunit, char **ctype, size_t ndim,
+ int linearmatrix);
char *
gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension);
@@ -107,7 +117,6 @@ gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm
*wcs);
-
/*************************************************************
*********** Distortions ***********
*************************************************************/
@@ -149,6 +158,9 @@ gal_wcs_clean_errors(struct wcsprm *wcs);
void
gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs);
+void
+gal_wcs_to_cd(struct wcsprm *wcs);
+
double
gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
diff --git a/lib/options.c b/lib/options.c
index e029072..73f33e5 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <string.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/git.h>
#include <gnuastro/txt.h>
#include <gnuastro/fits.h>
@@ -477,6 +478,54 @@ gal_options_read_searchin(struct argp_option *option, char
*arg,
void *
+gal_options_read_wcslinearmatrix(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk)
+{
+ char *str;
+ uint8_t value;
+ if(lineno==-1)
+ {
+ /* The output must be an allocated string (will be 'free'd later). */
+ value=*(uint8_t *)(option->value);
+ switch(value)
+ {
+ case GAL_WCS_LINEAR_MATRIX_PC: gal_checkset_allocate_copy("pc", &str);
+ break;
+ case GAL_WCS_LINEAR_MATRIX_CD: gal_checkset_allocate_copy("cd", &str);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+ "to fix the problem. %u is not a recognized WCS rotation "
+ "matrix code", __func__, PACKAGE_BUGREPORT, value);
+ }
+ return str;
+ }
+ else
+ {
+ /* If the option is already set, just return. */
+ if(option->set) return NULL;
+
+ /* Read the value. */
+ if( !strcmp(arg, "pc") ) value = GAL_WCS_LINEAR_MATRIX_PC;
+ else if( !strcmp(arg, "cd") ) value = GAL_WCS_LINEAR_MATRIX_CD;
+ else
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "'%s' (value "
+ "to '%s' option) couldn't be recognized as a known "
+ "WCS rotation matrix. Acceptable values are 'pc' "
+ "or 'cd'", arg, option->name);
+ *(uint8_t *)(option->value)=value;
+
+ /* For no un-used variable warning. This function doesn't need the
+ pointer.*/
+ return junk=NULL;
+ }
+}
+
+
+
+
+
+void *
gal_options_read_tableformat(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *junk)
{
diff --git a/lib/wcs.c b/lib/wcs.c
index f1d8144..40436de 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -53,14 +53,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/* Static functions on for this file. */
-static void
-gal_wcs_to_cd(struct wcsprm *wcs);
-
-
-
-
-
/*************************************************************
*********** Read WCS ***********
*************************************************************/
@@ -149,8 +141,8 @@ wcs_read_correct_pc_cd(struct wcsprm *wcs)
Don't call this function within a thread or use a mutex.
*/
struct wcsprm *
-gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t hendwcs,
- int *nwcs)
+gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, size_t hstartwcs,
+ size_t hendwcs, int *nwcs)
{
/* Declaratins: */
int sumcheck;
@@ -361,6 +353,9 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs,
size_t hendwcs,
}
}
+ /* If the user wants a CD linear matrix, do the conversion here. */
+ if(linearmatrix==GAL_WCS_LINEAR_MATRIX_CD) gal_wcs_to_cd(wcs);
+
/* Clean up and return. */
status=0;
if (fits_free_memory(fullheader, &status) )
@@ -374,8 +369,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs,
size_t hendwcs,
struct wcsprm *
-gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
- size_t hendwcs, int *nwcs)
+gal_wcs_read(char *filename, char *hdu, int linearmatrix,
+ size_t hstartwcs, size_t hendwcs, int *nwcs)
{
int status=0;
fitsfile *fptr;
@@ -389,7 +384,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
fptr=gal_fits_hdu_open_format(filename, hdu, 0);
/* Read the WCS information: */
- wcs=gal_wcs_read_fitsptr(fptr, hstartwcs, hendwcs, nwcs);
+ wcs=gal_wcs_read_fitsptr(fptr, linearmatrix, hstartwcs,
+ hendwcs, nwcs);
/* Close the FITS file and return. */
fits_close_file(fptr, &status);
@@ -403,7 +399,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
struct wcsprm *
gal_wcs_create(double *crpix, double *crval, double *cdelt,
- double *pc, char **cunit, char **ctype, size_t ndim)
+ double *pc, char **cunit, char **ctype,
+ size_t ndim, int linearmatrix)
{
size_t i;
int status;
@@ -441,6 +438,10 @@ gal_wcs_create(double *crpix, double *crval, double *cdelt,
error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
wcs_errmsg[status]);
+ /* If a CD matrix is desired make it. */
+ if(linearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
+ gal_wcs_to_cd(wcs);
+
/* Return the output WCS. */
return wcs;
}
@@ -497,15 +498,21 @@ void
gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
{
char *wcsstr;
- int tpvdist, status=0, nkeyrec;
-
- /* Prepare the main rotation matrix. Note that for TPV distortion, WCSLIB
- versions 7.3 and before couldn't deal with the CDELT keys, so to be
- safe, in such cases, we'll remove the effect of CDELT in the
- 'gal_wcs_to_cd' function. */
- tpvdist=wcs->lin.disseq && !strcmp(wcs->lin.disseq->dtype[1], "TPV");
- if( tpvdist ) gal_wcs_to_cd(wcs);
- else gal_wcs_decompose_pc_cdelt(wcs);
+ int cdfordist, status=0, nkeyrec;
+
+ /* For the TPV, TNX and ZPX distortions, WCSLIB can't deal with the CDELT
+ keys properly and its better to use the CD matrix instead, so we'll
+ use the 'gal_wcs_to_cd' function. */
+ cdfordist = ( wcs->lin.disseq
+ && ( !strcmp( wcs->lin.disseq->dtype[1], "TPV")
+ || !strcmp(wcs->lin.disseq->dtype[1], "TNX")
+ || !strcmp(wcs->lin.disseq->dtype[1], "ZPX") ) );
+
+ /* Finalize the linear transformation matrix. Note that some programs may
+ have worked on the WCS. So even if 'altlin' is already 2, we'll just
+ ensure that the final matrix is CD here. */
+ if(wcs->altlin==2 || cdfordist) gal_wcs_to_cd(wcs);
+ else gal_wcs_decompose_pc_cdelt(wcs);
/* Clean up small errors in the PC matrix and CDELT values. */
gal_wcs_clean_errors(wcs);
@@ -523,33 +530,33 @@ gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm
*wcs)
status=0;
/* WCSLIB is going to write PC+CDELT keywords in any case. But when we
- have a TPV distortion, it is cleaner to use a CD matrix. Also,
- including and before version 7.3, WCSLIB wouldn't convert coordinates
- properly if the PC matrix is used with the TPV distortion. So to help
- users with WCSLIB 7.3 or earlier, we need to conver the PC matrix to
- CD. 'gal_wcs_to_cd' function already made sure that CDELT=1, so
- effectively the CD matrix and PC matrix are equivalent, we just need
- to convert the keyword names and delete the CDELT keywords. Note that
- zero-valued PC/CD elements may not be present, so we'll manually set
- 'status' to zero and not let CFITSIO crash.*/
+ have a TPV, TNX or ZPX distortion, it is cleaner to use a CD matrix
+ (WCSLIB can't convert coordinates properly if the PC matrix is used
+ with the TPV distortion). So to help users avoid the potential
+ problems with WCSLIB. 'gal_wcs_to_cd' function already made sure that
+ CDELTi=1.0, so effectively the CD matrix and PC matrix are
+ equivalent, we just need to convert the keyword names and delete the
+ CDELT keywords. Note that zero-valued PC/CD elements may not be
+ present, so we'll manually set 'status' to zero to avoid CFITSIO from
+ crashing. */
if(wcs->altlin==2)
{
+ status=0; fits_delete_str(fptr, "CDELT1", &status);
+ status=0; fits_delete_str(fptr, "CDELT2", &status);
status=0; fits_modify_name(fptr, "PC1_1", "CD1_1", &status);
status=0; fits_modify_name(fptr, "PC1_2", "CD1_2", &status);
status=0; fits_modify_name(fptr, "PC2_1", "CD2_1", &status);
status=0; fits_modify_name(fptr, "PC2_2", "CD2_2", &status);
- status=0; fits_delete_str(fptr, "CDELT1", &status);
- status=0; fits_delete_str(fptr, "CDELT2", &status);
+ if(wcs->naxis==3)
+ {
+ status=0; fits_delete_str(fptr, "CDELT3", &status);
+ status=0; fits_modify_name(fptr, "PC1_3", "CD1_3", &status);
+ status=0; fits_modify_name(fptr, "PC2_3", "CD2_3", &status);
+ status=0; fits_modify_name(fptr, "PC3_1", "CD3_1", &status);
+ status=0; fits_modify_name(fptr, "PC3_2", "CD3_2", &status);
+ status=0; fits_modify_name(fptr, "PC3_3", "CD3_3", &status);
+ }
status=0;
- fits_write_comment(fptr, "The CD matrix is used instead of the "
- "PC+CDELT due to conflicts with TPV distortion "
- "in WCSLIB 7.3 (released on 2020/06/03) and "
- "ealier. By default Gnuastro will write "
- "PC+CDELT matrices because the rotation (PC) and "
- "pixel-scale (CDELT) are separate; providing "
- "more physically relevant metadata for human "
- "readers (PC+CDELT is also the default format "
- "of WCSLIB).", &status);
}
}
@@ -1307,7 +1314,7 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
/* Set the WCS structure to use the CD matrix. */
-static void
+void
gal_wcs_to_cd(struct wcsprm *wcs)
{
size_t i, j, naxis;
@@ -1588,8 +1595,10 @@ gal_wcs_coverage(char *filename, char *hdu, size_t
*ondim,
size_t i, ndim, *dsize=NULL, numrows;
double *x=NULL, *y=NULL, *z=NULL, *min, *max, *center, *width;
- /* Read the desired WCS. */
- wcs=gal_wcs_read(filename, hdu, 0, 0, &nwcs);
+ /* Read the desired WCS (note that the linear matrix is irrelevant here,
+ we'll just select PC because its the default WCS mode. */
+ wcs=gal_wcs_read(filename, hdu, GAL_WCS_LINEAR_MATRIX_PC,
+ 0, 0, &nwcs);
/* If a WCS doesn't exist, return NULL. */
if(wcs==NULL) return 0;
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 1f8b3ee: All programs: --wcslinearmatrix to select PCi_j and CDi_j in output,
Mohammad Akhlaghi <=