photometry_analysis_Python_code/photometry_analysis_Python_code.py at main · BoazRonZohar/photometry_analysis_Python_code · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 28 18:49:04 2025

@author: Lenovo
"""

import os
import astropy
import numpy as np
import pandas as pd
from astropy.io import fits
import matplotlib.pyplot as plt
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry
from astropy.stats import sigma_clipped_stats
from scipy.ndimage import gaussian_filter

# Function to calculate FWHM for each source
def compute_fwhm(data, x, y, size=10):
    """Measure FWHM around a light source."""
    # Check if the coordinates are within the image bounds
    x_min, x_max = int(x-size), int(x+size)
    y_min, y_max = int(y-size), int(y+size)

    # Ensure sub-image is within the bounds of the data array
    if x_min < 0 or y_min < 0 or x_max >= data.shape[1] or y_max >= data.shape[0]:
        print(f"Skipping source at ({x}, {y}) due to out-of-bounds sub-image.")
        return None

    # Extract sub-image centered around the source
    sub_image = data[y_min:y_max, x_min:x_max]
    smoothed = gaussian_filter(sub_image, sigma=2)  # Smoothing to reduce noise
    peak = np.max(smoothed)
    half_max = peak / 2

    # Find the width at half-max
    above_half_max = smoothed > half_max
    indices = np.argwhere(above_half_max)
    if indices.size > 0:
        min_x, max_x = indices[:, 1].min(), indices[:, 1].max()
        min_y, max_y = indices[:, 0].min(), indices[:, 0].max()
        fwhm_x = max_x - min_x
        fwhm_y = max_y - min_y
        return np.mean([fwhm_x, fwhm_y])
    return None

# Function to read a FITS file and perform source detection and measurements
def process_fits(filename, band):
    """Load a FITS image and detect and measure light sources."""
    hdul = fits.open(filename)
    data = hdul[0].data
    hdul.close()

    # Calculate background and standard deviation
    mean, median, std = sigma_clipped_stats(data, sigma=3.0)
    threshold = 5.0 * std  # Threshold for source detection

    # Source detection using DAOStarFinder
    daofind = DAOStarFinder(fwhm=3.0, threshold=threshold)
    sources = daofind(data - median)

    # Noise filtering: Keep only sources with sufficient peak amplitude
    sources = sources[sources['peak'] > 10 * std]

    # List to store results
    results = []
    for source in sources:
        x, y = source['xcentroid'], source['ycentroid']
        fwhm = compute_fwhm(data, x, y)
        if fwhm is not None:
            radius = 1.5 * fwhm  # Aperture radius based on FWHM
            aperture = CircularAperture((x, y), r=radius)
            annulus_inner_radius = radius * 1.5
            annulus_outer_radius = radius * 2.5  # Annulus for background measurement
            annulus = CircularAnnulus((x, y), r_in=annulus_inner_radius, r_out=annulus_outer_radius)

            # Perform photometry
            phot_table = aperture_photometry(data, [aperture, annulus])

            # Calculate background-subtracted flux
            background_mean = phot_table['aperture_sum_1'][0] / annulus.area  # Background mean
            background_subtracted_flux = phot_table['aperture_sum_0'][0] - background_mean * aperture.area

            # Apply a condition to discard unrealistic flux values
            if background_subtracted_flux < 0:  # Discard sources with negative flux values
                continue

            # Append the result to the list
            results.append([x, y, fwhm, radius, background_subtracted_flux, band, annulus_inner_radius, annulus_outer_radius])

    return results

# Function to generate unique CSV filename
def get_unique_filename(base_filename):
    index = 1
    # Check if the file already exists
    while os.path.exists(f"{base_filename}_{index}.csv"):
        index += 1
    return f"{base_filename}_{index}.csv"

# File path for the "B" band image (integrated from your specified path)
file_path = r"C:\Users\Lenovo\Desktop\M101 LCO  2024\1 Meter\BVR 1 Meter 23.5.2023\B\elp1m006-fa07-20230523-0074-e91.fits"

# Process the FITS image and obtain results
all_results = process_fits(file_path, "B")

# Convert results to a pandas DataFrame
df = pd.DataFrame(all_results, columns=["X", "Y", "FWHM", "Aperture Radius", "Flux", "Band", "Annulus Inner Radius", "Annulus Outer Radius"])

# Generate a unique filename for the CSV file
csv_filename = get_unique_filename("photometry_results_B")

# Save the DataFrame to the CSV file with the generated unique filename
df.to_csv(csv_filename, index=False)

print(f"Data saved to {csv_filename}")