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,107 @@
+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 astropy.units as u
+from astropy.coordinates import SkyCoord
+
+from microlensing_photometry.logistics import image_tools
+
+def find_images_shifts(reference,image,image_fraction =0.25, upsample_factor=1):
+    """
+    Estimate the shifts (X,Y) between two images. Generally a good idea to do only a fraction of the field of view
+    centred in the middle, where the effect of rotation is minimal.
+
+    Parameters
+    ----------
+    reference : array,  an image acting as the reference
+    image : array,  the image we want to align
+    image_fraction : float, the fraction of image around the center we want to analyze
+    upsample_factor : float, the degree of upsampling, if one wants subpixel accuracy
+
+
+    Returns
+    -------
+    shiftx : float, the shift in pixels in the x direction
+    shifty : float, 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. First, find shifts in X,Y between the image stars catalog and
+    a model image of the Gaia catalog. Then compute the full WCS solution using ransac and a affine transform.
+
+    Parameters
+    ----------
+    image : array, the image to refine the WCS solution
+    stars_image : array, the x,y positions of stars in the image
+    image_wcs : astropy.wcs, the original astropy WCS solution
+    gaia_catalog : astropy.Table, the entire gaia catalog
+    star_limit : int, the limit number of stars to use
+
+    Returns
+    -------
+    new_wcs : astropy.wcs, an updated astropy WCS object
+    """
+
+    skycoords = SkyCoord(ra=gaia_catalog['ra'].data[:star_limit],
+                      dec=gaia_catalog['dec'].data[:star_limit],
+                      unit=(u.degree, u.degree), frame='icrs')
+
+    fluxes = [1]*len(gaia_catalog['phot_g_mean_flux'].data)
+    star_pix = image_wcs.world_to_pixel(skycoords)
+
+    stars_positions = np.array(star_pix).T
+
+    model_gaia_image = image_tools.build_image(stars_positions, fluxes, image.shape,
+                                          image_fraction=1)
+
+    model_image = image_tools.build_image(stars_image[:,:2], [1]*len(stars_image), image.shape, image_fraction=1)
+
+    shiftx, shifty = find_images_shifts(model_gaia_image, model_image, image_fraction=0.25, upsample_factor=1)
+
+    dists = sspa.distance.cdist(stars_image[:star_limit,:2],
+                                np.c_[star_pix[0][:star_limit] - shiftx, star_pix[1][:star_limit] - shifty])
+    mask = dists < 10
+    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=5,
+                                   max_trials=300)
+
+    new_wcs = utils.fit_wcs_from_points(pts2[:star_limit][inliers].T, skycoords[cols][inliers])
+
+    #breakpoint()
+    ### 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/infrastructure/__init__.py

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

microlensing_photometry/infrastructure/logs.py

@@ -0,0 +1,105 @@
+import logging
+from os import path, remove, makedirs
+from sys import exit
+from astropy.time import Time
+from datetime import datetime
+import glob
+
+
+def start_log(log_dir, log_name):
+    """
+    Function to initialize a log file for a single entity.  The new file will automatically overwrite any
+    previously-existing logfile with the same name
+
+    This function also configures the log file to provide timestamps for
+    all entries.
+
+    Parameters
+    ----------
+    log_dir : str, the path of the log file
+    log_name : str, the name of the log file
+
+    Returns
+    -------
+    log : logger, an open logger object
+    """
+
+    # Console output not captured, though code remains for testing purposes
+    console = False
+
+    if path.isdir(log_dir) == False:
+        makedirs(log_dir)
+
+    log_file = path.join(log_dir, log_name + '.log')
+
+    if path.isfile(log_file) == True:
+        remove(log_file)
+
+    # To capture the logging stream from the whole script, create
+    # a log instance together with a console handler.  
+    # Set formatting as appropriate.
+    log = logging.getLogger('analyst_'+log_name)
+
+    if len(log.handlers) == 0:
+        log.setLevel(logging.INFO)
+        file_handler = logging.FileHandler(log_file)
+        file_handler.setLevel(logging.INFO)
+
+        if console == True:
+            console_handler = logging.StreamHandler()
+            console_handler.setLevel(logging.INFO)
+
+        formatter = logging.Formatter(fmt='%(asctime)s %(message)s', \
+                                      datefmt='%Y-%m-%dT%H:%M:%S')
+        file_handler.setFormatter(formatter)
+
+        if console == True:
+            console_handler.setFormatter(formatter)
+
+        log.addHandler(file_handler)
+        if console == True:
+            log.addHandler(console_handler)
+
+    log.info('Process start with opening log for ' + log_name + '\n')
+
+    return log
+
+
+def ifverbose(log, setup, string):
+    """Function to write to a logfile only if the verbose parameter in the
+    metadata is set to True"""
+
+    ### Not sure here, I comment
+
+    #if log != None and setup.verbosity >= 1:
+    #    log.info(string)
+
+    #if setup.verbosity == 2:
+
+    #    try:
+
+    #        print(string)
+
+    #    except IOError:
+
+    #        pass
+
+
+def close_log(log):
+    '''
+    Function that closes a log.
+
+    Parameters
+    ----------
+    log: logger, a logger instance to close
+    '''
+
+    #log.info( 'Processing complete\n' )
+    #logging.shutdown()
+
+    for handler in log.handlers:
+        log.info('Processing complete.\n')
+        if isinstance(handler, logging.FileHandler):
+            handler.close()
+        log.removeFilter(handler)
+

microlensing_photometry/logistics/GaiaTools/GaiaCatalog.py

@@ -0,0 +1,44 @@
+from astroquery.gaia import Gaia
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+import numpy as np
+
+from microlensing_photometry.logistics import vizier_tools
+
+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 centered at ra,dec and 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_catalog = vizier_tools.search_vizier_for_sources(ra, dec, radius, 'Gaia-EDR3', row_limit=-1,
+                                  coords='degree', log=None, debug=False)
+        mask = np.isfinite(gaia_catalog['phot_g_mean_flux'])
+        sub_gaia_catalog = gaia_catalog[mask]
+        sub_gaia_catalog =  sub_gaia_catalog[sub_gaia_catalog['phot_g_mean_flux'].argsort()[::-1],]
+        sub_gaia_catalog.write(catalog_name, format='ascii')
+
+        gaia_catalog = sub_gaia_catalog
+
+    return  gaia_catalog

microlensing_photometry/logistics/GaiaTools/__init__.py

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

microlensing_photometry/logistics/__init__.py

@@ -0,0 +1,3 @@
+from . import GaiaTools
+from . import vizier_tools
+from . import image_tools

microlensing_photometry/logistics/image_tools.py

@@ -0,0 +1,46 @@
+import numpy as np
+
+from microlensing_photometry.photometry import psf
+
+def build_image(star_positions, fluxes, image_shape, image_fraction = 0.25, star_limit = 1000):
+    """
+    Construct an image with fake stars on top of a null background. Only the star inside the image fraction are computed
+
+    Parameters
+    ----------
+    image : array, the image to refine the WCS solution
+    stars_image : array, the x,y positions of stars in the image
+    image_wcs : astropy.wcs, the original astropy WCS solution
+    gaia_catalog : astropy.Table, the entire gaia catalog
+    star_limit : int, the limit number of stars to use
+
+    Returns
+    -------
+    model_image : array, a model image of the field
+    """
+
+    leny, lenx = (np.array(image_shape) * image_fraction).astype(int)
+    mask = ((np.abs(star_positions[:,1] - image_shape[0] / 2) < leny)
+            & (np.abs(star_positions[:,0] - image_shape[1] / 2) < lenx))
+
+    sub_star_positions = np.array(star_positions)[mask][:star_limit]
+
+    XX, YY = np.indices((21, 21))
+
+    model_gaussian = psf.Gaussian2d(1, XX.max()/2, XX.max()/2, 3, 3, XX, YY)
+
+    model_image = np.zeros(image_shape)
+
+    for ind in range(len(sub_star_positions)):
+
+        try:
+
+            model_image[sub_star_positions[ind,1].astype(int) - int(XX.max()/2):sub_star_positions[ind,1].astype(int) + int(XX.max()/2+1),
+            sub_star_positions[ind,0].astype(int) - int(XX.max()/2):sub_star_positions[ind,0].astype(int) + int(XX.max()/2+1)] += model_gaussian * fluxes[ind]
+
+        except:
+
+            # skip borders
+            pass
+
+    return model_image