Notebook

Useful Functionalities¶

Since Astropy is not a successor of IRAF, you may not easily find some useful functionalites in astropy-related packages, which were just a "basic" set of IRAF functions. These include IMCOMBINE with the option offset=wcs.

Image Combine with Offset¶

In ccdproc, there is a functionality to combine images, such as function combine. You can also do the combiner as in the example of a section in ccdproc official manual. A more friendly way of expressing this is:

ccd1 = CCDData.read('fits/fits1.fits', unit='adu')
ccd2 = CCDData.read('fits/fits2.fits', unit='adu')
ccd3 = CCDData.read('fits/fits3.fits', unit='adu')
comb = Combiner([ccd1, ccd2, ccd3])

The median/average/etc combine can be done as in the manual, and the clipped mask array can also be generated (see the manual). One thing to stress is the "WCS offset combination", which is imcombine img1 img2 output=img_comb offset=wcs in IRAF:

from ccdproc import wcs_project
from astropy.wcs import WCS
wcs_origin = WCS(ccd1.header)
reprojected = []
for img in my_list_of_images:
    new_image = wcs_project(img, target_wcs)
    reprojected.append(new_image)

The following shows a more practical example with four non-sidereal mode tracking images taken from Okayama Astrophysical Observatory:

In [3]:

import numpy as np
from ccdproc import CCDData, wcs_project, combine
from astropy.wcs import WCS

# To silence the VerifyWarning, do the following:
from astropy.utils.exceptions import AstropyWarning
import warnings
warnings.simplefilter('ignore', category=AstropyWarning)

#%%
# Load filename list
prefix = 'reproj_Tutorial/'
filelist = np.loadtxt(prefix+'fits.list', dtype=bytes).astype(str)

# make an empty list for original and reprojected images
before = []
reproj = []

# Specify the "target" wcs that other images to be re-projected
hdr_std = CCDData.read(prefix+filelist[0], unit='adu').header
wcs_std = WCS(hdr_std)

# Re-project all the images
for fname in filelist:
    ccd = CCDData.read(prefix+fname, unit='adu')
    before.append(ccd)
    reproj.append(wcs_project(ccd, wcs_std))

# Average combine and save
avg = combine(reproj, output_file='test_avg.fits', method='average')

In the following, I showd the original image which was used for the target wcs setting (image of index 0), the ratio between this original image and average combined image, and the average combined imgage for comparison.

In [10]:

from matplotlib import pyplot as plt

fig, ax = plt.subplots(1, 3, figsize=(12, 6))
im0 = ax[0].imshow(before[0].data, vmin=250, vmax=500, origin='lower')
im1 = ax[1].imshow(before[0].data/avg.data, vmin=0.9, vmax=1.2, origin='lower')
im2 = ax[2].imshow(avg.data, vmin=250, vmax=500, origin='lower')

ax[0].set_title('image 1')
ax[1].set_title('image 1 / avg_combined\n(image1 = wcs standard)', y=1.05)
ax[2].set_title('avg_combined')

plt.setp(ax[1].get_yticklabels(), visible=False)
plt.setp(ax[2].get_yticklabels(), visible=False)
plt.setp((ax[0], ax[1], ax[2]),
         xlim=[0, 1024], ylim=[0, 1024],
         xticks=np.arange(0, 1024, 250),
         yticks=np.arange(0, 1024, 250))

cbar_ax0 = fig.add_axes([0.12, 0.15, 0.23, 0.02]) # left, bottom, width, height
cbar_ax1 = fig.add_axes([0.40, 0.15, 0.23, 0.02]) # left, bottom, width, height
cbar_ax2 = fig.add_axes([0.68, 0.15, 0.23, 0.02]) # left, bottom, width, height
fig.colorbar(im0, cax=cbar_ax0, ticks=np.arange(250, 501, 50),
             orientation='horizontal', label='ADU')
fig.colorbar(im1, cax=cbar_ax1, ticks=np.arange(0.9, 1.21, 0.1),
             orientation='horizontal',
             label='ratio')
fig.colorbar(im2, cax=cbar_ax2, ticks=np.arange(250, 501, 50),
             orientation='horizontal',
             label='ADU')

plt.show()

Question: What will happen when you do median combine? These four are the non-sidereal tracking mode image, so you should be able to find some point-like asteroid near the center of each image. But does that asteroid visible in median combined image? Why or why not?
Question: Why do you think the middle column of the image has values exactly 1 for pixels $y \gtrsim 950$ ?
TIP: You may notice that the result is different from IRAF IMCOMBINE task. I think this is because ccdproc has not yet fully implemented such functionality, and I asked the developers for that. Let's be tay tuned for future improvements.

Query Stars (`astroquery`)¶

The astroquery version 0.3.5 documentation explains all the basic ideas of how to use astroquery. Installation can be done by:

conda install -c astropy astroquery

An example is as follows:

What I want to do:

Query the URAT1 catalog stars which contains the APASS r-band magnitude
- The stars should have "exact match" with APASS (mfa=1).
- The magnitude should be between 14 to 16 with error between 0.0 to 0.15 magnitude.
- The names of the columns and their meanings are explained here (Type "URAT" and search)
- I want to query all the available stars.

Code:

from astroquery.vizier import Vizier as V
from astropy import units as u
from astropy.io import fits

# Prepare for quearyinig Vizier with filters
# up to infinitely many rows. By default, this is 50.
viz_filt = V(columns=['RAJ2000', 'DEJ2000', 'mfa', 'rmag', 'e_rmag'] ,
             column_filters={'mfa':'=1', 
                             'rmag':'{0} .. {1}'.format(14, 16),
                             'e_rmag':'0.0 .. 0.15'})
viz_filt.ROW_LIMIT = -1

# Load FITS image
hdul = fits.open('fits/image.fits')
header = hdul[0].header
coords = SkyCoord(header['RA'] + header['DEC'], unit=(u.hourangle, u.deg))

# Query!
result = viz_filt.query_region(coords, 
                               radius=10*u.arcmin, 
                               catalog=["URAT1"])
# Save data
data = result['I/329/urat1']

Query SSSBs (Asteroids / Comets)¶

There is a useful module developed by M. Mommert here. You can install on terminal via

pip install callhorizons

and import it on python via

import callhorizons

An example is

import callhorizons
from astropy.io import fits
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import astropy.units as u
#%%
# make filelist and load filelist
filelist= np.loadtxt('QY1fits.list', dtype=bytes).astype(str)

# open fits as hdu (header-data unit) list
hdulist = fits.open(filelist[0])
header  = hdulist[0].header     # zero-th extension for non MEF

# Set time for query
print(header['date-obs'])  # It's in isot in my case, e.g., 2016-08-02T20:23:56.61
time = Time(header['date-obs'], format='isot')
time +=  0.5 * header['exptime'] * u.s  # add half of exposure (seconds)

# Query! We need: (1) target name, (2) time, (3) observatory
dq = callhorizons.query('331471')
dq.set_discreteepochs(time.jd)   # input discrete times in jd unit
dq.get_ephemerides(observatory_code=511) # observatory code
# Many more options available such as skip_daylight.

# Check:
print(dq.fields) # gives all available data in "dq"
print(dq['RA'])  # gives the RA of the target in degrees

# save RA and DEC:
RA  = dq['RA'][0]  # We had only one discrete epochs, so [0] is needed.
DEC = dq['DEC'][0] # These are in degrees.

# Set wcs and change RA, DEC into image XY:
wcs = WCS(header)
XY = SkyCoord(RA, DEC, unit='deg').to_pixel(wcs)
X  = str( round(float(XY[0]), 2) ) 
Y  = str( round(float(XY[1]), 2) )
print('ImageXY = ({0}, {1})'.format(X, Y))
# ImageXY = (394.69, 607.56)

Useful Functionalities¶

Image Combine with Offset¶

Query Stars (astroquery)¶

Query SSSBs (Asteroids / Comets)¶

Extract Contour¶

Query Stars (`astroquery`)¶