FITS WCS rotation angles

crota.py

#!/usr/bin/env python3
import sys
import numpy as np
from astropy.wcs import wcs

w = wcs.WCS(sys.argv[1])

sx = np.sqrt(w.wcs.cd[0, 0]**2+w.wcs.cd[0, 1]**2)*3600.0
sy = np.sqrt(w.wcs.cd[1, 0]**2+w.wcs.cd[1, 1]**2)*3600.0
ang1 = np.arctan2(w.wcs.cd[0, 1], w.wcs.cd[0, 0])*180.0/np.pi
ang2 = np.arctan2(-w.wcs.cd[1, 0], w.wcs.cd[1, 1])*180.0/np.pi

print("scale: %gx%g arcsec, rotation: %g, %g deg" %(sx, sy, ang1, ang2))