[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 20f1d41: Fits: new --wcscoordsys option to con
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 20f1d41: Fits: new --wcscoordsys option to convert coordinate systems |
|
Date: |
Thu, 25 Mar 2021 21:15:18 -0400 (EDT) |
branch: master
commit 20f1d41d7ce559753ff98736ce0543b6615e4f4a
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Fits: new --wcscoordsys option to convert coordinate systems
Until now there was no easy way to convert the coordinates of an image's
WCS. I had previously asked about this from Mark Calabretta (author of
WCSLIB) about this and he kindly added the 'wcsccs' function for this job
to WCSLIB 7.5 (released on March 21st, 2021).
With this commit, the 'gal_wcs_coordsys_convert' wrapper function (over
'wcsccs') has been added to Gnuastro's library for this job (besides a
function to read the user's desired coordinate system and another function
to set the pole coordinates for the various conversions). The pole
coordinates were found (inspired from the examples in WCSLIB's manual) from
the NED coordinate calculator.
With the library functions in place, it was very easy to add the necessary
steps in the Fits program for the new '--wcscoordsys' option. It will parse
the input WCS structure, find its coordinate system, and convert it to the
user's desired system.
---
NEWS | 20 +++
bin/fits/args.h | 13 ++
bin/fits/keywords.c | 22 ++-
bin/fits/main.h | 2 +
bin/fits/ui.c | 13 +-
bin/fits/ui.h | 1 +
configure.ac | 17 +++
doc/gnuastro.texi | 76 +++++++++-
lib/Makefile.am | 1 +
lib/gnuastro/wcs.h | 27 ++++
lib/wcs.c | 400 ++++++++++++++++++++++++++++++++++++++++++++++++++++
11 files changed, 580 insertions(+), 12 deletions(-)
diff --git a/NEWS b/NEWS
index f8a3110..dc4ac45 100644
--- a/NEWS
+++ b/NEWS
@@ -53,6 +53,16 @@ See the end of the file for license conditions.
- counts-to-jy: Convert counts to Janskys through a zero point based
on AB magnitudes.
+ Fits:
+ --wcscoordsys: convert the WCS coordinate system of the input into any
+ recognized coordinate system (currently supports: equatorial (J2000,
+ B1950), ecliptic (J2000, B1950), Galactic and Supergalactic. For
+ example if 'image.fits' is in galactic coordinates, you can use this
+ command to convert its WCS to equatorial (J2000):
+ astfits image.fits --wcscoordsys=eq-j2000
+ This option only works with WCSLIB 7.5 and above (released in March
+ 2021), otherwise it will not work (abort with an informative warning).
+
MakeCatalog:
- Newly added measurement columns:
--upperlimitsb: upper-limit surface brightness for the given label
@@ -117,6 +127,16 @@ See the end of the file for license conditions.
- GAL_ARITHMETIC_OP_ACOSH: Inverse hyperbolic cosine.
- GAL_ARITHMETIC_OP_ATANH: Inverse hyperbolic tangent.
- GAL_ARITHMETIC_OP_COUNTS_TO_JY: Convert counts to Janskys.
+ - WCS coordinate system identifiers:
+ - GAL_WCS_COORDSYS_EQB1950: 1950.0 (Besselian-year) equatorial coords.
+ - GAL_WCS_COORDSYS_EQJ2000: 2000.0 (Julian-year) equatorial coords.
+ - GAL_WCS_COORDSYS_ECB1950: 1950.0 (Besselian-year) ecliptic coords.
+ - GAL_WCS_COORDSYS_ECJ2000: 2000.0 (Julian-year) ecliptic coords.
+ - GAL_WCS_COORDSYS_GALACTIC: Galactic coordinates.
+ - GAL_WCS_COORDSYS_SUPERGALACTIC: Supergalactic coordinates.
+ - gal_wcs_coordsys_from_string: WCS coordinate system from string.
+ - gal_wcs_coordsys_identify: Parse WCS struct to find coordinate system.
+ - gal_wcs_coordsys_convert: Convert the coordinate system of the WCS.
** Removed features
diff --git a/bin/fits/args.h b/bin/fits/args.h
index af1d90c..ffe56a9 100644
--- a/bin/fits/args.h
+++ b/bin/fits/args.h
@@ -341,6 +341,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "wcscoordsys",
+ UI_KEY_WCSCOORDSYS,
+ "STR",
+ 0,
+ "Convert WCS coordinate system.",
+ UI_GROUP_KEYWORD,
+ &p->wcscoordsys,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
diff --git a/bin/fits/keywords.c b/bin/fits/keywords.c
index 852858f..4a1d483 100644
--- a/bin/fits/keywords.c
+++ b/bin/fits/keywords.c
@@ -438,7 +438,7 @@ keywords_date_to_seconds(struct fitsparams *p, fitsfile
*fptr)
static void
-keywords_distortion_wcs(struct fitsparams *p)
+keywords_wcs_convert(struct fitsparams *p)
{
int nwcs;
size_t ndim, *insize;
@@ -467,13 +467,13 @@ keywords_distortion_wcs(struct fitsparams *p)
0, 0, &nwcs);
if(inwcs==NULL)
error(EXIT_FAILURE, 0, "%s (hdu %s): doesn't have any WCS structure "
- "for converting its distortion",
+ "for converting its coordinate system or distortion",
p->input->v, p->cp.hdu);
/* In case there is no dataset and the conversion is between TPV to SIP,
we need to set a default size and use that for the conversion, but we
also need to warn the user. */
- if(data==NULL)
+ if(p->wcsdistortion && data==NULL)
{
if( !p->cp.quiet
&& gal_wcs_distortion_identify(inwcs)==GAL_WCS_DISTORTION_TPV
@@ -492,14 +492,22 @@ keywords_distortion_wcs(struct fitsparams *p)
else dsize=data->dsize;
/* Do the conversion. */
- outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid, dsize);
+ if(p->wcscoordsys)
+ outwcs=gal_wcs_coordsys_convert(inwcs, p->coordsysid);
+ else if(p->wcsdistortion)
+ outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid, dsize);
+ else
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The requested mode for this function is not "
+ "recognized", __func__, PACKAGE_BUGREPORT);
/* Set the output filename. */
if(p->cp.output)
output=p->cp.output;
else
{
- if( asprintf(&suffix, "-%s.fits", p->wcsdistortion)<0 )
+ if( asprintf(&suffix, "-%s.fits",
+ p->wcsdistortion ? p->wcsdistortion : p->wcscoordsys)<0 )
error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
output=gal_checkset_automatic_output(&p->cp, p->input->v, suffix);
}
@@ -1035,8 +1043,8 @@ keywords(struct fitsparams *p)
}
/* Convert the input's distortion to the desired output distortion. */
- if(p->wcsdistortion)
- keywords_distortion_wcs(p);
+ if(p->wcsdistortion || p->wcscoordsys)
+ keywords_wcs_convert(p);
/* Return. */
return r;
diff --git a/bin/fits/main.h b/bin/fits/main.h
index 9c66221..c0254c2 100644
--- a/bin/fits/main.h
+++ b/bin/fits/main.h
@@ -78,12 +78,14 @@ struct fitsparams
uint8_t *verify; /* Verify the CHECKSUM and DATASUM keys. */
char *copykeys; /* Range of keywords to copy in output. */
char *datetosec; /* Convert FITS date to seconds. */
+ char *wcscoordsys; /* Name of new WCS coordinate system. */
char *wcsdistortion; /* WCS distortion to write in output. */
uint8_t quitonerror; /* Quit if an error occurs. */
uint8_t colinfoinstdout; /* Print column info in output. */
/* Internal: */
int mode; /* Operating on HDUs or keywords. */
+ int coordsysid; /* ID of desired coordinate system.*/
int distortionid; /* ID of desired distortion. */
long copykeysrange[2]; /* Start and end of copy. */
gal_fits_list_key_t *write_keys; /* Keys to write in the header. */
diff --git a/bin/fits/ui.c b/bin/fits/ui.c
index fb741fb..b8e9436 100644
--- a/bin/fits/ui.c
+++ b/bin/fits/ui.c
@@ -325,7 +325,7 @@ ui_read_check_only_options(struct fitsparams *p)
if( p->date || p->comment || p->history || p->asis || p->keyvalue
|| p->delete || p->rename || p->update || p->write || p->verify
|| p->printallkeys || p->copykeys || p->datetosec
- || p->wcsdistortion )
+ || p->wcscoordsys || p->wcsdistortion )
{
/* Check if a HDU is given. */
if(p->cp.hdu==NULL)
@@ -340,15 +340,20 @@ ui_read_check_only_options(struct fitsparams *p)
/* Keyword-related options that must be called alone. */
checkkeys = ( (p->keyvalue!=NULL)
+ (p->datetosec!=NULL)
+ + (p->wcscoordsys!=NULL)
+ (p->wcsdistortion!=NULL) );
if( ( checkkeys
&& ( p->date || p->comment || p->history || p->asis
|| p->delete || p->rename || p->update || p->write
|| p->verify || p->printallkeys || p->copykeys ) )
|| checkkeys>1 )
- error(EXIT_FAILURE, 0, "'--keyvalue', '--datetosec' and "
- "'--wcsdistortion' cannot currently be called with "
- "any other option");
+ error(EXIT_FAILURE, 0, "'--keyvalue', '--datetosec', "
+ "'--wcscoordsys' and '--wcsdistortion' cannot "
+ "currently be called with any other option");
+
+ /* Give an ID to recognized coordinate systems. */
+ if(p->wcscoordsys)
+ p->coordsysid=gal_wcs_coordsys_from_string(p->wcscoordsys);
/* Identify the requested distortion. Note that this also acts as a
sanity check because it will crash with an error if the given
diff --git a/bin/fits/ui.h b/bin/fits/ui.h
index 4ee5333..61f70af 100644
--- a/bin/fits/ui.h
+++ b/bin/fits/ui.h
@@ -77,6 +77,7 @@ enum option_keys_enum
UI_KEY_SKYCOVERAGE,
UI_KEY_OUTHDU,
UI_KEY_COPYKEYS,
+ UI_KEY_WCSCOORDSYS,
UI_KEY_PRIMARYIMGHDU,
UI_KEY_WCSDISTORTION,
};
diff --git a/configure.ac b/configure.ac
index 173763d..bd2e865 100644
--- a/configure.ac
+++ b/configure.ac
@@ -555,6 +555,12 @@ AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_OBSFIX],
[$has_wcslib_obsfix],
[WCSLIB comes with OBSFIX macro])
AC_SUBST(HAVE_WCSLIB_OBSFIX, [$has_wcslib_obsfix])
+# If the WCS library has the 'wcsccs' function.
+AC_CHECK_LIB([wcs], [wcsccs], [has_wcslib_wcsccs=1],
+ [has_wcslib_wcsccs=0; anywarnings=yes], [-lcfitsio -lm])
+AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_WCSCCS], [$has_wcslib_wcsccs],
+ [WCSLIB comes with wcsccs])
+AC_SUBST(HAVE_WCSLIB_WCSCCS, [$has_wcslib_wcsccs])
# If the pthreads library has 'pthread_barrier_wait'.
AC_CHECK_LIB([pthread], [pthread_barrier_wait], [has_pthread_barrier=1],
@@ -1115,6 +1121,17 @@ AS_IF([test x$enable_guide_message = xyes],
AS_ECHO([" operations you can ignore this warning."])
AS_ECHO([]) ])
+ AS_IF([test "x$has_wcslib_wcsccs" = "x0"],
+ [dependency_notice=yes
+ AS_ECHO([" - WCSLIB
(https://www.atnf.csiro.au/people/mcalabre/WCS) version"])
+ AS_ECHO([" on this system doesn't support conversion of
coordinate systems"])
+ AS_ECHO([" (through the 'wcsccs' function that was
introduced in WCSLIB 7.5, "])
+ AS_ECHO([" March 2021). For example converting from
equatorial J2000 to"])
+ AS_ECHO([" Galactic coordinates). This build won't crash but
the related"])
+ AS_ECHO([" functionalities in Gnuastro will be disabled. If
you don't need"])
+ AS_ECHO([" such operations you can ignore this warning."])
+ AS_ECHO([]) ])
+
AS_IF([test "x$has_libjpeg" = "xno"],
[dependency_notice=yes
AS_ECHO([" - libjpeg (http://ijg.org), could not be linked
with in your library"])
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 61b9503..e3ad925 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9506,7 +9506,46 @@ In this case (following the GNU C Library), this option
will make the following
This is a very useful option for operations on the FITS date values, for
example sorting FITS files by their dates, or finding the time difference
between two FITS files.
The advantage of working with the Unix epoch time is that you don't have to
worry about calendar details (for example the number of days in different
months, or leap years, etc).
-@item --wcsdistortion STR
+@item --wcscoordsys=STR
+@cindex Galactic coordinate system
+@cindex Ecliptic coordinate system
+@cindex Equatorial coordinate system
+@cindex Supergalactic coordinate system
+@cindex Coordinate system: Galactic
+@cindex Coordinate system: Ecliptic
+@cindex Coordinate system: Equatorial
+@cindex Coordinate system: Supergalactic
+Convert the coordinate system of the image's world coordinate system (WCS) to
the given coordinate system (@code{STR}) and write it into the file given to
@option{--output} (or an automatically named file if no @option{--output} has
been given).
+
+For example with the command below, @file{img-eq.fits} will have an identical
dataset (pixel values) as @file{image.fits}.
+However, the WCS coordinate system of @file{img-eq.fits} will be the
equatorial coordinate system in the Julian calendar epoch 2000 (which is the
most common epoch used today).
+Fits will automatically extract the current coordinate system of
@file{image.fits} and as long as its one of the recognized coordinate systems
listed below, it will do the conversion.
+
+@example
+$ astfits image.fits --coordsys=eq-j2000 --output=img-eq.fits
+@end example
+
+The currently recognized coordinate systems are listed below (the most common
one today is @code{eq-j2000}):
+
+@table @code
+@item eq-j2000
+2000.0 (Julian-year) equatorial coordinates.
+@item eq-b1950
+1950.0 (Besselian-year) equatorial coordinates.
+@item ec-j2000
+2000.0 (Julian-year) ecliptic coordinates.
+@item ec-b1950
+1950.0 (Besselian-year) ecliptic coordinates.
+@item galactic
+Galactic coordinates.
+@item supergalactic
+Supergalactic coordinates.
+@end table
+
+The Equatorial and Ecliptic coordinate systems are defined by the mean equator
and equinox epoch: either the Besselian year 1950.0, or the Julian year 2000.
+For more on their difference and links for further reading about epochs in
astronomy, please see the description in
@url{https://en.wikipedia.org/wiki/Epoch_(astronomy), Wikipedia}.
+
+@item --wcsdistortion=STR
@cindex WCS distortion
@cindex Distortion, WCS
@cindex SIP WCS distortion
@@ -25606,6 +25645,26 @@ 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_COORDSYS_EQB1950
+@deffnx Macro GAL_WCS_COORDSYS_EQJ2000
+@deffnx Macro GAL_WCS_COORDSYS_ECB1950
+@deffnx Macro GAL_WCS_COORDSYS_ECJ2000
+@deffnx Macro GAL_WCS_COORDSYS_GALACTIC
+@deffnx Macro GAL_WCS_COORDSYS_SUPERGALACTIC
+@deffnx Macro GAL_WCS_COORDSYS_INVALID
+@cindex Galactic coordinate system
+@cindex Ecliptic coordinate system
+@cindex Equatorial coordinate system
+@cindex Supergalactic coordinate system
+@cindex Coordinate system: Galactic
+@cindex Coordinate system: Ecliptic
+@cindex Coordinate system: Equatorial
+@cindex Coordinate system: Supergalactic
+Recognized WCS coordinate systems in Gnuastro.
+@code{EQ} and @code{EC} stand for the Equatorial and Ecliptic coordinate
systems.
+In the equatorial and ecliptic coordinates, @code{B1950} stands for the
Besselian 1950 epoch and @code{J2000} stands for the Julian 2000 epoch.
+@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.
@@ -25735,6 +25794,21 @@ 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_coordsys_from_string (char @code{*coordsys})
+Convert the given string to Gnuastro's integer-based WCS coordinate system
identifier (one of the @code{GAL_WCS_COORDSYS_*}, listed above).
+The expected strings can be seen in the description of the
@option{--wcscoordsys} option of the Fits program, see @ref{Keyword inspection
and manipulation}.
+@end deftypefun
+
+@deftypefun int gal_wcs_coordsys_identify (struct wcsprm @code{*wcs})
+Read the given WCS structure and return its coordinate system as one of
Gnuastro's WCS coordinate system identifiers (the macros
@code{GAL_WCS_COORDSYS_*}, listed above).
+@end deftypefun
+
+@deftypefun {struct wcsprm *} gal_wcs_coordsys_convert (struct wcsprm
@code{*inwcs}, int @code{coordsysid})
+Return a newly allocated WCS structure with the @code{coordsysid} coordinate
system identifier.
+The Gnuastro WCS distortion identifiers are defined in the
@code{GAL_WCS_COORDSYS_*} macros mentioned above.
+Since the returned dataset is newly allocated, if you don't need the original
dataset after this, use the WCSLIB library function @code{wcsfree} to free the
input, for example @code{wcsfree(inwcs)}.
+@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/Makefile.am b/lib/Makefile.am
index 5b427c6..4b4b297 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -150,6 +150,7 @@ gnuastro/config.h: Makefile $(internaldir)/config.h.in
-e 's|@HAVE_WCSLIB_DIS_H[@]|$(HAVE_WCSLIB_DIS_H)|g' \
-e 's|@HAVE_WCSLIB_MJDREF[@]|$(HAVE_WCSLIB_MJDREF)|g' \
-e 's|@HAVE_WCSLIB_OBSFIX[@]|$(HAVE_WCSLIB_OBSFIX)|g' \
+ -e 's|@HAVE_WCSLIB_WCSCCS[@]|$(HAVE_WCSLIB_WCSCCS)|g' \
-e 's|@HAVE_WCSLIB_VERSION[@]|$(HAVE_WCSLIB_VERSION)|g' \
-e 's|@HAVE_PTHREAD_BARRIER[@]|$(HAVE_PTHREAD_BARRIER)|g' \
-e 's|@RESTRICT_REPLACEMENT[@]|$(RESTRICT_REPLACEMENT)|g' \
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index beef0bc..1ae01b9 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -70,6 +70,19 @@ enum gal_wcs_distortions
GAL_WCS_DISTORTION_WAT, /* The WAT polynomial distortion. */
};
+/* Macros to identify coordinate system for convesions. */
+enum gal_wcs_coordsys
+{
+ GAL_WCS_COORDSYS_INVALID, /* Invalid (=0 by C standard). */
+
+ GAL_WCS_COORDSYS_EQB1950, /* Equatorial B1950 */
+ GAL_WCS_COORDSYS_EQJ2000, /* Equatorial J2000 */
+ GAL_WCS_COORDSYS_ECB1950, /* Ecliptic B1950 */
+ GAL_WCS_COORDSYS_ECJ2000, /* Ecliptic J2000 */
+ GAL_WCS_COORDSYS_GALACTIC, /* Galactic */
+ GAL_WCS_COORDSYS_SUPERGALACTIC, /* Super-galactic */
+};
+
/* Macros to identify the type of distortion for conversions. */
enum gal_wcs_linear_matrix
{
@@ -121,6 +134,20 @@ gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm
*wcs);
*********** Distortions ***********
*************************************************************/
int
+gal_wcs_coordsys_from_string(char *coordsys);
+
+int
+gal_wcs_coordsys_identify(struct wcsprm *inwcs);
+
+struct wcsprm *
+gal_wcs_coordsys_convert(struct wcsprm *inwcs, int coordsysid);
+
+
+
+/*************************************************************
+ *********** Distortions ***********
+ *************************************************************/
+int
gal_wcs_distortion_from_string(char *distortion);
char *
diff --git a/lib/wcs.c b/lib/wcs.c
index 40436de..77d711b 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -629,6 +629,406 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
+
+/*************************************************************
+ *********** Coordinate system ***********
+ *************************************************************/
+int
+gal_wcs_coordsys_from_string(char *coordsys)
+{
+ if( !strcmp(coordsys,"eq-j2000") ) return GAL_WCS_COORDSYS_EQJ2000;
+ else if( !strcmp(coordsys,"eq-b1950") ) return GAL_WCS_COORDSYS_EQB1950;
+ else if( !strcmp(coordsys,"ec-j2000") ) return GAL_WCS_COORDSYS_ECJ2000;
+ else if( !strcmp(coordsys,"ec-b1950") ) return GAL_WCS_COORDSYS_ECB1950;
+ else if( !strcmp(coordsys,"galactic") ) return GAL_WCS_COORDSYS_GALACTIC;
+ else if( !strcmp(coordsys,"supergalactic") )
+ return GAL_WCS_COORDSYS_SUPERGALACTIC;
+ else
+ error(EXIT_FAILURE, 0, "WCS coordinate system name '%s' not "
+ "recognized, currently recognized names are 'eq-j2000', "
+ "'eq-b1950', 'galactic' and 'supergalactic'", coordsys);
+
+ /* Control should not reach here. */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+ "problem. Control should not reach the end of this function",
+ __func__, PACKAGE_BUGREPORT);
+ return GAL_WCS_COORDSYS_INVALID;
+}
+
+
+
+
+/* Identify the coordinate system of the WCS. */
+int
+gal_wcs_coordsys_identify(struct wcsprm *wcs)
+{
+ /* Equatorial (we are keeping the dash ('-') to make sure it is a
+ standard). */
+ if ( !strncmp(wcs->ctype[0], "RA---", 5)
+ && !strncmp(wcs->ctype[1], "DEC--", 5) )
+ {
+ if ( !strncmp(wcs->radesys, "FK4", 3) )
+ return GAL_WCS_COORDSYS_EQB1950;
+ else if ( !strncmp(wcs->radesys, "FK5", 3) )
+ return GAL_WCS_COORDSYS_EQJ2000;
+ else
+ error(EXIT_FAILURE, 0, "%s: the '%s' value for 'RADESYS' is "
+ "not yet implemented! Please contact us at %s to "
+ "implement it", __func__, wcs->radesys, PACKAGE_BUGREPORT);
+ }
+
+ /* Ecliptic. */
+ else if ( !strncmp(wcs->ctype[0], "ELON-", 5)
+ && !strncmp(wcs->ctype[1], "ELAT-", 5) )
+ if ( !strncmp(wcs->radesys, "FK4", 3) )
+ return GAL_WCS_COORDSYS_ECB1950;
+ else if ( !strncmp(wcs->radesys, "FK5", 3) )
+ return GAL_WCS_COORDSYS_ECJ2000;
+ else
+ error(EXIT_FAILURE, 0, "%s: the '%s' value for 'RADESYS' is "
+ "not yet implemented! Please contact us at %s to "
+ "implement it", __func__, wcs->radesys, PACKAGE_BUGREPORT);
+
+ /* Galactic. */
+ else if ( !strncmp(wcs->ctype[0], "GLON-", 5)
+ && !strncmp(wcs->ctype[1], "GLAT-", 5) )
+ return GAL_WCS_COORDSYS_GALACTIC;
+
+ /* SuperGalactic. */
+ else if ( !strncmp(wcs->ctype[0], "SLON-", 5)
+ && !strncmp(wcs->ctype[1], "SLAT-", 5) )
+ return GAL_WCS_COORDSYS_SUPERGALACTIC;
+
+ /* Other. */
+ else
+ error(EXIT_FAILURE, 0, "%s: the CTYPE values '%s' and '%s' are "
+ "not yet implemented! Please contact us at %s to "
+ "implement it", __func__, wcs->ctype[0], wcs->ctype[1],
+ PACKAGE_BUGREPORT);
+
+ /* Control should not reach here. */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+ "problem. Control should not reach the end of this function",
+ __func__, PACKAGE_BUGREPORT);
+ return GAL_WCS_COORDSYS_INVALID;
+}
+
+
+
+
+
+/* Set the pole coordinates (current values taken from the WCSLIB
+ manual.
+ lng2p1: pole of input (1) system in output (2) system's logitude.
+ lat2p1: pole of input (1) system in output (2) system's latitude.
+ lng1p2: pole of output (2) system in input (1) system's longitude.
+
+ Values from NED (inspired by WCSLIB manual's example).
+ https://ned.ipac.caltech.edu/coordinate_calculator
+
+ longi (deg) latit (deg) OUTPUT INPUT
+ ----- ----- ------ -----
+ (------------, -----------) B1950 equ. coords. of B1950 equ. pole.
+ (180.31684301, 89.72174782) J2000 equ. coords. of B1950 equ. pole.
+ (90.000000000, 66.55421111) B1950 ecl. coords. of B1950 equ. pole.
+ (90.699521110, 66.56068919) J2000 ecl. coords. of B1950 equ. pole.
+ (123.00000000, 27.40000000) Galactic coords. of B1950 equ. pole.
+ (26.731537070, 15.64407736) Supgalactic coords. of B1950 equ. pole.
+
+ (359.68621044, 89.72178502) B1950 equ. coords. of J2000 equ. pole.
+ (------------, -----------) J2000 equ. coords. of J2000 equ. pole.
+ (89.300755510, 66.55417728) B1950 ecl. coords. of J2000 equ. pole.
+ (90.000000000, 66.56070889) J2000 ecl. coords. of J2000 equ. pole.
+ (122.93200023, 27.12843056) Galactic coords. of J2000 equ. pole.
+ (26.450516650, 15.70886131) Supgalactic coords. of J2000 equ. pole.
+
+ (270.00000000, 66.55421111) B1950 equ. coords. of B1950 ecl. pole.
+ (269.99920697, 66.55421892) J2000 equ. coords. of B1950 ecl. pole.
+ (------------, -----------) B1950 ecl. coords. of B1950 ecl. pole.
+ (267.21656404, 89.99350237) J2000 ecl. coords. of B1950 ecl. pole.
+ (96.376479150, 29.81195400) Galactic coords. of B1950 ecl. pole.
+ (33.378919140, 38.34766498) Supgalactic coords. of B1950 ecl. pole.
+
+ (270.00099211, 66.56069675) B1950 equ. coords. of J2000 ecl. pole.
+ (270.00000000, 66.56070889) J2000 equ. coords. of J2000 ecl. pole.
+ (86.517962160, 89.99350236) B1950 ecl. coords. of J2000 ecl. pole.
+ (------------, -----------) J2000 ecl. coords. of J2000 ecl. pole.
+ (96.383958840, 29.81163604) Galactic coords. of J2000 ecl. pole.
+ (33.376119480, 38.34154959) Supgalactic coords. of J2000 ecl. pole.
+
+ (192.25000000, 27.40000000) B1950 equ. coords. of Galactic pole.
+ (192.85949646, 27.12835323) J2000 equ. coords. of Galactic pole.
+ (179.32094769, 29.81195400) B1950 ecl. coords. of Galactic pole.
+ (180.02317894, 29.81153742) J2000 ecl. coords. of Galactic pole.
+ (------------, -----------) Galactic coords. of Galactic pole.
+ (90.000000000, 6.320000000) Supgalactic coords. of Galactic pole.
+
+ (283.18940711, 15.64407736) B1950 equ. coords. of SupGalactic pole.
+ (283.75420420, 15.70894043) J2000 equ. coords. of SupGalactic pole.
+ (286.26975051, 38.34766498) B1950 ecl. coords. of SupGalactic pole.
+ (286.96654469, 38.34158720) J2000 ecl. coords. of SupGalactic pole.
+ (47.370000000, 6.320000000) Galactic coords. of SupGalactic pole.
+ (------------, -----------) Supgalactic coords. of SupGalactic pole.
+ */
+static void
+wcs_coordsys_insys_pole_in_outsys(int insys, int outsys, double *lng2p1,
+ double *lat2p1, double *lng1p2)
+{
+ switch( insys )
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=180.31684301; *lat2p1=89.72174782; *lng1p2=359.68621044;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=90.000000000; *lat2p1=66.55421111; *lng1p2=270.00000000;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=90.699521110; *lat2p1=66.56068919; *lng1p2=270.00099211;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=123.00000000; *lat2p1=27.40000000; *lng1p2=192.25000000;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=26.731537070; *lat2p1=15.64407736; *lng1p2=283.18940711;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input EQB1950)", __func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=359.68621044; *lat2p1=89.72178502; *lng1p2=180.31684301;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=89.300755510; *lat2p1=66.55417728; *lng1p2=269.99920697;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=90.000000000; *lat2p1=66.56070889; *lng1p2=270.00000000;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=122.93200023; *lat2p1=27.12843056; *lng1p2=192.85949646;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=26.450516650; *lat2p1=15.70886131; *lng1p2=283.75420420;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input EQJ2000)", __func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ case GAL_WCS_COORDSYS_ECB1950:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=270.00000000; *lat2p1=66.55421111; *lng1p2=90.000000000;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=269.99920697; *lat2p1=66.55421892; *lng1p2=89.300755510;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=267.21656404; *lat2p1=89.99350237; *lng1p2=86.517962160;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=96.383958840; *lat2p1=29.81163604; *lng1p2=179.32094769;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=33.378919140; *lat2p1=38.34766498; *lng1p2=286.26975051;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input ECB1950)", __func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=270.00099211; *lat2p1=66.56069675; *lng1p2=90.699521110;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=270.00000000; *lat2p1=66.56070889; *lng1p2=90.000000000;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=86.517962160; *lat2p1=89.99350236; *lng1p2=267.21656404;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=96.383958840; *lat2p1=29.81163604; *lng1p2=180.02317894;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=33.376119480; *lat2p1=38.34154959; *lng1p2=286.96654469;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input ECJ2000)", __func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=192.25000000; *lat2p1=27.40000000; *lng1p2=123.00000000;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=192.85949646; *lat2p1=27.12835323; *lng1p2=122.93200023;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=179.32094769; *lat2p1=29.81195400; *lng1p2=96.376479150;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=180.02317894; *lat2p1=29.81153742; *lng1p2=96.383958840;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=90.000000000; *lat2p1=6.320000000; *lng1p2=47.370000000;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input GALACTIC)", __func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ switch( outsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *lng2p1=283.18940711; *lat2p1=15.64407736; *lng1p2=26.731537070;
return;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *lng2p1=283.75420420; *lat2p1=15.70894043; *lng1p2=26.450516650;
return;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *lng2p1=286.26975051; *lat2p1=38.34766498; *lng1p2=33.378919140;
return;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *lng2p1=286.96654469; *lat2p1=38.34158720; *lng1p2=33.376119480;
return;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *lng2p1=47.370000000; *lat2p1=6.320000000; *lng1p2=90.000000000;
return;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN;
return;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'outsys' (input SUPERGALACTIC)",
__func__,
+ PACKAGE_BUGREPORT, outsys);
+ }
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'insys'", __func__,
+ PACKAGE_BUGREPORT, insys);
+ }
+
+}
+
+
+
+
+
+static void
+wcs_coordsys_ctypes(int coordsys, char **clng, char **clat, char **radesys)
+{
+ switch( coordsys)
+ {
+ case GAL_WCS_COORDSYS_EQB1950:
+ *clng="RA"; *clat="DEC"; *radesys="FK4"; break;
+ case GAL_WCS_COORDSYS_EQJ2000:
+ *clng="RA"; *clat="DEC"; *radesys="FK5"; break;
+ case GAL_WCS_COORDSYS_ECB1950:
+ *clng="ELON"; *clat="ELAT"; *radesys="FK4"; break;
+ case GAL_WCS_COORDSYS_ECJ2000:
+ *clng="ELON"; *clat="ELAT"; *radesys="FK5"; break;
+ case GAL_WCS_COORDSYS_GALACTIC:
+ *clng="GLON"; *clat="GLAT"; *radesys=NULL; break;
+ case GAL_WCS_COORDSYS_SUPERGALACTIC:
+ *clng="SLON"; *clat="SLAT"; *radesys=NULL; break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The code '%d' isn't a recognized WCS "
+ "coordinate system ID for 'coordsys'", __func__,
+ PACKAGE_BUGREPORT, coordsys);
+ }
+}
+
+
+
+
+/* Convert the coordinate system. */
+struct wcsprm *
+gal_wcs_coordsys_convert(struct wcsprm *wcs, int outcoordsys)
+{
+ int incoordsys;
+ char *alt=NULL; /* Only concerned with primary wcs. */
+ double equinox=0.0f; /* To preserve current value. */
+ struct wcsprm *out=NULL;
+ char *clng, *clat, *radesys;
+ double lng2p1=NAN, lat2p1=NAN, lng1p2=NAN;
+
+
+ /* Just incase the input is a NULL pointer. */
+ if(wcs==NULL) return NULL;
+
+ /* Get the input's coordinate system and see if it should be converted at
+ all or not (if the output coordinate system is different). If its the
+ same, just copy the input and return. */
+ incoordsys=gal_wcs_coordsys_identify(wcs);
+ if(incoordsys==outcoordsys)
+ {
+ out=gal_wcs_copy(wcs);
+ return out;
+ }
+
+ /* Find the necessary pole coordinates. Note that we have already
+ accounted for the fact that the input and output coordinate systems
+ may be the same above, so the NaN outputs will never occur here. */
+ wcs_coordsys_insys_pole_in_outsys(incoordsys, outcoordsys,
+ &lng2p1, &lat2p1, &lng1p2);
+
+ /* Find the necessary CTYPE names of the output. */
+ wcs_coordsys_ctypes(outcoordsys, &clng, &clat, &radesys);
+
+ /* Convert the WCS's coordinate system (if 'wcsccs' is available). */
+#if GAL_CONFIG_HAVE_WCSLIB_WCSCCS
+ out=gal_wcs_copy(wcs);
+ wcsccs(out, lng2p1, lat2p1, lng1p2, clng, clat, radesys, equinox, alt);
+#else
+
+ /* Just to avoid compiler warnings for 'equinox' and 'alt'. */
+ if(alt) lng2p1+=equinox;
+
+ /* Print error message and abort. */
+ error(EXIT_FAILURE, 0, "%s: the 'wcsccs' function isn't available "
+ "in the version of WCSLIB that this Gnuastro was built with "
+ "('wcsccs' was first available in WCSLIB 7.5, released on "
+ "March 2021). Therefore, Gnuastro can't preform the WCS "
+ "coordiante system conversion in the WCS. Please update your "
+ "WCSLIB and re-build Gnuastro with it to use this feature. "
+ "You can follow the instructions here to install the latest "
+ "version of WCSLIB:\n"
+ " https://www.gnu.org/software/gnuastro/manual/html_node/";
+ "WCSLIB.html\n"
+ "And then re-build Gnuastro as described here:\n"
+ " https://www.gnu.org/software/gnuastro/manual/";
+ "html_node/Quick-start.html\n\n",
+ __func__);
+#endif
+
+ /* Return. */
+ return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/*************************************************************
*********** Distortions ***********
*************************************************************/
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 20f1d41: Fits: new --wcscoordsys option to convert coordinate systems,
Mohammad Akhlaghi <=