Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix make_fitswcs_transform error when handling header with no PC #142

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Fix make_fitswcs_transform error when handling header with no PC #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
713abc8
Fixes from flake8 check.
pllim Mar 20, 2018
75e97e9
make_fitswcs_transform to transform dict hdr properly
pllim Mar 20, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions gwcs/tests/test_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,30 @@ def test_unknown_ctype():
assert_allclose(b, expected[1], atol=10**-8)


# https://github.com/spacetelescope/gwcs/issues/139
def test_transform_hdr_dict():
"""2D test case from Ginga."""
header = {
'BITPIX': -32, 'BLANK': -32768, 'BUNIT': 'ADU', 'EXTEND': False,
'CD1_1': -5.611e-05, 'CD1_2': 0.0, 'CD2_1': 0.0, 'CD2_2': 5.611e-05,
'CDELT1': -5.611e-05, 'CDELT2': 5.611e-05,
'CRPIX1': 5276.0, 'CRPIX2': 25.0,
'CRVAL1': 299.91736667, 'CRVAL2': 22.68769444,
'CTYPE1': 'RA---TAN', 'CTYPE2': 'DEC--TAN',
'CUNIT1': 'degree', 'CUNIT2': 'degree', 'EQUINOX': 2000.0,
'NAXIS': 2, 'NAXIS1': 2272, 'NAXIS2': 4273,
'RA': '19:59:40.168', 'RA2000': '19:59:40.168',
'DEC': '+22:41:15.70', 'DEC2000': '+22:41:15.70', 'RADECSYS': 'FK5',
'SIMPLE': True, 'WCS-ORIG': 'SUBARU Toolkit'}
w = gwutils.make_fitswcs_transform(header, dict_to_fits=True)
xy_ans = np.array([120, 100])
radec_deg_ans = (300.2308791294835, 22.691653517073615)

# TODO: Check with Nadia
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nden , I don't understand

  • why w() result can't agree for default rtol=1e-7
  • why w.inverse() would give me original XY + 1 and not XY

Am I using it wrong?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

GWCS by default uses 0-based pixel coordinates. CRPIX values are 1-based. So if you want to compare the inverse transform, it's going to return 0-based values.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Though maybe I'm misunderstanding the problem?

Copy link
Contributor Author

@pllim pllim Mar 21, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am just generally confused... What I am seeing is this:

ra, deg = w(x, y)

but

w.inverse(ra, dec) = x + 1, y + 1

All I want to know is if this is the expected behavior.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jdavies-st is correct. When you are mixing FITS and gwcs you have to take care of 0- vs 1- based coordinates. The gwcs transform is 0-based. However, the values of crpix are 1-based. If you change the CRPIX values to be 0-based, i.e. subtract 1 from them , you will get agreement (you'll get better agreement in the forward direction too).

BTW, the header has incorrect FITS WCS to start with. It is not allowed to mix CD matrix and CDELT values and in this case the CDELT values are ignored and assumed to be 1.

I don't want to encourage changes and "fixes" to these functions as I think it's time to implement full conversion between the two objects.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, so I have to -1 from the keyword value, not inputs. Gotcha.

So, should I close this and wait for your "full conversion" then?

assert_allclose(w(*xy_ans), radec_deg_ans, rtol=1e-5)
assert_allclose(w.inverse(*radec_deg_ans), xy_ans + 1)


def test_get_axes():
wcsinfo = {'CTYPE': np.array(['MRSAL1A', 'MRSBE1A', 'WAVE'])}
cel, spec, other = gwutils.get_axes(wcsinfo)
Expand Down
71 changes: 46 additions & 25 deletions gwcs/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,11 @@
import functools
import numpy as np
from astropy.modeling import models as astmodels
from astropy.modeling.models import Mapping
from astropy.modeling import core, projections
from astropy.modeling import core, projections, is_separable
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing import is a real bug found by my flake8 checker. If you want, I can put this in as separate PR. Though none of your existing tests caught it, so probably not important?

from astropy.io import fits
from astropy import coordinates as coords
from astropy import units as u

from astropy.utils.decorators import deprecated


# these ctype values do not include yzLN and yzLT pairs
sky_pairs = {"equatorial": ["RA", "DEC"],
"ecliptic": ["ELON", "ELAT"],
Expand Down Expand Up @@ -89,7 +85,8 @@ def _get_values(units, *args):
----------
units : str or `~astropy.units.Unit`
Units of the wcs object.
The input values are converted to ``units`` before the values are returned.
The input values are converted to ``units`` before
the values are returned.
"""
val = []
values = []
Expand All @@ -116,7 +113,8 @@ def _compute_lon_pole(skycoord, projection):
Compute the longitude of the celestial pole of a standard frame in the
native frame.

This angle then can be used as one of the Euler angles (the other two being skyccord)
This angle then can be used as one of the Euler angles
(the other two being skycoord)
to rotate the native frame into the standard frame ``skycoord.frame``.

Parameters
Expand Down Expand Up @@ -153,7 +151,8 @@ def _compute_lon_pole(skycoord, projection):
else:
lon_pole = 180
else:
raise UnsupportedProjectionError("Projection {0} is not supported.".format(projection))
raise UnsupportedProjectionError(
"Projection {0} is not supported.".format(projection))
if unit is not None:
lon_pole = lon_pole * unit
return lon_pole
Expand All @@ -166,7 +165,8 @@ def get_projcode(wcs_info):
return None
projcode = wcs_info['CTYPE'][sky_axes[0]][5:8].upper()
if projcode not in projections.projcodes:
raise UnsupportedProjectionError('Projection code %s, not recognized' % projcode)
raise UnsupportedProjectionError(
'Projection code %s, not recognized' % projcode)
#projcode = None
return projcode

Expand Down Expand Up @@ -291,18 +291,23 @@ def get_axes(header):


def _is_skysys_consistent(ctype, sky_inmap):
""" Determine if the sky axes in CTYPE mathch to form a standard celestial system."""
"""
Determine if the sky axes in CTYPE mathch to form a
standard celestial system.
"""

for item in sky_pairs.values():
if ctype[sky_inmap[0]] == item[0]:
if ctype[sky_inmap[1]] != item[1]:
raise ValueError(
"Inconsistent ctype for sky coordinates {0} and {1}".format(*ctype))
"Inconsistent ctype for sky coordinates "
"{0} and {1}".format(*ctype))
break
elif ctype[sky_inmap[1]] == item[0]:
if ctype[sky_inmap[0]] != item[1]:
raise ValueError(
"Inconsistent ctype for sky coordinates {0} and {1}".format(*ctype))
"Inconsistent ctype for sky coordinates "
"{0} and {1}".format(*ctype))
sky_inmap = sky_inmap[::-1]
break

Expand All @@ -319,7 +324,7 @@ def _is_skysys_consistent(ctype, sky_inmap):
}


def make_fitswcs_transform(header):
def make_fitswcs_transform(header, dict_to_fits=False):
"""
Create a basic FITS WCS transform.
It does not include distortions.
Expand All @@ -329,13 +334,23 @@ def make_fitswcs_transform(header):
header : astropy.io.fits.Header or dict
FITS Header (or dict) with basic WCS information

dict_to_fits : bool
Convert given header (if dict) to FITS header
before further processing.

"""
if not isinstance(header, (fits.Header, dict)):
raise TypeError("Expected a FITS Header or a dict.")

if isinstance(header, dict) and dict_to_fits:
header = fits.Header(header)

if isinstance(header, fits.Header):
# This converts CD matrix into PC for linear transformation below.
wcs_info = read_wcs_from_header(header)
elif isinstance(header, dict):
wcs_info = header
else:
raise TypeError("Expected a FITS Header or a dict.")
wcs_info = header

transforms = []
wcs_linear = fitswcs_linear(wcs_info)
transforms.append(wcs_linear)
Expand Down Expand Up @@ -404,7 +419,7 @@ def fitswcs_linear(header):

if not wcs_info['has_cd']:
# Do not compute scaling since CDELT* = 1 if CD is present.
scaling_models = [astmodels.Scale(scale, name='cdelt' + str(i + 1)) \
scaling_models = [astmodels.Scale(scale, name='cdelt' + str(i + 1))
for i, scale in enumerate(cdelt)]

scaling = functools.reduce(lambda x, y: x & y, scaling_models)
Expand Down Expand Up @@ -443,7 +458,8 @@ def fitswcs_nonlinear(header):
# TODO: write "def compute_lonpole(projcode, l)"
# Set a defaul tvalue for now
thetap = 180
n2c = astmodels.RotateNative2Celestial(phip, lonp, thetap, name="crval")
n2c = astmodels.RotateNative2Celestial(
phip, lonp, thetap, name="crval")
transforms.append(n2c)
if transforms:
return functools.reduce(core._model_oper('|'), transforms)
Expand Down Expand Up @@ -497,21 +513,23 @@ def separable_axes(wcsobj, start_frame=None, end_frame=None):
Computes the separability of axes in ``end_frame``.

Returns a 1D boolean array of size frame.naxes where True means
the axis is completely separable and False means the axis is nonseparable
from at least one other axis.
the axis is completely separable and False means the axis is
nonseparable from at least one other axis.

Parameters
----------
wcsobj : `~gwcs.wcs.WCS`
WCS object
start_frame : `~gwcs.coordinate_frames.CoordinateFrame`
A frame in the WCS pipeline.
The transform between start_frame and the end frame is used to compute the
The transform between start_frame and the end frame
is used to compute the
mapping inputs: outputs.
If None the input_frame is used as start_frame.
end_frame : `~gwcs.coordinate_frames.CoordinateFrame`
A frame in the WCS pipeline.
The transform between start_frame and the end frame is used to compute the
The transform between start_frame and the end frame
is used to compute the
mapping inputs: outputs.
If None wcsobj.output_frame is used.

Expand All @@ -525,15 +543,18 @@ def separable_axes(wcsobj, start_frame=None, end_frame=None):
start_frame = wcsobj.input_frame
else:
if start_frame not in wcsobj.available_frames:
raise ValueError("Unrecognized frame {0}".format(start_frame))
raise ValueError(
"Unrecognized frame {0}".format(start_frame))
if end_frame is None:
end_frame = wcsobj.output_frame
else:
if end_frame not in wcsobj.available_frames:
raise ValueError("Unrecognized frame {0}".format(end_frame))
raise ValueError(
"Unrecognized frame {0}".format(end_frame))
transform = wcsobj.get_transform(start_frame, end_frame)
else:
raise ValueError("A starting frame is needed to determine separability of axes.")
raise ValueError("A starting frame is needed to determine "
"separability of axes.")

sep = is_separable(transform)
return [sep[ax] for ax in end_frame.axes_order]