Python astropy.coordinates.ICRS Examples

def test_skycoord_to_pixel_distortions(mode):

    # Import astropy.coordinates here to avoid circular imports
    from astropy.coordinates import SkyCoord

    header = get_pkg_data_filename('data/sip.fits')
    with pytest.warns(FITSFixedWarning):
        wcs = WCS(header)

    ref = SkyCoord(202.50 * u.deg, 47.19 * u.deg, frame='icrs')

    xp, yp = skycoord_to_pixel(ref, wcs, mode=mode)

    # WCS is in FK5 so we need to transform back to ICRS
    new = pixel_to_skycoord(xp, yp, wcs, mode=mode).transform_to('icrs')

    assert_allclose(new.ra.degree, ref.ra.degree)
    assert_allclose(new.dec.degree, ref.dec.degree) 

def test_matching_function_3d_and_sky():
    from astropy.coordinates import ICRS
    from astropy.coordinates.matching import match_coordinates_3d, match_coordinates_sky

    cmatch = ICRS([4, 2.1]*u.degree, [0, 0]*u.degree, distance=[1, 5] * u.kpc)
    ccatalog = ICRS([1, 2, 3, 4]*u.degree, [0, 0, 0, 0]*u.degree, distance=[1, 1, 1, 5] * u.kpc)

    idx, d2d, d3d = match_coordinates_3d(cmatch, ccatalog)
    npt.assert_array_equal(idx, [2, 3])

    assert_allclose(d2d, [1, 1.9] * u.deg)
    assert np.abs(d3d[0].to_value(u.kpc) - np.radians(1)) < 1e-6
    assert np.abs(d3d[1].to_value(u.kpc) - 5*np.radians(1.9)) < 1e-5

    idx, d2d, d3d = match_coordinates_sky(cmatch, ccatalog)
    npt.assert_array_equal(idx, [3, 1])

    assert_allclose(d2d, [0, 0.1] * u.deg)
    assert_allclose(d3d, [4, 4.0000019] * u.kpc) 

def test_matching_function():
    from astropy.coordinates import ICRS
    from astropy.coordinates.matching import match_coordinates_3d
    # this only uses match_coordinates_3d because that's the actual implementation

    cmatch = ICRS([4, 2.1]*u.degree, [0, 0]*u.degree)
    ccatalog = ICRS([1, 2, 3, 4]*u.degree, [0, 0, 0, 0]*u.degree)

    idx, d2d, d3d = match_coordinates_3d(cmatch, ccatalog)
    npt.assert_array_equal(idx, [3, 1])
    npt.assert_array_almost_equal(d2d.degree, [0, 0.1])
    assert d3d.value[0] == 0

    idx, d2d, d3d = match_coordinates_3d(cmatch, ccatalog, nthneighbor=2)
    assert np.all(idx == 2)
    npt.assert_array_almost_equal(d2d.degree, [1, 0.9])
    npt.assert_array_less(d3d.value, 0.02) 

def assert_equal(cls, old, new):
        assert isinstance(old, ICRS)
        assert isinstance(new, ICRS)
        assert u.allclose(new.ra, old.ra)
        assert u.allclose(new.dec, old.dec) 

def from_tree(cls, node, ctx):
        angle = Angle(QuantityType.from_tree(node['ra']['wrap_angle'], ctx))
        wrap_angle = Angle(angle)
        ra = Longitude(
            node['ra']['value'],
            unit=node['ra']['unit'],
            wrap_angle=wrap_angle)
        dec = Latitude(node['dec']['value'], unit=node['dec']['unit'])

        return ICRS(ra=ra, dec=dec) 

def test_skycoord_ra_dec(tmpdir):

    ra = Longitude([1, 2, 3], unit=u.deg)  # Could also use Angle
    dec = np.array([4.5, 5.2, 6.3]) * u.deg  # Astropy Quantity
    c = SkyCoord(ra, dec, frame='icrs')
    tree = dict(coord=c)
    assert_roundtrip_tree(tree, tmpdir)

    c = SkyCoord(frame=ICRS, ra=ra, dec=dec, obstime='2001-01-02T12:34:56')
    tree = dict(coord=c)
    assert_roundtrip_tree(tree, tmpdir) 

def test_scalar_skycoord(tmpdir):

    c = SkyCoord(10, 20, unit="deg")  # defaults to ICRS frame
    tree = dict(coord=c)
    assert_roundtrip_tree(tree, tmpdir) 

def test_icrs_compound(tmpdir):

    icrs = ICRS(ra=[0, 1, 2]*units.deg, dec=[3, 4, 5]*units.deg)

    tree = {'coord': icrs}

    assert_roundtrip_tree(tree, tmpdir) 

def test_icrs_nodata(tmpdir):
    tree = {'coord': ICRS()}

    assert_roundtrip_tree(tree, tmpdir) 

def test_icrs_basic(tmpdir):
    wrap_angle = Angle(1.5, unit=units.rad)
    ra = Longitude(25, unit=units.deg, wrap_angle=wrap_angle)
    dec = Latitude(45, unit=units.deg)

    tree = {'coord': ICRS(ra=ra, dec=dec)}

    assert_roundtrip_tree(tree, tmpdir) 

def wcsapi_to_celestial_frame(wcs):
    for cls, _, kwargs, *_ in wcs.world_axis_object_classes.values():
        if issubclass(cls, SkyCoord):
            return kwargs.get('frame', ICRS())
        elif issubclass(cls, BaseCoordinateFrame):
            return cls(**kwargs) 

def test_concatenate():
    from astropy.coordinates import FK5, SkyCoord, ICRS
    from astropy.coordinates.funcs import concatenate

    # Just positions
    fk5 = FK5(1*u.deg, 2*u.deg)
    sc = SkyCoord(3*u.deg, 4*u.deg, frame='fk5')

    res = concatenate([fk5, sc])
    np.testing.assert_allclose(res.ra, [1, 3]*u.deg)
    np.testing.assert_allclose(res.dec, [2, 4]*u.deg)

    with pytest.raises(TypeError):
        concatenate(fk5)

    with pytest.raises(TypeError):
        concatenate(1*u.deg)

    # positions and velocities
    fr = ICRS(ra=10*u.deg, dec=11.*u.deg,
              pm_ra_cosdec=12*u.mas/u.yr,
              pm_dec=13*u.mas/u.yr)
    sc = SkyCoord(ra=20*u.deg, dec=21.*u.deg,
                  pm_ra_cosdec=22*u.mas/u.yr,
                  pm_dec=23*u.mas/u.yr)

    res = concatenate([fr, sc])

    with pytest.raises(ValueError):
        concatenate([fr, fk5])

    fr2 = ICRS(ra=10*u.deg, dec=11.*u.deg)
    with pytest.raises(ValueError):
        concatenate([fr, fr2]) 

def test_constellations(recwarn):
    from astropy.coordinates import ICRS, FK5, SkyCoord
    from astropy.coordinates.funcs import get_constellation

    inuma = ICRS(9*u.hour, 65*u.deg)

    n_prewarn = len(recwarn)
    res = get_constellation(inuma)
    res_short = get_constellation(inuma, short_name=True)
    assert len(recwarn) == n_prewarn  # neither version should not make warnings

    assert res == 'Ursa Major'
    assert res_short == 'UMa'
    assert isinstance(res, str) or getattr(res, 'shape', None) == tuple()

    # these are taken from the ReadMe for Roman 1987
    ras = [9, 23.5, 5.12, 9.4555, 12.8888, 15.6687, 19, 6.2222]
    decs = [65, -20, 9.12, -19.9, 22, -12.1234, -40, -81.1234]
    shortnames = ['UMa', 'Aqr', 'Ori', 'Hya', 'Com', 'Lib', 'CrA', 'Men']

    testcoos = FK5(ras*u.hour, decs*u.deg, equinox='B1950')
    npt.assert_equal(get_constellation(testcoos, short_name=True), shortnames)

    # test on a SkyCoord, *and* test Boötes, which is special in that it has a
    # non-ASCII character
    bootest = SkyCoord(15*u.hour, 30*u.deg, frame='icrs')
    boores = get_constellation(bootest)
    assert boores == 'Boötes'
    assert isinstance(boores, str) or getattr(boores, 'shape', None) == tuple() 

def test_match_catalog_nounit():
    from .. import ICRS, CartesianRepresentation
    from ..matching import match_coordinates_sky

    i1 = ICRS([[1], [2], [3]], representation_type=CartesianRepresentation)
    i2 = ICRS([[1], [2], [4, 5]], representation_type=CartesianRepresentation)
    i, sep, sep3d = match_coordinates_sky(i1, i2)
    assert_allclose(sep3d, [1]*u.dimensionless_unscaled) 

def test_matching_method():
    from astropy.coordinates import ICRS, SkyCoord
    from astropy.utils import NumpyRNGContext
    from astropy.coordinates.matching import match_coordinates_3d, match_coordinates_sky

    with NumpyRNGContext(987654321):
        cmatch = ICRS(np.random.rand(20) * 360.*u.degree,
                      (np.random.rand(20) * 180. - 90.)*u.degree)
        ccatalog = ICRS(np.random.rand(100) * 360. * u.degree,
                       (np.random.rand(100) * 180. - 90.)*u.degree)

    idx1, d2d1, d3d1 = SkyCoord(cmatch).match_to_catalog_3d(ccatalog)
    idx2, d2d2, d3d2 = match_coordinates_3d(cmatch, ccatalog)

    npt.assert_array_equal(idx1, idx2)
    assert_allclose(d2d1, d2d2)
    assert_allclose(d3d1, d3d2)

    # should be the same as above because there's no distance, but just make sure this method works
    idx1, d2d1, d3d1 = SkyCoord(cmatch).match_to_catalog_sky(ccatalog)
    idx2, d2d2, d3d2 = match_coordinates_sky(cmatch, ccatalog)

    npt.assert_array_equal(idx1, idx2)
    assert_allclose(d2d1, d2d2)
    assert_allclose(d3d1, d3d2)

    assert len(idx1) == len(d2d1) == len(d3d1) == 20 

def test_skycoord_to_pixel(mode):

    # Import astropy.coordinates here to avoid circular imports
    from astropy.coordinates import SkyCoord

    header = get_pkg_data_contents('data/maps/1904-66_TAN.hdr', encoding='binary')
    wcs = WCS(header)

    ref = SkyCoord(0.1 * u.deg, -89. * u.deg, frame='icrs')

    xp, yp = skycoord_to_pixel(ref, wcs, mode=mode)

    # WCS is in FK5 so we need to transform back to ICRS
    new = pixel_to_skycoord(xp, yp, wcs, mode=mode).transform_to('icrs')

    assert_allclose(new.ra.degree, ref.ra.degree)
    assert_allclose(new.dec.degree, ref.dec.degree)

    # Make sure you can specify a different class using ``cls`` keyword
    class SkyCoord2(SkyCoord):
        pass

    new2 = pixel_to_skycoord(xp, yp, wcs, mode=mode,
                             cls=SkyCoord2).transform_to('icrs')

    assert new2.__class__ is SkyCoord2
    assert_allclose(new2.ra.degree, ref.ra.degree)
    assert_allclose(new2.dec.degree, ref.dec.degree) 

def test_wcs_to_celestial_frame_correlated():

    # Regression test for a bug that caused wcs_to_celestial_frame to fail when
    # the celestial axes were correlated with other axes.

    # Import astropy.coordinates here to avoid circular imports
    from astropy.coordinates.builtin_frames import ICRS

    mywcs = WCS(naxis=3)
    mywcs.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'FREQ'
    mywcs.wcs.cd = np.ones((3, 3))
    mywcs.wcs.set()
    frame = wcs_to_celestial_frame(mywcs)
    assert isinstance(frame, ICRS) 

def _celestial_frame_to_wcs_builtin(frame, projection='TAN'):

    # Import astropy.coordinates here to avoid circular imports
    from astropy.coordinates import BaseRADecFrame, FK4, FK4NoETerms, FK5, ICRS, ITRS, Galactic

    # Create a 2-dimensional WCS
    wcs = WCS(naxis=2)

    if isinstance(frame, BaseRADecFrame):

        xcoord = 'RA--'
        ycoord = 'DEC-'
        if isinstance(frame, ICRS):
            wcs.wcs.radesys = 'ICRS'
        elif isinstance(frame, FK4NoETerms):
            wcs.wcs.radesys = 'FK4-NO-E'
            wcs.wcs.equinox = frame.equinox.byear
        elif isinstance(frame, FK4):
            wcs.wcs.radesys = 'FK4'
            wcs.wcs.equinox = frame.equinox.byear
        elif isinstance(frame, FK5):
            wcs.wcs.radesys = 'FK5'
            wcs.wcs.equinox = frame.equinox.jyear
        else:
            return None
    elif isinstance(frame, Galactic):
        xcoord = 'GLON'
        ycoord = 'GLAT'
    elif isinstance(frame, ITRS):
        xcoord = 'TLON'
        ycoord = 'TLAT'
        wcs.wcs.radesys = 'ITRS'
        wcs.wcs.dateobs = frame.obstime.utc.isot
    else:
        return None

    wcs.wcs.ctype = [xcoord + '-' + projection, ycoord + '-' + projection]

    return wcs 

def test_caching_components_and_classes():

    # Make sure that when we change the WCS object, the classes and components
    # are updated (we use a cache internally, so we need to make sure the cache
    # is invalidated if needed)

    wcs = WCS_SIMPLE_CELESTIAL

    assert wcs.world_axis_object_components == [('celestial', 0, 'spherical.lon.degree'),
                                                ('celestial', 1, 'spherical.lat.degree')]

    assert wcs.world_axis_object_classes['celestial'][0] is SkyCoord
    assert wcs.world_axis_object_classes['celestial'][1] == ()
    assert isinstance(wcs.world_axis_object_classes['celestial'][2]['frame'], ICRS)
    assert wcs.world_axis_object_classes['celestial'][2]['unit'] is u.deg

    wcs.wcs.radesys = 'FK5'

    frame = wcs.world_axis_object_classes['celestial'][2]['frame']
    assert isinstance(frame, FK5)
    assert frame.equinox.jyear == 2000.

    wcs.wcs.equinox = 2010

    frame = wcs.world_axis_object_classes['celestial'][2]['frame']
    assert isinstance(frame, FK5)
    assert frame.equinox.jyear == 2010. 

def correct_vel(ra, dec, vel, vlsr=vlsr0):
    """Corrects the proper motion for the speed of the Sun
    Arguments:
        ra - RA in deg
        dec -- Declination in deg
        pmra -- pm in RA in mas/yr
        pmdec -- pm in declination in mas/yr
        dist -- distance in kpc
    Returns:
        (pmra,pmdec) the tuple with the proper motions corrected for the Sun's motion
    """
    
    C=acoo.ICRS(ra=ra*auni.deg,dec=dec*auni.deg,
                    radial_velocity=vel*auni.km/auni.s,
                    distance=np.ones_like(vel)*auni.kpc,
                    pm_ra_cosdec=np.zeros_like(vel)*auni.mas/auni.year,
                    pm_dec=np.zeros_like(vel)*auni.mas/auni.year)
    #frame = acoo.Galactocentric (galcen_vsun = np.array([ 11.1, vlsr+12.24, 7.25])*auni.km/auni.s)
    kw = dict(galcen_v_sun = acoo.CartesianDifferential(np.array([ 11.1, vlsr+12.24, 7.25])*auni.km/auni.s))                                 
    frame = acoo.Galactocentric (**kw)
    Cg = C.transform_to(frame)
    Cg1 = acoo.Galactocentric(x=Cg.x, y=Cg.y, z=Cg.z,
                              v_x=Cg.v_x*0,
                              v_y=Cg.v_y*0,
                              v_z=Cg.v_z*0, **kw)
    C1=Cg1.transform_to(acoo.ICRS)
    return np.asarray(((C.radial_velocity-C1.radial_velocity)/(auni.km/auni.s)).decompose()) 

def set_frame(self, new_frame):
        """
            Set a new frame for the coordinates (the default is ICRS J2000)
            
            :param new_frame: a coordinate frame from astropy
            :return: (none)
            """
        assert new_frame.lower() in ['icrs', 'galactic', 'fk5', 'fk4', 'fk4_no_e' ]
                
        self._frame = new_frame 

def set_frame(self, new_frame):
        """
        Set a new frame for the coordinates (the default is ICRS J2000)

        :param new_frame: a coordinate frame from astropy
        :return: (none)
        """
        assert isinstance(new_frame, BaseCoordinateFrame)

        self._frame = new_frame 

def _setup(self):

        self._frame = ICRS() 

def azel2radec(
    az_deg: "ndarray", el_deg: "ndarray", lat_deg: "ndarray", lon_deg: "ndarray", time: datetime, *, use_astropy: bool = True
) -> typing.Tuple["ndarray", "ndarray"]:
    """
    viewing angle (az, el) to sky coordinates (ra, dec)

    Parameters
    ----------
    az_deg : "ndarray"
         azimuth [degrees clockwize from North]
    el_deg : "ndarray"
             elevation [degrees above horizon (neglecting aberration)]
    lat_deg : "ndarray"
              observer latitude [-90, 90]
    lon_deg : "ndarray"
              observer longitude [-180, 180] (degrees)
    time : datetime.datetime or str
           time of observation
    use_astropy : bool, optional
                 default use astropy.

    Returns
    -------
    ra_deg : "ndarray"
         ecliptic right ascension (degress)
    dec_deg : "ndarray"
         ecliptic declination (degrees)
    """

    if use_astropy and Time is not None:

        obs = EarthLocation(lat=lat_deg * u.deg, lon=lon_deg * u.deg)

        direc = AltAz(location=obs, obstime=Time(str2dt(time)), az=az_deg * u.deg, alt=el_deg * u.deg)

        sky = SkyCoord(direc.transform_to(ICRS()))

        return sky.ra.deg, sky.dec.deg

    return vazel2radec(az_deg, el_deg, lat_deg, lon_deg, time) 

def strip_telescope_header(header, simplify_wcs=True):
    """
    Strip non-JWST keywords that confuse `jwst.datamodels.util.open`.

    Parameters
    ----------
    header : `~astropy.io.fits.Header`
        Input FITS header.

    """
    import astropy.wcs as pywcs

    new_header = header.copy()

    if 'TELESCOP' in new_header:
        if new_header['TELESCOP'] != 'JWST':
            keys = ['TELESCOP', 'FILTER', 'DETECTOR', 'INSTRUME']
            for key in keys:
                if key in header:
                    new_header.remove(key)

    if simplify_wcs:
        # Make simple WCS header
        orig_wcs = pywcs.WCS(new_header)
        new_header = orig_wcs.to_header()

        new_header['EXTNAME'] = 'SCI'
        new_header['RADESYS'] = 'ICRS'
        new_header['CDELT1'] = -new_header['PC1_1']
        new_header['CDELT2'] = new_header['PC2_2']
        new_header['PC1_1'] = -1
        new_header['PC2_2'] = 1

    return new_header 

def test_unit_representation_subclass():

    class Longitude180(Longitude):
        def __new__(cls, angle, unit=None, wrap_angle=180*u.deg, **kwargs):
            self = super().__new__(cls, angle, unit=unit, wrap_angle=wrap_angle,
                                   **kwargs)
            return self

    class UnitSphericalWrap180Representation(UnitSphericalRepresentation):
        attr_classes = OrderedDict([('lon', Longitude180),
                                    ('lat', Latitude)])

    class SphericalWrap180Representation(SphericalRepresentation):
        attr_classes = OrderedDict([('lon', Longitude180),
                                    ('lat', Latitude),
                                    ('distance', u.Quantity)])

        _unit_representation = UnitSphericalWrap180Representation

    class MyFrame(ICRS):
        default_representation = SphericalWrap180Representation
        frame_specific_representation_info = {
            'spherical': [
                RepresentationMapping('lon', 'ra'),
                RepresentationMapping('lat', 'dec')]
        }
        frame_specific_representation_info['unitsphericalwrap180'] = \
            frame_specific_representation_info['sphericalwrap180'] = \
            frame_specific_representation_info['spherical']

    @frame_transform_graph.transform(FunctionTransform,
                                     MyFrame, astropy.coordinates.ICRS)
    def myframe_to_icrs(myframe_coo, icrs):
        return icrs.realize_frame(myframe_coo._data)

    f = MyFrame(10*u.deg, 10*u.deg)
    assert isinstance(f._data, UnitSphericalWrap180Representation)
    assert isinstance(f.ra, Longitude180)

    g = f.transform_to(astropy.coordinates.ICRS)
    assert isinstance(g, astropy.coordinates.ICRS)
    assert isinstance(g._data, UnitSphericalWrap180Representation)

    frame_transform_graph.remove_transform(MyFrame,
                                           astropy.coordinates.ICRS,
                                           None) 

def test_kdtree_storage(functocheck, args, defaultkdtname, bothsaved):
    from astropy.coordinates import ICRS

    def make_scs():
        cmatch = ICRS([4, 2.1]*u.degree, [0, 0]*u.degree, distance=[1, 2]*u.kpc)
        ccatalog = ICRS([1, 2, 3, 4]*u.degree, [0, 0, 0, 0]*u.degree, distance=[1, 2, 3, 4]*u.kpc)
        return cmatch, ccatalog

    cmatch, ccatalog = make_scs()
    functocheck(cmatch, ccatalog, *args, storekdtree=False)
    assert 'kdtree' not in ccatalog.cache
    assert defaultkdtname not in ccatalog.cache

    cmatch, ccatalog = make_scs()
    functocheck(cmatch, ccatalog, *args)
    assert defaultkdtname in ccatalog.cache
    assert 'kdtree' not in ccatalog.cache

    cmatch, ccatalog = make_scs()
    functocheck(cmatch, ccatalog, *args, storekdtree=True)
    assert 'kdtree' in ccatalog.cache
    assert defaultkdtname not in ccatalog.cache

    cmatch, ccatalog = make_scs()
    assert 'tislit_cheese' not in ccatalog.cache
    functocheck(cmatch, ccatalog, *args, storekdtree='tislit_cheese')
    assert 'tislit_cheese' in ccatalog.cache
    assert defaultkdtname not in ccatalog.cache
    assert 'kdtree' not in ccatalog.cache
    if bothsaved:
        assert 'tislit_cheese' in cmatch.cache
        assert defaultkdtname not in cmatch.cache
        assert 'kdtree' not in cmatch.cache
    else:
        assert 'tislit_cheese' not in cmatch.cache

    # now a bit of a hacky trick to make sure it at least tries to *use* it
    ccatalog.cache['tislit_cheese'] = 1
    cmatch.cache['tislit_cheese'] = 1
    with pytest.raises(TypeError) as e:
        functocheck(cmatch, ccatalog, *args, storekdtree='tislit_cheese')
    assert 'KD' in e.value.args[0] 

def _wcs_to_celestial_frame_builtin(wcs):

    # Import astropy.coordinates here to avoid circular imports
    from astropy.coordinates import (FK4, FK4NoETerms, FK5, ICRS, ITRS,
                                     Galactic, SphericalRepresentation)

    # Import astropy.time here otherwise setup.py fails before extensions are compiled
    from astropy.time import Time

    if wcs.wcs.lng == -1 or wcs.wcs.lat == -1:
        return None

    radesys = wcs.wcs.radesys

    if np.isnan(wcs.wcs.equinox):
        equinox = None
    else:
        equinox = wcs.wcs.equinox

    xcoord = wcs.wcs.ctype[wcs.wcs.lng][:4]
    ycoord = wcs.wcs.ctype[wcs.wcs.lat][:4]

    # Apply logic from FITS standard to determine the default radesys
    if radesys == '' and xcoord == 'RA--' and ycoord == 'DEC-':
        if equinox is None:
            radesys = "ICRS"
        elif equinox < 1984.:
            radesys = "FK4"
        else:
            radesys = "FK5"

    if radesys == 'FK4':
        if equinox is not None:
            equinox = Time(equinox, format='byear')
        frame = FK4(equinox=equinox)
    elif radesys == 'FK4-NO-E':
        if equinox is not None:
            equinox = Time(equinox, format='byear')
        frame = FK4NoETerms(equinox=equinox)
    elif radesys == 'FK5':
        if equinox is not None:
            equinox = Time(equinox, format='jyear')
        frame = FK5(equinox=equinox)
    elif radesys == 'ICRS':
        frame = ICRS()
    else:
        if xcoord == 'GLON' and ycoord == 'GLAT':
            frame = Galactic()
        elif xcoord == 'TLON' and ycoord == 'TLAT':
            # The default representation for ITRS is cartesian, but for WCS
            # purposes, we need the spherical representation.
            frame = ITRS(representation_type=SphericalRepresentation,
                         obstime=wcs.wcs.dateobs or None)
        else:
            frame = None

    return frame 

def hdu_to_imagemodel(in_hdu):
    """
    Workaround for initializing a `jwst.datamodels.ImageModel` from a
    normal FITS ImageHDU that could contain HST header keywords and
    unexpected WCS definition.

    TBD

    Parameters
    ----------
    in_hdu : `astropy.io.fits.ImageHDU`


    Returns
    -------
    img : `jwst.datamodels.ImageModel`


    """
    from astropy.io.fits import ImageHDU, HDUList
    from astropy.coordinates import ICRS

    from jwst.datamodels import util
    import gwcs

    hdu = ImageHDU(data=in_hdu.data, header=in_hdu.header)

    new_header = strip_telescope_header(hdu.header)

    hdu.header = new_header

    # Initialize data model
    img = util.open(HDUList([hdu]))

    # Initialize GWCS
    tform = gwcs.wcs.utils.make_fitswcs_transform(new_header)
    hwcs = gwcs.WCS(forward_transform=tform, output_frame=ICRS())  # gwcs.CelestialFrame())
    sh = hdu.data.shape
    hwcs.bounding_box = ((-0.5, sh[0]-0.5), (-0.5, sh[1]-0.5))

    # Put gWCS in meta, where blot/drizzle expect to find it
    img.meta.wcs = hwcs

    return img