How to use the jwst.datamodels.util function in jwst

Anything `jwst.datamodels.util.open` can accept for initialization.

    Returns
    -------
    with_wcs : `jwst.datamodels.ImageModel`
        Image model with full `~gwcs` in `with_wcs.meta.wcs`.

    """
    from jwst.datamodels import util
    from jwst.assign_wcs import AssignWcsStep

    # from jwst.stpipe import crds_client
    # from jwst.assign_wcs import assign_wcs

    # HDUList -> jwst.datamodels.ImageModel
    img = util.open(input)

    # AssignWcs to pupulate img.meta.wcsinfo
    step = AssignWcsStep()
    with_wcs = step.process(img)

    # Above should be more robust to get all of the necessary ref files
    #dist_file = crds_client.get_reference_file(img, 'distortion')
    #reference_files = {'distortion': dist_file}
    #with_wcs = assign_wcs.load_wcs(img, reference_files=reference_files)

    return with_wcs

"""
    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