Skip to content

Instantly share code, notes, and snippets.

@mwcraig
Created November 24, 2019 22:04
Show Gist options
  • Save mwcraig/33b9fab4edbffefdeac4a63be4b60b5b to your computer and use it in GitHub Desktop.
Save mwcraig/33b9fab4edbffefdeac4a63be4b60b5b to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from astropy.io import fits\n",
"from astropy.wcs import WCS\n",
"from astropy.nddata import CCDData"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## This notebook is checking what astropy does when `PC` and `CD` are included in a FITS header\n",
"\n",
"The file below has both."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"bad_file = 'TIC_83210867.01-S001-R089-C001-rp-as_pc.new'"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"with fits.open(bad_file) as f:\n",
" raw_bad = f[0]"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def print_cd_pc(header):\n",
" \"\"\"\n",
" Print which of PC and CD are in the header.\n",
" \"\"\"\n",
" cd_pc = [k for k in header.keys() if k.startswith('CD') or k.startswith('PC')]\n",
" for k in sorted(cd_pc):\n",
" print(f'{k:6} = {header[k]:e}')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CD1_1 = 1.522566e-04\n",
"CD1_2 = -3.573089e-05\n",
"CD2_1 = 3.575750e-05\n",
"CD2_2 = 1.522348e-04\n",
"PC1_1 = 1.522621e-04\n",
"PC1_2 = -3.564429e-05\n",
"PC2_1 = 3.566437e-05\n",
"PC2_2 = 1.521763e-04\n"
]
}
],
"source": [
"print_cd_pc(raw_bad.header)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yep, both are there and they are not equal to each other. If they were, there would be no isue."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Read this file as a `CCDData`, which will convert the header to a WCS."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
"the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
]
}
],
"source": [
"ccdd = CCDData.read(bad_file, format='fits')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is nothing at all special about this test point except that I know it is not `CRPIX` in this header."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"test_point = [4000, 4000]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"# This will be handy for deleting either CD or PC in a moment\n",
"CDs = [k for k in raw_bad.header if k.startswith('CD')]\n",
"PCs = [k for k in raw_bad.header if k.startswith('PC')]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make a wcs without PC and confirm that it has no PC"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CD1_1 = 1.522566e-04\n",
"CD1_2 = -3.573089e-05\n",
"CD2_1 = 3.575750e-05\n",
"CD2_2 = 1.522348e-04\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
"the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
]
},
{
"data": {
"text/plain": [
"WCS Keywords\n",
"\n",
"Number of WCS axes: 2\n",
"CTYPE : 'RA---TAN-SIP' 'DEC--TAN-SIP' \n",
"CRVAL : 297.459006197 34.8301203697 \n",
"CRPIX : 2055.0 2048.5 \n",
"CD1_1 CD1_2 : 0.000152256648923 -3.57308903932e-05 \n",
"CD2_1 CD2_2 : 3.57575046574e-05 0.000152234782375 \n",
"NAXIS : 4109 4096"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"head_no_pcs = raw_bad.header.copy()\n",
"for pc in PCs:\n",
" del head_no_pcs[pc]\n",
"\n",
"print_cd_pc(head_no_pcs)\n",
"wcs_no_pcs = WCS(head_no_pcs)\n",
"wcs_no_pcs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make a wcs without CD and confirm that it has no CD"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PC1_1 = 1.522621e-04\n",
"PC1_2 = -3.564429e-05\n",
"PC2_1 = 3.566437e-05\n",
"PC2_2 = 1.521763e-04\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
"the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
]
},
{
"data": {
"text/plain": [
"WCS Keywords\n",
"\n",
"Number of WCS axes: 2\n",
"CTYPE : 'RA---TAN-SIP' 'DEC--TAN-SIP' \n",
"CRVAL : 297.459006197 34.8301203697 \n",
"CRPIX : 2055.0 2048.5 \n",
"PC1_1 PC1_2 : 0.000152262101123 -3.56442880876e-05 \n",
"PC2_1 PC2_2 : 3.56643743388e-05 0.000152176346785 \n",
"CDELT : 1.0 1.0 \n",
"NAXIS : 4109 4096"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"head_no_cds = raw_bad.header.copy()\n",
"for cd in CDs:\n",
" del head_no_cds[cd]\n",
"\n",
"print_cd_pc(head_no_cds)\n",
"wcs_no_cds = WCS(head_no_cds)\n",
"wcs_no_cds"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compare the CCDData-generated coordinates with these two headers\n",
"\n",
"This comparison ignores the SIP distortions because I want to look at just the transform, not the distortion."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WCS from CCDData [array(297.73642431), array(35.19632622)]\n",
"WCS no CD, just PC [array(297.73642431), array(35.19632622)]\n",
"WCS no PC, just CD [array(297.73620541), array(35.19662204)]\n"
]
}
],
"source": [
"print('WCS from CCDData', ccdd.wcs.wcs_pix2world(*test_point, 0))\n",
"print('WCS no CD, just PC', wcs_no_cds.wcs_pix2world(*test_point, 0))\n",
"print('WCS no PC, just CD', wcs_no_pcs.wcs_pix2world(*test_point, 0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conclusion: astropy uses PC, not CD, if both are present"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment