<jupyter_start><jupyter_text>The image viewer widget is optional...<jupyter_code>iw = ImageWidget()
iw
iw.load_fits('kelt-16-S001-R001-C100-r.fit')
#iw.load_fits('kelt-16-b-S001-R001-C130-r.fit')
iw.center_on(SkyCoord(313.92459686592696, 31.78329193408578, unit='degree'), pixel_coords_offset=0)
iw.center_on(SkyCoord.from_name('kelt-16'))<jupyter_output><empty_output><jupyter_text>Required stuff starts here  Read one image to use detect sourcesThis could be any image; I chose one roughly in the middle of the sequence. Ideally maybe we'd choose sources that were only in all of the images, but this will work for now.<jupyter_code>ccd = CCDData.read('kelt-16-S001-R001-C100-r.fit')
def faster_sigma_clip_stats(data, sigma=5, iters=5, axis=None):
    data = data.copy()
    for _ in range(iters):
        central = bn.nanmedian(data, axis=axis)
        try:
            central = central[:, np.newaxis]
        except (ValueError, IndexError):
            pass
                
        std_dif = 1.4826 * bn.nanmedian(np.abs(data - central))
        clips = np.abs(data - central) / std_dif > sigma
        if clips.sum() == 0:
            break
        data[clips] = np.nan
    return bn.nanmean(data, axis=axis), bn.nanmedian(data, axis=axis), bn.nanstd(data, axis=axis)<jupyter_output><empty_output><jupyter_text>The background light from the sky needs to be removed before looking for sources. We'll do that by calculating the median. "Sigma clipping" is the term for excluding data that is far from the average...here we do it mostly to exclude the "bright" pixels (i.e. the stars) from our estimate of the background.<jupyter_code>men, med, std = faster_sigma_clip_stats(ccd.data, sigma=5)
men, med, std<jupyter_output>C:\Anaconda3\envs\ast366\lib\site-packages\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in greater
  # This is added back by InteractiveShellApp.init_path()
C:\Anaconda3\envs\ast366\lib\site-packages\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in greater
  # This is added back by InteractiveShellApp.init_path()
C:\Anaconda3\envs\ast366\lib\site-packages\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in greater
  # This is added back by InteractiveShellApp.init_path()<jupyter_text>This sets up the source detection. The FWHM turns out to be key...making it too small results in a single star being detected as two separate sources.Stars must be brighter than the threshold to count as sources. Making the number higher gives you fewer detected sources, lower gives you more. There is no "magic" number.<jupyter_code>dao = DAOStarFinder(threshold=10 * std, fwhm=8, exclude_border=True)<jupyter_output><empty_output><jupyter_text>Actually detect the stars...<jupyter_code>stars = dao(ccd - med)<jupyter_output><empty_output><jupyter_text>Only include this if you have an image viewer aboveIt is handy to see where the detected sources are....<jupyter_code>iw.reset_markers()
iw.marker = {'type': 'circle', 'color': 'lightgreen', 'radius': 10}
iw.add_markers(stars, x_colname='xcentroid', y_colname='ycentroid', pixel_coords_offset=0)<jupyter_output><empty_output><jupyter_text>Pick aperture and annulus sizesProbably the best way to do this for now is to upload one of the kelt-16 images to the notebook we used in class to plot star profiles. You want an aperture about 2×FWHM, an annulus with  an inner radius at least 10 pixels larger than that, and outer radius about 15 pixels larger than the inner.<jupyter_code>aperture_rad = 12
inner_annulus = aperture_rad + 15
outer_annulus = inner_annulus + 15<jupyter_output><empty_output><jupyter_text>Set up the aperture objects used later by photutils to do the photometry.<jupyter_code>aps = CircularAperture([stars['xcentroid'], stars['ycentroid']], r=aperture_rad)
anuls = CircularAnnulus([stars['xcentroid'], stars['ycentroid']], inner_annulus, outer_annulus)<jupyter_output><empty_output><jupyter_text>Set pixels that are above point where CCD becomes non-linear to invalid valueThis will ensure that the aperture sum is an invalid value if one or more of the pixels in the aperture is non-linear. *Change the value below to what you think is appropriate -- I'm going from memory.*<jupyter_code>max_adu = 45000
ccd.data[ccd.data > max_adu] = np.nan<jupyter_output><empty_output><jupyter_text>Set up sky coordinates of target listWe'll want to the list of stars in RA/Dec so we can find them in each image. `star_locs` is waht is used later on to place the apertures<jupyter_code>star_locs = ccd.wcs.all_pix2world(stars['xcentroid'], stars['ycentroid'], 0)<jupyter_output><empty_output><jupyter_text>Check for stars closer than 2 × aperture radius and remove from list<jupyter_code>star_coords = SkyCoord(ra=star_locs[0], dec=star_locs[1], frame='icrs', unit='degree')
idxc, d2d, d3d = star_coords.match_to_catalog_sky(star_coords, nthneighbor=2)
d2d.min()
too_close = d2d < (aperture_rad * 2 * 0.563 * u.arcsec)<jupyter_output><empty_output><jupyter_text>Include this only if you have the image viewer -- it shows which stars are being removed<jupyter_code>iw.marker = {'type': 'circle', 'color': 'red', 'radius': 20}
iw.add_markers(stars[too_close], x_colname='xcentroid', y_colname='ycentroid', pixel_coords_offset=0)
star_locs = (star_locs[0][~too_close], star_locs[1][~too_close])<jupyter_output><empty_output><jupyter_text>Make up a unique ID for each star<jupyter_code>star_ids = np.arange(len(star_locs[0])) + 1<jupyter_output><empty_output><jupyter_text>Define a couple of convenience functions  Calculate average pixel values in annulus, rejecting extreme values<jupyter_code>def clipped_sky_per_pix_stats(data, annulus, sigma=5, iters=5):
    # Get a list of masks from the annuli
    # Use the 'center' method because then pixels are either in or out. To use
    # 'partial' or 'exact' we would need to do a weighted sigma clip and I'm not sure 
    # how to do that.
    masks = annulus.to_mask(method='center')
    
    anul_list = []
    for mask in masks:
        # Multiply the mask times the data
        to_clip = mask.multiply(data.data, fill_value=np.nan)
        anul_list.append(to_clip.flatten())
    # Convert the list to an array for doing the sigma clip
    anul_array = np.array(anul_list)
    # Turn all zeros into np.nan...
    anul_array[anul_array == 0] = np.nan
    avg_sky_per_pix, med_sky_per_pix, std_sky_per_pix = faster_sigma_clip_stats(anul_array, 
                                                                                sigma=sigma, 
                                                                                iters=iters,
                                                                                axis=1
                                                                               )

    return (avg_sky_per_pix * data.unit, med_sky_per_pix * data.unit, std_sky_per_pix * data.unit)<jupyter_output><empty_output><jupyter_text>Add more columns to the data table<jupyter_code>def add_to_photometry_table(phot, ccd, annulus, apertures, fname='', 
                            star_ids=None, gain=None):
    phot.rename_column('aperture_sum_0', 'aperture_sum')
    phot['aperture_sum'].unit = u.adu
    phot.rename_column('aperture_sum_1', 'annulus_sum')
    star_locs = ccd.wcs.all_pix2world(phot['xcenter'], phot['ycenter'], 0)
    star_coords = SkyCoord(ra=star_locs[0], dec=star_locs[1], frame='icrs', unit='degree')
    phot['RA'] = star_coords.ra
    phot['Dec'] = star_coords.dec
    print('        ...calculating clipped sky stats')
    avg_sky_per_pix, med_sky_per_pix, std_sky_per_pix = clipped_sky_per_pix_stats(ccd, annulus)
    print('        ...DONE calculating clipp sky stats')
    phot['sky_per_pix_avg'] = avg_sky_per_pix
    phot['sky_per_pix_med'] = med_sky_per_pix
    phot['sky_per_pix_std'] = std_sky_per_pix
    phot['aperture'] = apertures.r * u.pixel
    phot['aperture_area'] = apertures.area() # * u.pixel * u.pixel
    phot['annulus_inner'] = annulus.r_in * u.pixel
    phot['annulus_outer'] = annulus.r_out * u.pixel
    phot['annulus_area'] = annulus.area() #* u.pixel * u.pixel
    phot['exposure'] = [ccd.header['exposure']] * len(phot) * u.second
    phot['date-obs'] = [ccd.header['DATE-OBS']] * len(phot)
    night = Time(ccd.header['DATE-OBS'], scale='utc')
    night.format = 'mjd'
    phot['night'] = np.int(np.floor(night.value - 0.5))
    phot['aperture_net_flux'] = phot['aperture_sum'] - phot['aperture_area'] * phot['sky_per_pix_avg']
    
    if gain is not None:
        phot['mag_inst_{}'.format(ccd.header['filter'])] = -2.5 * np.log10(gain * phot['aperture_net_flux'].value / phot['exposure'].value)

    metadata_to_add = ['AIRMASS', 'FILTER']
    for meta in metadata_to_add:
        phot[meta.lower()] = [ccd.header[meta]] * len(phot)    
    if fname:
        phot['file'] = fname
    if star_ids is not None:
        phot['star_id'] = star_ids<jupyter_output><empty_output><jupyter_text>We are ready to do photometry on all of the images Get the images in this folder...<jupyter_code>ifc = ImageFileCollection('.')<jupyter_output><empty_output><jupyter_text>Loop over all images and do photometry on them. + Change the object name if you need to.+ If there are images you want to skip put them in a different folder<jupyter_code>phots = []
missing_stars = []
for a_ccd, fname in ifc.ccds(object='kelt-16b', return_fname=True):
    print('on image ', fname)
    # Convert RA/Dec to pixel coordinates for this image
    pix_coords = a_ccd.wcs.all_world2pix(star_locs[0], star_locs[1], 0)
    xs, ys = pix_coords
    
    # Remove anything that is too close to the edges/out of frame
    padding = 3 * aperture_rad
    out_of_bounds = (xs < padding) | (xs > (a_ccd.shape[0] - padding)) | (ys < padding) | (ys > (a_ccd.shape[1] - padding))
    in_bounds = ~out_of_bounds
    
    # Find centroids of each region around star that is in_bounds
    xs_in = xs[in_bounds]
    ys_in = ys[in_bounds]
    print('    ...finding centroids')
    xcen, ycen = centroid_sources(a_ccd.data, xs_in, ys_in, box_size=2 * aperture_rad + 1)

    # Calculate offset between centroid in this image and the positions based on 
    # input RA/Dec. Later we wil set the magnitude of those with large differences
    # to an invalid value (maybe).
    center_diff = np.sqrt((xs_in - xcen)**2 + (ys_in - ycen)**2)

    #Set up apertures and annuli based on the centroids in this image.
    aps = CircularAperture((xcen, ycen), r=aperture_rad)
    anuls = CircularAnnulus((xcen, ycen), inner_annulus, outer_annulus)

    # Set any clearly bad values to NaN
    a_ccd.data[a_ccd.data > max_adu] = np.nan
    print('    ...doing photometry')
    # Do the photometry...
    pho = aperture_photometry(a_ccd.data, (aps, anuls), mask=a_ccd.mask, method='center')

    # We may have some stars we did not do photometry for because those stars were out of bounds.
    # Add the ones we missed to the list of missing
    missed = star_ids[out_of_bounds]
    missing_stars.append(missed)

    # Add all the extra goodies to the table
    print('    ...adding extra columns')
    add_to_photometry_table(pho, a_ccd, anuls, aps, fname=fname, star_ids=star_ids[in_bounds],
                            gain=1.47)
    # And add the final table to the list of tables
    phots.append(pho)<jupyter_output>WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING:astropy:FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system 
the RADECSYS keyword is deprecated, use RADESYSa.<jupyter_text>Combine all of the individual photometry tables into one<jupyter_code>all_phot = vstack(phots)<jupyter_output><empty_output><jupyter_text>Eliminate any stars that are missing from one or more imagesThis makes life a little easier later...<jupyter_code>uniques = set()
for miss in missing_stars:
    uniques.update(set(miss))

actually_bad = sorted([u for u in uniques if u in all_phot['star_id']])
len(uniques), len(actually_bad)
all_phot.add_index('star_id')
bad_rows = all_phot.loc_indices[actually_bad]
all_phot.remove_indices('star_id')
all_phot.remove_rows(sorted(bad_rows))<jupyter_output><empty_output><jupyter_text>Write out the file!<jupyter_code>all_phot
all_phot.write('all_the_photometry.fits')
all_phot.write('all_the_photometry.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10731589147286821, 'fraction_numerical': 0.014211886304909561, 'mean_word_length': 4.02637987012987, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '385700', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4370, 'n_tokens_neox': 4099, 'n_words': 1107, 'path': 'photometry-on-directory.ipynb', 'repo_name': 'mwcraig/aavso2018', 'size': '12384'}