#!/usr/bin/env python
# coding: utf-8

# # Photometry Module
# 
# **Lecturer:** Chris Copperwheat<br>
# **Jupyter Notebook Author:** Kishalay De & Cameron Hummels
# 
# This is a Jupyter notebook lesson taken from the GROWTH Summer School 2020.  For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2020.html
# 
# ## Objective
# Measure photometric fluxes from astronomical ultraviolet, optical, infrared image data.
# 
# ## Key steps
# - Calibrate images to derive relationship between counts and brightness on sky.
# - Use aperture photometry tools to calculate brightness of sources.
# - Use Point Spread Function photometry tools to calculate brightness of sources.
# 
# ## Required dependencies
# 
# See GROWTH school webpage for detailed instructions on how to install these modules and packages.  Nominally, you should be able to install the python modules with `pip install <module>`.  The external astromatic packages are easiest installed using package managers (e.g., `rpm`, `apt-get`).
# 
# ### Python modules
# * python 3
# * astropy
# * numpy
# * matplotlib
# * astroquery
# * photutils
# 
# ### External packages
# * SExtractor https://www.astromatic.net/software
# * PSFex https://www.astromatic.net/software
# 
# Let's start by importing a few necessary modules first.

# In[ ]:


import numpy as np
import numpy.ma as ma
import os
import astropy.units as u
from astropy.io import ascii
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.stats import sigma_clipped_stats, sigma_clip
import subprocess


# ## Test dependencies
# 
# In order for this jupyter notebook to function correctly, we must have some external software installed, as described above.  The following step assures that these are installed properly before getting to the rest of the content of this lesson.

# In[ ]:


def test_dependency(dep, alternate_name=None):
    """
    Test external dependency by trying to run it as a subprocess
    """
    try:
        subprocess.check_output(dep, stderr=subprocess.PIPE, shell=True)
        print("%s is installed properly as %s. OK" % (dep, dep))
        return 1
    except subprocess.CalledProcessError:
        try:
            subprocess.check_output(alternate_name, stderr=subprocess.PIPE, shell=True)
            print("%s is installed properly as %s. OK" % (dep, alternate_name))
            return 1
        except subprocess.CalledProcessError:
            print("===%s/%s IS NOT YET INSTALLED PROPERLY===" % (dep, alternate_name))
            return 0
    
dependencies = [('sextractor', 'sex'), ('psfex', 'PSFEx')]
i = 0
for dep_name1, dep_name2 in dependencies:
    i += test_dependency(dep_name1, dep_name2)
print("%i out of %i external dependencies installed properly.\n" % (i, len(dependencies)))
if i != len(dependencies):
    print("Please correctly install these programs before continuing by following the instructions in README.md.")
else:
    print("You are ready to continue.") 


# Let's plot a reduced image from the previous module.

# In[ ]:


from astropy.io import fits
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import os

# Move to the data directory for our analysis
os.chdir('data')
imageName = 'aC0_20181013-174714-557.wcs.fits.proc.cr.fits'
f = fits.open(imageName)
data = f[0].data  #This is the image array
header = f[0].header

#Compute some image statistics for scaling the image plot
mean, median, sigma = sigma_clipped_stats(data)

#plot the image with some reasonable scale
plt.figure(figsize=(10,10))
plt.imshow(data, vmin=median-3*sigma, vmax=median+3*sigma)
plt.show()


# ## Querying an external catalog
# 
# The first step is to obtain the known magnitudes of sources in the field in order to subsequently compare them to the measured instrumental magnitudes. We thus have to refer to some external catalog for sources in the field.  For this tutorial, we will refer the Panstarrs 1 (PS1) catalog. We will query PS1 using astropy's query interface _astroquery_ to query the Vizier server and get a list of sources. First, we have to know the coordinates of the center of this field from the image header.

# In[ ]:


#strong the image WCS into an object
w = WCS(header)

#Get the RA and Dec of the center of the image
[raImage, decImage] = w.all_pix2world(data.shape[0]/2, data.shape[1]/2, 1)

#Set the box size to search for catalog stars
boxsize = 30 # arcminutes

#Magnitude cut-offs of sources to be cross-matched against
maxmag = 18


# Since PS1 is much deeper than a single image from the GROWTH India telescope, we also put a cut-off on the magnitude range of sources we want to calibrate against. For this specific example, a limiting magnitude of 18 mag in g band should suffice. We also want to reject transients and variable sources in PS1 (need a cut-off on the number of detections in the PS1 survey), and only select sources that have $>3 \sigma$ detections in PS1, so we will incorporate those filters into the PS1 query.

# In[ ]:


from astroquery.vizier import Vizier
#Vizier.VIZIER_SERVER = 'vizier.ast.cam.ac.uk'

catNum = 'II/349'#This is the catalog number of PS1 in Vizier
print('\nQuerying Vizier %s around RA %.4f, Dec %.4f with a radius of %.4f arcmin'%(catNum, raImage, decImage, boxsize))

try:
    #You can set the filters for the individual columns (magnitude range, number of detections) inside the Vizier query
    v = Vizier(columns=['*'], column_filters={"gmag":"<%.2f"%maxmag, "Nd":">6", "e_gmag":"<1.086/3"}, row_limit=-1)
    Q = v.query_region(SkyCoord(ra = raImage, dec = decImage, unit = (u.deg, u.deg)), radius = str(boxsize)+'m', catalog=catNum, cache=False)
    #query vizier around (ra, dec) with a radius of boxsize
    print(Q[0])
except:
    print('I cannnot reach the Vizier database. Is the internet working?')


# Just in case the Vizier query did not go through (this can happen, depending on the load on the Vizier servers), read in the source table in the file ps1_v641cyg.tab. No need to do this step if your query did go through.

# In[ ]:


from astropy.table import Table

Q = [Table.read('ps1_v641cyg.tab')]

print(Q[0])


# The query should have retrieved about 10,000 known sources in this field. We will now filter stars that are away from the edges of the image (same as in the case of the SExtractor detected sources). Note that you can also put additional filters in the catalog sources by just adding more conditions to the the commands below.
# 
# Using the same method as we used to figure out the center of the image, use the catalog coordinates from the astropy table to get image coordinates of the sources. Then create a table called good_cat_stars to only store stars that within the pixel range (500, 3500) in both x and y directions in the image.

# In[ ]:





# Let's get the image coordinates of the filtered good catalog stars and plot them on the image.

# In[ ]:


ps1_imCoords = w.all_world2pix(good_cat_stars['RAJ2000'],good_cat_stars['DEJ2000'], 1)

fig = plt.figure(figsize=(10,10))
ax = fig.gca()
plt.imshow(data, vmin=median-3*sigma, vmax=median+3*sigma)
circles = [plt.Circle((ps1_imCoords[0][i], ps1_imCoords[1][i]), radius = 5, edgecolor='r', facecolor='None') for i in range(len(ps1_imCoords[0]))]
for c in circles:
    ax.add_artist(c)
    
plt.show()


# ## Aperture photometry
# 
# To start with aperture photometry, we have to identify point sources in this image and compare their brightness to a known catalog. We will use a popular code known as 'Source Extractor' (a.k.a. SExtractor) to detect point sources in the image.

# ### Source detection
# 
# SExtractor is a very powerful code, and can do a large range of complex source detection tasks, but needs to be appropriately configured via its configuration and parameter file to specify the kind of sources we want to detect and the parameters we want to measure.
# 
# In this case, you have a configuration file available in the directory designed to detect sources above a 10$\sigma$ threshold from the background and measure it's flux. We'll go over the details, but let's see how SExtractor runs first.

# In[ ]:


configFile = 'photomCat.sex'
catalogName = imageName+'.cat'
paramName = 'photomCat.param'
try:
    command = 'sextractor -c %s %s -CATALOG_NAME %s -PARAMETERS_NAME %s' % (configFile, imageName, catalogName, paramName)
    print('Executing command: %s' % command)
    rval = subprocess.run(command.split(), check=True)
except subprocess.CalledProcessError as err:
    print('Could not run sextractor with exit error %s'%err)


# Here, we are running SExtractor with a configuration file specified with the '-c' flag and asking it to output the list of sources it deteced to a file called imageName+'.cat' in a format known as the FITS LDAC. The ouput parameters to be produced by SExtractor are in the file 'photomCat.param' and specified by the flag 'PARAMETERS_NAME'. We will now read in the source catalog produced by SExtractor into an astropy table. Before that, here is a function to read FITS LDAC tables with astropy.

# In[ ]:


def get_table_from_ldac(filename, frame=1):
    """
    Load an astropy table from a fits_ldac by frame (Since the ldac format has column 
    info for odd tables, giving it twce as many tables as a regular fits BinTableHDU,
    match the frame of a table to its corresponding frame in the ldac file).
    
    Parameters
    ----------
    filename: str
        Name of the file to open
    frame: int
        Number of the frame in a regular fits file
    """
    from astropy.table import Table
    if frame>0:
        frame = frame*2
    tbl = Table.read(filename, hdu=frame)
    return tbl


# In[ ]:


#This is a python wrapper for reading LDAC files produced by SExtractor
sourceTable = get_table_from_ldac(catalogName)
#Let's look at the contents of the table
print(sourceTable.colnames)
print(sourceTable)


# Note the columns in the output table -- important ones are 'XWIN_IMAGE' and 'YWIN_IMAGE' that are the centroids of the detected stars in image coordinates while 'ALPHAWIN_J2000' and 'DELTAWIN_J2000' are the same in world coordinates. We also have the background subtracted flux of the sources measured in fixed apertures ('FLUX_APER') and the magnitudes in fixed apertures ('MAG_APER').

# In order to derive a good photometric solution, we want to find 'clean' sources that are not blended with other objects and are away from the edges of the image, and get their instrumental magnitudes. The SExtractor output has a 'FLAGS' column that indicates unblended sources away from bad pixels with FLAGS = 0 (read the SExtractor manual to find out what other values of the flags mean). In order to get point-like sources, we also want the FWHM of the source to be smaller than about 2 arc seconds that are away from the edges of the image. So let's select the sources that satisfy this condition.

# In[ ]:


#filter on the sources to select the ones satisfying our criteria
cleanSources = sourceTable[(sourceTable['FLAGS']==0) & (sourceTable['FWHM_WORLD'] < 2) & (sourceTable['XWIN_IMAGE']<3500) & (sourceTable['XWIN_IMAGE']>500) &(sourceTable['YWIN_IMAGE']<3500) &(sourceTable['YWIN_IMAGE']>500)]
print(cleanSources)


# Now overlay these detected sources on the reduced image to see how sextractor performed

# In[ ]:





# ### Catalog cross-matching
# 
# So far, we have the catalog of sources detected in the image, and an external catalog of sources we want to compare against. The next step is to cross-match the two catalogs and find associations between sources detected in the image and sources on the sky.

# In[ ]:


sourceCatCoords = SkyCoord(ra=cleanSources['ALPHAWIN_J2000'], dec=cleanSources['DELTAWIN_J2000'], frame='icrs', unit='degree')
ps1CatCoords = SkyCoord(ra=good_cat_stars['RAJ2000'], dec=good_cat_stars['DEJ2000'], frame='icrs', unit='degree')

#Now cross match sources

#Set the cross-match distance threshold to 0.6 arcsec, or just about one pixel
photoDistThresh = 0.6
idx_image, idx_ps1, d2d, d3d = ps1CatCoords.search_around_sky(sourceCatCoords, photoDistThresh*u.arcsec)
#idx_image are indexes into sourceCatCoords for the matched sources, while idx_ps1 are indexes into ps1CatCoords for the matched sources

print('Found %d good cross-matches'%len(idx_image))


# ### Zero-point derivation
# 
# Now that we have cross-matched sources between the image and PS1, we can start deriving zero-points. Recall that the instrumental magnitude of a source is related to its physical magnitude via the zero-point relation
# 
# <h2> <center>
# 
# $m_{app}$ = ZP + $m_{ins}$ <br><br>
# 
# $m_{ins} = -2.5 \, \textrm{log} \, (\textrm{Total flux})$
# 
# </center></h2>
# 
# It should be easy to see that the instrumental magniutde of a source depends on the way you measure the flux, and this is exactly where the different ways of photometry come into play.
# 
# As we've already discussed, aperture photometry involves measuring the flux of a source in a fixed aperture of some radius. The amount of flux you measure is also dependent on the size of the aperture you choose. The larger the aperture, the more flux you capture from the wings of the PSF. At the same time, a larger aperture might also have contamination from neighboring sources, so you have to be careful.
# 
# In order to be consistent with flux measurements of your catalog calibrators and your target, you have to measure fluxes in the same apertures so that you are 'missing' the same fraction of flux between your calibrators and target source. As the zero-point depends on the difference between the catalog and instrumental magnitudes, the zero-point you derive will also be a function of the size of the aperture you choose -- so an aperture specific zero-point will allow us to account for these issues. Let's see an example of this.

# The SExtractor catalogs we got from the image actually also contains aperture photometry for the detected sources in the 'MAG_APER' keyword, so most of our work is actually already done. In this case, the SExtractor parameter file already specified the diameter of the apertures we wanted to measure: 4.0 pixels to 13.0 pixels in steps of 1.0 pixel. 
# 
# __Exercise__: Go ahead and look at the SExtractor configuration file and find where these apertures are specified.
# 
# To see how the measured instrumental magnitude changes between apertures of different sizes, we can plot the differences between the magnitude measured in the largest aperture (13.0 pixels) and every other aperture and plot the differences. 

# In[ ]:


#Apertures in the SExtractor configuration file
aperture_diameter = np.arange(4, 14)

#For each aperture, we are going to compute the magniutde difference between the largest pixel aperture and that specific aperture for every source in cross-matched catalog
magDiff = np.ma.zeros((len(aperture_diameter), len(idx_image)))
for j in range(len(aperture_diameter)):
    magDiff[j] = sigma_clip(cleanSources['MAG_APER'][:,9][idx_image] - cleanSources['MAG_APER'][:,j][idx_image])
#Here, magDiff is a 2D array contaning the difference magnitudes for each source and aperture
    
#Now, let's plot the magnitude differences 
plt.figure(figsize = (8,8))
plt.plot(aperture_diameter, magDiff, 'r.', markersize=1)
plt.xlabel('Aperture diameter (pixels)', fontsize=15)
plt.ylabel('Largest aperture mag - aperture mag', fontsize = 15)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.show()


# As apparent, the magnitudes become systematically fainter for smaller apertures. Recall that the larger the aperture, the more flux and the lower magnitude. This curve is often called as the __Curve of Growth__.
# 
# Now, let's overplot the instrumental (for any given aperture) and PS1 magnitudes for the good cross-matched sources.

# In[ ]:


plt.figure(figsize=(8,8))
#Plotting instrumental magnitude for aperture sizes of 5.0, 6.0 and 7.0 pixels
plt.plot(good_cat_stars['gmag'][idx_ps1], cleanSources['MAG_APER'][:,2][idx_image] , 'r.', label='5 pixel')
plt.plot(good_cat_stars['gmag'][idx_ps1], cleanSources['MAG_APER'][:,3][idx_image] , 'g.', label='6 pixel')
plt.plot(good_cat_stars['gmag'][idx_ps1], cleanSources['MAG_APER'][:,4][idx_image] , 'b.', label='7 pixel')
plt.ylim(-16, -7.5)
plt.xlabel('PS1 magnitude', fontsize=15)
plt.ylabel('Instrumental magnitude', fontsize=15)
plt.legend(fontsize=15)
plt.show()


# That is a good correlation between the PS1 and instrumental magnitudes, also showing the small systematic offsets between the magnitudes of different aperture sizes. 
# 
# We are ready to derive a zero-point for each aperture. From the earlier equation, the zero-point is just the difference between the catalog and instrumental magnitudes, so we want to measure the average offset between these magnitudes, including rejection of bad outliers that might affect the photometry. The results will be stored in an array of dictionaries containing the mean, median and standard deviation of the zero-point.

# In[ ]:


zeroPoints = []
for i in range(len(aperture_diameter)):
    #Array of differences between the catalog and instrumental magnitudes
    offsets = ma.array(good_cat_stars['gmag'][idx_ps1] - cleanSources['MAG_APER'][:,i][idx_image])
    #Compute sigma clipped statistics
    zero_mean, zero_med, zero_std = sigma_clipped_stats(offsets)
    zeroDict = {'diameter': aperture_diameter[i], 'zp_mean': zero_mean, 'zp_median': zero_med, 'zp_std': zero_std}
    zeroPoints.append(zeroDict)
    print(zeroDict)


# As apparent, the zero-point changes with the size of the aperture. So, when you want to derive photometry, you should measure fluxes in an aperture and apply the zero-point correction for the same aperture to get accurate photometry. This correction for the aperture size is known as an __Aperture correction__. 

# We will apply the respective zero-point corrections and compute uncertainties by adding the measured magnitude uncertainties and zero-point uncertainties in quadrature. We will use an astropy based package for photometry called photutils. photutils contains a large number of useful tools for performing photometry on images

# In[ ]:


from astropy import units as u
from photutils import SkyCircularAperture, SkyCircularAnnulus, aperture_photometry

ra = 324.7915750
dec = 46.7254686
#Specify the position of the target in sky coordinates
position = SkyCoord(ra = ra, dec = dec, unit = u.deg, frame = 'icrs')
#Now set up a set of apertures to measure flux from
aperture_radii = np.array([4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0]) / 2 #Aperutre radii in pixels 
apertures = [SkyCircularAperture(position, r = r * u.pix) for r in aperture_radii]
#Translate the apertures from sky coordinates to pixel coordinates
pix_apertures = [a.to_pixel(w) for a in apertures]
#Now compute the fluxes within these apertures on our data
phot_table = aperture_photometry(data, pix_apertures)
#Translating the output into usable format
for col in phot_table.colnames:
        phot_table[col].info.format = '%.8g'    
#print the output
print(phot_table)


# So far we've measured the flux inside a set of fixed size circular apertures on the image. We now need to compute the background around the source, so we will measure the flux in an annulus around the source coordinates. We will measure the background in annulus of radius 10.0 pixels and width of 10.0 pixels.

# In[ ]:


#Create the annulus aperture
anuRadius = 10.0
anuWidth = 10.0
annulus_aperture = SkyCircularAnnulus(position, r_in = anuRadius * u.pix, r_out = (anuRadius + anuWidth) * u.pix)
pix_annulus_aperture = annulus_aperture.to_pixel(w)
#Measuring the flux inside an aperture annulus
annulus_phot_table = aperture_photometry(data, pix_annulus_aperture)
for col in annulus_phot_table.colnames:
        annulus_phot_table[col].info.format = '%.8g'    
#print the output
print(annulus_phot_table)


# Now we are ready to measure the flux of the source by subtracting the background flux from the flux measured on the source.

# In[ ]:


#Estimate the background flux per unit pixel area
bkg_mean = annulus_phot_table['aperture_sum'] / pix_annulus_aperture.area
for i in range(len(pix_apertures)):
    #estimate the background flux for this specific aperture size
    bkg_flux = bkg_mean * pix_apertures[i].area
    #subtract the background flux for this aperture from the source flux
    source_flux = phot_table['aperture_sum_%d'%i] - bkg_flux
    source_mag = zeroPoints[i]['zp_median'] - 2.5 * np.log10(source_flux)
    print('Found source magnitude of %.2f for aperture of radius %d pixels'%(source_mag, aperture_radii[i]))


# __Question__: Think of a way to estimate the statistical error of this measurement assuming noise in the pixels is Poisson. Also recall that images are usually recorded in units of ADU (digital units), while the Poisson noise is applicable only to electron counts. How do you convert ADUs to electrons?

# Note that the aperture photometry magnitudes are completely consistent with each other within the errors of the measurements (the physical magnitude of a source should not depend on the aperture you used to measure the flux), suggesting that the aperture corrections worked fine. Also note that the errors become progressively larger for large apertures since we are accumulating more noise from the background. 
# 
# __Question__ : Can you think of another way we could have obtained by the aperture photomtery magnitude for this source? Think of what was contained in the SExtractor catalogs -- how would you extract the magnitudes for the source from the SExtractor catalog?

# ## PSF Photometry
# 
# Aperture photometry is fairly straightforward and works well for sources that are relatively isolated (there are no other sources in the aperture or the background). However, this may not be the case for sources in crowded fields (e.g. the galactic plane), where finding a source-free aperture can be difficult. This is where PSF photometry comes in.

# ### Creating a PSF model
# 
# The idea behind PSF photometry is to create a model for the Point Spread Function of stars in the field by fitting a large number of stars. The models are then iteratively subtracted from the data and refined to minimize the residuals. There are a number of PSF photometry packages in use today, with varying degrees of manual innvolvement and automation. For this module, we will use a package known as PSFEx, also from the Astromatic set of codes. 
# 
# PSFEx takes a SExtractor catalog as input. The catalog should be in the FITS LDAC format and contain cut-outs of the detected stars in the field, on which it performs the model fitting. For this module, we actually already produced the star cut-outs when we ran SExtractor the first time -- the 'PARAMETERS_NAME' file had an entry called VIGNET, specifying the size of the cut-out around each star. We will now run PSFEx on this catalog and try to produce a PSF model for the image. We again have an input configuration file called 'psfex_conf.psfex' that specified the run-time parameters for PSFEx.

# In[ ]:


psfConfigFile = 'psfex_conf.psfex'
try:
    command = 'psfex -c %s %s' % (psfConfigFile, catalogName)
    print('Executing command: %s' % command)
    rval = subprocess.run(command.split(), check=True)
except subprocess.CalledProcessError as err:
    print('Could not run psfex with exit error %s'%err)


# Now plot the PSF model produced by PSFEx. The best-fit Moffat model is stored as a FITS file in the output image 'moffat_' + imageName + '.fits'.

# In[ ]:





# It is also instructive to have a look at the radial profile of the best-fit PSF model. Let's make a radial profile and plot it.

# In[ ]:


psfImageCenter = [(psfModelData.shape[0]-1)/2, (psfModelData.shape[1]-1)/2]
y, x = np.indices(psfModelData.shape)
r = np.sqrt((x-psfImageCenter[0])**2 + (y-psfImageCenter[1])**2)
r = r.astype(np.int)

tbin = np.bincount(r.ravel(), psfModelData.ravel())
nr = np.bincount(r.ravel())
radialprofile = tbin/nr

plt.figure(figsize=(6,6))
plt.plot(range(len(radialprofile)), radialprofile, 'k.', markersize=15)
plt.xlabel('Radial distance (pixels)', fontsize=15)
plt.ylabel('PSF Amplitude', fontsize=15)
plt.xlim(-1,20)
plt.show()


# __Question__: Figure out why the above code is able to produce the radial profile of the PSF.

# ### Computing PSF-fit photometry for point sources
# 
# We can now feed this PSF model back to SExtractor to make it perform PSF model fitting on each of the detected sources and produce PSF-fit fluxes for each source. The full PSF model is stored in a file called imageName + '.psf'. This can take a bit of time depending on how fast your machine is.

# In[ ]:


psfName = imageName + '.psf'
psfcatalogName = imageName+'.psf.cat'
psfparamName = 'photomPSF.param' #This is a new set of parameters to be obtained from SExtractor, including PSF-fit magnitudes
try:
    #We are supplying SExtactor with the PSF model with the PSF_NAME option
    command = 'sextractor -c %s %s -CATALOG_NAME %s -PSF_NAME %s -PARAMETERS_NAME %s' % (configFile, imageName, psfcatalogName, psfName, psfparamName)
    print("Executing command: %s" % command)
    rval = subprocess.run(command.split(), check=True)
except subprocess.CalledProcessError as err:
    print('Could not run sextractor with exit error %s'%err)


# Let's read in the new PSF-fit SExtractor catalog as before and check it's contents. You should find additional new columns like 'MAG_POINTSOURCE' and 'MAGERR_POINTSOURCE' -- these are the new PSF model fit magnitudes measured with SExtractor.

# In[ ]:


psfsourceTable = get_table_from_ldac(psfcatalogName)
#Let's look at the contents of the table
print(psfsourceTable.colnames)
print(psfsourceTable)


# We are now ready to repeat effectively the same exercise as before -- compute an image zero-point based on PSF-fitted magnitudes. 'XMODEL_IMAGE' and 'YMODEL_IMAGE' are the centroids of the stars based on the PSF fitting. The 'FLAGS_MODEL' column can be used to find sources where the PSF-fitting photometry did not work well (the flag is set to 0 for good sources). We will repeat the cross-match between the good PSF-fit sources and the PS1 catalog.

# In[ ]:


#Selecting the clean sources away from image edges as before 
cleanPSFSources = psfsourceTable[(psfsourceTable['FLAGS']==0) & (psfsourceTable['FLAGS_MODEL']==0)  & (psfsourceTable['FWHM_WORLD'] < 2) & (psfsourceTable['XMODEL_IMAGE']<3500) & (psfsourceTable['XMODEL_IMAGE']>500) &(psfsourceTable['YMODEL_IMAGE']<3500) &(psfsourceTable['YMODEL_IMAGE']>500)]


# ### Catalog cross-matching
# 
# As we did before, perform a cross-match between the clean PSF-fit sources and the PS1 catalog

# In[ ]:





# Now, make a scatter plot of the PSF-fit magnitudes against the PS1 magnitudes for the matched sources

# In[ ]:





# ### Zero-point derivation
# 
# Note how the PSF-fit magnitudes deviate from a straight line at the bright end -- this is because the PSF starts becoming saturated on top of very bright stars, leading to poor PSF fits. These will be automatically clipped when we do the 3 sigma clipping for deriving the zero-point. As in the case of aperture photometry, compute a PSF photometry zero-point by sigma-clipping on the offsets between the PSF-fit and PS1 magnitudes. Print the results.

# In[ ]:


psfoffsets = ma.array(good_cat_stars['gmag'][idx_psfps1] - cleanPSFSources['MAG_POINTSOURCE'][idx_psfimage])
#Compute sigma clipped statistics
zero_psfmean, zero_psfmed, zero_psfstd = sigma_clipped_stats(psfoffsets)
print('PSF Mean ZP: %.2f\nPSF Median ZP: %.2f\nPSF STD ZP: %.2f'%(zero_psfmean, zero_psfmed, zero_psfstd))


# In order to obtain a PSF-fit magnitude for our favorite source, we can directly read off the SExtractor PSF-fit photometry table for the coordinates of our source by cross-matching the table coordinates to the coordinates of our source.

# In[ ]:


ra = 324.7915750
dec = 46.7254686

v641cyg_coords = SkyCoord(ra=[ra], dec=[dec], frame='icrs', unit='degree')
idx_v641cyg, idx_cleanpsf_v641cyg, d2d, d3d = psfsourceCatCoords.search_around_sky(v641cyg_coords, photoDistThresh*u.arcsec)
print('Found the source at index %d'%idx_cleanpsf_v641cyg[0])


# Now obtain the PSF-fit magnitude of the source from SExtractor catalog and print the result.

# In[ ]:





# The aperture-corrected magnitudes and PSF-fit magnitudes should be consistent with each other.