simulating HST WCS

hst_wcs.py

# Simulate WFC3 WCS
import numpy as np
from astropy.io import fits
from stwcs import updatewcs

primary_header_kw = {'PA_V3': 81.872513, 'instrume': 'WFC3', 'detector': 'IR',
                     'idctab': 'w3m18525i_idc.fits', 'filter': 'F160W',
                     'telescop': 'HST', 'date-obs': '2009-08-07'}

sci_header_kw = {'CRPIX1': 562, 'CRPIX2': 562,
                 'CRVAL1': 47.3705 , 'CRVAL2': 61.59299543,
                 'CTYPE1': 'RA---TAN', 'CTYPE2': 'DEC--TAN', 'RADESYS': 'ICRS',
                 'inherit': True, 'EXTNAME': 'SCI'}

hdul = fits.HDUList()
phdu = fits.PrimaryHDU()
data = np.zeros((1024, 1024))
scihdu1 = fits.ImageHDU(data=data)
scihdu2 = fits.ImageHDU(data=data)

for key, val in sci_header_kw.items():
    scihdu1.header[key] = val

for key, val in sci_header_kw.items():
    scihdu2.header[key] = val

for key, val in primary_header_kw.items():
    phdu.header[key] = val

scihdu1.header['extver'] = 1
scihdu2.header['extver'] = 2

plate_scale = 0.12825
pav3 = np.deg2rad(primary_header_kw['PA_V3'])
cd = np.array([[np.cos(pav3), np.sin(pav3)], [-np.sin(pav3), np.cos(pav3)]])
cd /= 3600 / plate_scale

scihdu1.header['CD1_1'] = cd[0, 0]
scihdu1.header['CD1_2'] = cd[0, 1]
scihdu1.header['CD2_1'] = cd[1, 0]
scihdu1.header['CD2_2'] = cd[1, 1]

scihdu2.header['CD1_1'] = cd[0, 0]
scihdu2.header['CD1_2'] = cd[0, 1]
scihdu2.header['CD2_1'] = cd[1, 0]
scihdu2.header['CD2_2'] = cd[1, 1]


hdul.append(phdu)
hdul.append(scihdu1)
hdul.append(scihdu2)
hdul.writeto('wfc3.fits', clobber=True)

# Update the WCS  makes small corrections to the primary WCS and attaches SIP
from stwcs import updatewcs
from stwcs.wcsutil import HSTWCS
updatewcs.updatewcs('wfc3.fits', checkfiles=False)

# Look at the wcs
w1 = HSTWCS('wfc3.fits', ext=1)
w1.printwcs()