Initial commit for aperture photometry

microlensing_photometry/__init__.py

@@ -0,0 +1 @@
+from . import logistics

microlensing_photometry/astrometry/__init__.py

@@ -0,0 +1 @@
+from . import wcs

microlensing_photometry/astrometry/wcs.py

@@ -0,0 +1,93 @@
+from skimage.registration import phase_cross_correlation
+import numpy as np
+from astropy.coordinates import SkyCoord
+import scipy.spatial as sspa
+from skimage.measure import ransac
+from skimage import transform as tf
+from astropy.wcs import WCS,utils
+
+import microlensing_photometry.logistics.GaiaTools.GaiaCatalog as GC
+
+def find_images_shifts(reference,image,image_fraction =0.25, upsample_factor=1):
+    """
+    Estimate the shifts between two images. Generally a good idea to do a fraction in the middle,
+    where the effect of rotation is minimal.
+
+    Parameters
+    ----------
+    reference : an image acting as the reference
+    image :  the image we want to align
+    image_fraction : the fraction of image around the center we want to analyze
+    upsample_factor : the degree of upsampling, if one wants subpixel accuracy
+
+
+    Returns
+    -------
+    shiftx : the shift in pixels in the x direction
+    shifty : the shift in pixels in the y direction
+    """
+
+    leny, lenx = (np.array(image.shape) * image_fraction).astype(int)
+    starty,startx = (np.array(image.shape)*0.5-[leny/2,lenx/2]).astype(int)
+
+    subref = reference.astype(float)[starty:starty+leny,startx:startx+lenx]
+    subimage = image.astype(float)[starty:starty+leny,startx:startx+lenx]
+    shifts, errors, phasediff = phase_cross_correlation(subref,subimage,
+                                                        normalization=None,
+                                                        upsample_factor=upsample_factor)
+
+    shifty,shiftx = shifts
+
+    return shiftx,shifty
+
+
+def refine_image_wcs(image, stars_image, image_wcs, gaia_catalog, star_limit = 1000):
+    """
+    Refine the WCS of an image with Gaia catalog
+
+    Parameters
+    ----------
+    image : an array with the image
+    stars_image :  the x,y positions of stars in the image
+    image_wcs : the original astropy WCS solution
+    gaia_catalog : the entire gaia catalog
+
+
+    Returns
+    -------
+    new_wcs : an updated astropy WCS object
+    """
+    mo_image,skycoords,star_pix = GC.build_Gaia_image(gaia_catalog, image_wcs, image_shape=image.shape)
+
+    shiftx, shifty = find_images_shifts(mo_image, image, image_fraction=0.25, upsample_factor=1)
+
+    dists = sspa.distance.cdist(stars_image[:star_limit],
+                                np.c_[star_pix[0][:star_limit] - shiftx, star_pix[1][:star_limit] - shifty])
+
+    mask = dists < 20
+    lines, cols = np.where(mask)
+
+    pts1 = np.c_[star_pix[0], star_pix[1]][:star_limit][cols]
+    pts2 = np.c_[stars_image[:,0], stars_image[:,1]][:star_limit][lines]
+    model_robust, inliers = ransac((pts2, pts1), tf.AffineTransform, min_samples=10, residual_threshold=3,
+                                   max_trials=300)
+
+    new_wcs = utils.fit_wcs_from_points(pts2[:star_limit][inliers].T, skycoords[cols][inliers])
+
+
+    ### might be valuable some looping here
+
+    # Refining???
+    # dists2 = distance.cdist(np.c_[sources['xcentroid'],sources['ycentroid']][:500],projected2[:500])
+    # mask2 = dists2<1
+    # lines2,cols2 = np.where(mask2)
+    # pts12 = np.c_[star_pix[0],star_pix[1]][:500][cols2]
+    # pts22 = np.c_[sources['xcentroid'],sources['ycentroid']][:500][lines2]
+
+    # model_robust2, inliers2 = ransac(( pts22,pts12), tf. AffineTransform,min_samples=10, residual_threshold=1, max_trials=300)
+    # projected22 = model_robust2.inverse(pts12)
+    # projected222 = model_robust2.inverse(np.c_[star_pix[0],star_pix[1]])
+
+    # print(shifts)
+
+    return new_wcs

microlensing_photometry/logistics/GaiaTools/GaiaCatalog.py

@@ -0,0 +1,105 @@
+from astroquery.gaia import Gaia
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+import numpy as np
+
+def collect_Gaia_catalog(ra,dec,radius=15,row_limit = 10000,catalog_name='Gaia_catalog.dat',
+                         catalog_path='./'):
+    """
+    Collect the Gaia catalog for the field ra,dec,radius.
+
+    Parameters
+    ----------
+    ra : right ascenscion in degree
+    dec : right ascenscion in degree
+    radius : radius in arcmin
+    row_limit : the maximum number of stars (default:10000)
+    catalog_name : the name of the catalog to save/load
+    catalog_path : the path where to find the catalog
+
+    Returns
+    -------
+    gaia_catalog : astropy.Table, sorted by magnitude (brighter to fainter)
+    """
+
+    try:
+
+        from astropy.table import Table
+
+        gaia_catalog = Table.read(catalog_name, format='ascii')
+
+    except:
+
+        Gaia.ROW_LIMIT = row_limit
+        coord = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree), frame='icrs')
+
+        radius = u.Quantity(radius /60., u.deg)
+
+        gaia_catalog = Gaia.query_object_async(coordinate=coord,radius = radius)
+        gaia_catalog =  gaia_catalog[ gaia_catalog['phot_g_mean_mag'].argsort(),]
+        gaia_catalog.write(catalog_name, format='ascii')
+
+
+    return  gaia_catalog
+
+def build_Gaia_image(gaia_catalog,image_wcs, image_shape=(4096,4096), star_limit = 1000):
+    """
+    Construct a model image of the Gaia catalog, useful for estimating intial X,Y shifts.
+
+    Parameters
+    ----------
+    gaia_catalog : an astropy table containing the Gaia catalog
+    image_wcs :  an astropy wcs object to convert ra,dec to x,y
+    image_shape : the shape of the image to build
+    star_limit : limit to the X brightest star to simulate
+
+    Returns
+    -------
+    gaia_image : a model image of the Gaia catalog
+    """
+
+    coords = SkyCoord(ra=gaia_catalog['ra'].data[:star_limit],
+                      dec=gaia_catalog['dec'].data[:star_limit],
+                      unit=(u.degree, u.degree), frame='icrs')
+
+    fluxes = 10**((22-gaia_catalog['phot_g_mean_mag'].data)/2.5)
+
+    star_pix = image_wcs.world_to_pixel(coords)
+    XX, YY = np.indices((13, 13))
+    model_gaussian = Gaussian2d(1, 6, 6, 2, 2, XX, YY)
+
+    model_image = np.zeros(image_shape)
+
+    for ind in range(len(star_pix[0][:star_limit])):
+        try:
+            model_image[star_pix[1].astype(int)[ind] - 6:star_pix[1].astype(int)[ind] + 7,
+            star_pix[0].astype(int)[ind] - 6:star_pix[0].astype(int)[ind] + 7] += model_gaussian*fluxes[ind]
+        except:
+            #skip borders
+            pass
+
+    return model_image,coords,star_pix
+
+def Gaussian2d(intensity,x_center,y_center,width_x,width_y,X_star,Y_star):
+    """
+    A 2d Gaussian model, to ingest stars in image model
+
+     Parameters
+     ----------
+     intensity : the intensity of the gaussian
+     x_center :  the gaussian center in x direction
+     y_center :  the gaussian center in y direction
+     width_x :  the gaussian sigma x
+     width_y :  the gaussian sigma y
+     X_star :  2D  array containing the X value
+     Y_star :  2D  array containing the Y value
+
+     Returns
+     -------
+     model : return the 2D gaussiam
+     """
+    model = intensity * np.exp(
+            -(((X_star -x_center) / width_x) ** 2 + \
+              ((Y_star - y_center) / width_y) ** 2) / 2)
+
+    return model

microlensing_photometry/logistics/GaiaTools/__init__.py

@@ -0,0 +1 @@
+from . import GaiaCatalog

microlensing_photometry/logistics/__init__.py

@@ -0,0 +1 @@
+from . import GaiaTools

microlensing_photometry/photometry/__init__.py

@@ -0,0 +1,2 @@
+from . import aperture_photometry
+from . import photometric_scale_factor

microlensing_photometry/photometry/aperture_photometry.py

@@ -0,0 +1,105 @@
+from photutils.aperture import aperture_photometry
+from photutils.aperture import CircularAnnulus, CircularAperture
+from photutils.aperture import ApertureStats
+from astropy.io import fits
+from astropy.wcs import WCS
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+import numpy as np
+
+from microlensing_photometry.astrometry import wcs as lcowcs
+
+class AperturePhotometryAnalyst(object):
+    """
+    Class of worker to perform WCS refinement and aperture photometry for one image
+
+    Attributes
+    ----------
+    image_path : a path+name to the data
+    gaia_catalog : the Gaia catalog of the field
+
+
+    Returns
+    -------
+    new_wcs : an updated astropy WCS object
+    """
+    def __init__(self,image_path,gaia_catalog):
+
+        self.image_path = image_path
+        self.image_layers = fits.open(self.image_path)
+        self.gaia_catalog = gaia_catalog
+
+        ### To fix with config files
+        self.image_data = self.image_layers[0].data
+        self.image_errors = self.image_layers[3].data
+        self.image_original_wcs = WCS(self.image_layers[0].header)
+
+        self.process_image()
+
+    def process_image(self):
+
+        self.find_star_catalog()
+        self.refine_wcs()
+        self.run_aperture_photometry()
+
+    def find_star_catalog(self):
+
+        if self.image_layers[1].header['EXTNAME']=='CAT':
+            #BANZAI image
+            self.star_catalog = np.c_[self.image_layers[1].data['x'],
+                                      self.image_layers[1].data['y']]
+        else:
+            ### run starfinder
+            pass
+
+    def refine_wcs(self):
+
+        wcs2 = lcowcs.refine_image_wcs(self.image_data , self.star_catalog,
+                                       self.image_original_wcs, self.gaia_catalog,
+                                       star_limit = 1000)
+
+        self.image_new_wcs = wcs2
+
+    def run_aperture_photometry(self):
+
+        skycoord = SkyCoord(ra=self.gaia_catalog['ra'], dec=self.gaia_catalog['dec'], unit=(u.degree, u.degree))
+        xx, yy = self.image_new_wcs.world_to_pixel(skycoord)
+        positions = np.c_[xx,yy]
+
+        phot_table = run_aperture_photometry(self.image_data, self.image_errors, positions, 3)
+        exptime = self.image_layers[0].header['EXPTIME']
+
+        phot_table['aperture_sum'] /= exptime
+        phot_table['aperture_sum_err'] /= exptime
+
+        self.aperture_photometry_table = phot_table
+
+def run_aperture_photometry(image, error, positions,radius):
+    """
+    Aperture photometry on a image
+
+    Parameters
+    ----------
+    image : the image data
+    error : the error data (2D)
+    positions: [X,Y] positions of the stars to place aperture
+    radius: aperture radius
+
+    Returns
+    -------
+    phot_table : the photometric table
+    """
+    aperture = CircularAperture(positions, r=radius)
+    annulus_aperture = CircularAnnulus(positions, r_in=radius+3, r_out=radius+5)
+
+
+    aperstats = ApertureStats(image, annulus_aperture)
+
+    bkg_mean = aperstats.mean
+
+    phot_table = aperture_photometry(image, aperture, error=error)
+    total_bkg = aperture.area * bkg_mean
+
+    phot_table['aperture_sum'] -= total_bkg
+
+    return phot_table

microlensing_photometry/photometry/photometric_scale_factor.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+
+def photometric_scale_factor_from_lightcurves(lcs):
+    """
+    Estimate the 16,50,84 percentiles of the photometric scale factor
+
+    Parameters
+    ----------
+    lcs : an array containing all the lightcurves
+
+    Returns
+    -------
+    pscales : the photometric scale factors
+    """
+    pscales = np.nanpercentile(lcs/np.nanmedian(lcs,axis=1)[:,None],[16,50,84],axis=0)
+
+    return pscales

pyproject.toml

@@ -14,3 +14,4 @@ fsspec = "2024.10.0"
 [build-system]
 requires = ["poetry-core"]
 build-backend = "poetry.core.masonry.api"
+