Skip to content

fix forced photometry bug that skips first isophote #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ConnorStoneAstro merged 2 commits into from
Dec 16, 2025
Merged

fix forced photometry bug that skips first isophote #25

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fd0f465
fix forced photometry bug that skips first isophote
ConnorStoneAstro Dec 16, 2025
ecf45b9
slight change to find header
ConnorStoneAstro Dec 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
163 changes: 63 additions & 100 deletions src/autoprof/pipeline_steps/Isophote_Extract.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@

__all__ = ("Isophote_Extract_Forced", "Isophote_Extract", "Isophote_Extract_Photutils")


def _Generate_Profile(IMG, results, R, parameters, options):

# Create image array with background and mask applied
Expand Down Expand Up @@ -86,9 +87,7 @@ def _Generate_Profile(IMG, results, R, parameters, options):
compare_interp = []
for i in range(len(R)):
if "ap_isoband_fixed" in options and options["ap_isoband_fixed"]:
isobandwidth = (
options["ap_isoband_width"] if "ap_isoband_width" in options else 0.5
)
isobandwidth = options["ap_isoband_width"] if "ap_isoband_width" in options else 0.5
else:
isobandwidth = R[i] * (
options["ap_isoband_width"] if "ap_isoband_width" in options else 0.025
Expand Down Expand Up @@ -126,12 +125,10 @@ def _Generate_Profile(IMG, results, R, parameters, options):
else 5
),
sigmaclip=options["ap_isoclip"] if "ap_isoclip" in options else False,
sclip_iterations=options["ap_isoclip_iterations"]
if "ap_isoclip_iterations" in options
else 10,
sclip_nsigma=options["ap_isoclip_nsigma"]
if "ap_isoclip_nsigma" in options
else 5,
sclip_iterations=(
options["ap_isoclip_iterations"] if "ap_isoclip_iterations" in options else 10
),
sclip_nsigma=options["ap_isoclip_nsigma"] if "ap_isoclip_nsigma" in options else 5,
)
else:
isisophoteband = True
Expand All @@ -144,32 +141,21 @@ def _Generate_Profile(IMG, results, R, parameters, options):
mask=mask,
more=True,
sigmaclip=options["ap_isoclip"] if "ap_isoclip" in options else False,
sclip_iterations=options["ap_isoclip_iterations"]
if "ap_isoclip_iterations" in options
else 10,
sclip_nsigma=options["ap_isoclip_nsigma"]
if "ap_isoclip_nsigma" in options
else 5,
sclip_iterations=(
options["ap_isoclip_iterations"] if "ap_isoclip_iterations" in options else 10
),
sclip_nsigma=options["ap_isoclip_nsigma"] if "ap_isoclip_nsigma" in options else 5,
)
isotot = np.sum(
_iso_between(dat, 0, R[i], parameters[i], results["center"], mask=mask)
)
isotot = np.sum(_iso_between(dat, 0, R[i], parameters[i], results["center"], mask=mask))
medflux = _average(
isovals[0],
options["ap_isoaverage_method"]
if "ap_isoaverage_method" in options
else "median",
options["ap_isoaverage_method"] if "ap_isoaverage_method" in options else "median",
)
scatflux = _scatter(
isovals[0],
options["ap_isoaverage_method"]
if "ap_isoaverage_method" in options
else "median",
options["ap_isoaverage_method"] if "ap_isoaverage_method" in options else "median",
)
if (
"ap_iso_measurecoefs" in options
and not options["ap_iso_measurecoefs"] is None
):
if "ap_iso_measurecoefs" in options and not options["ap_iso_measurecoefs"] is None:
if (
mask is None
and (not "ap_isoclip" in options or not options["ap_isoclip"])
Expand Down Expand Up @@ -203,9 +189,7 @@ def _Generate_Profile(IMG, results, R, parameters, options):
cogdirect.append(isotot)
else:
sb.append(
flux_to_sb(medflux, options["ap_pixscale"], zeropoint)
if medflux > 0
else 99.999
flux_to_sb(medflux, options["ap_pixscale"], zeropoint) if medflux > 0 else 99.999
)
sbE.append(
(2.5 * scatflux / (np.sqrt(len(isovals[0])) * medflux * np.log(10)))
Expand Down Expand Up @@ -327,14 +311,10 @@ def _Generate_Profile(IMG, results, R, parameters, options):
SBprof_data["ellip"] = list(parameters[p]["ellip"] for p in range(end_prof))
SBprof_data["ellip_e"] = list(parameters[p]["ellip err"] for p in range(end_prof))
SBprof_data["pa"] = list(parameters[p]["pa"] * 180 / np.pi for p in range(end_prof))
SBprof_data["pa_e"] = list(
parameters[p]["pa err"] * 180 / np.pi for p in range(end_prof)
)
SBprof_data["pa_e"] = list(parameters[p]["pa err"] * 180 / np.pi for p in range(end_prof))
SBprof_data["pixels"] = list(pixels)
SBprof_data["maskedpixels"] = list(maskedpixels)
SBprof_data[
"totflux_direct" if fluxunits == "intensity" else "totmag_direct"
] = list(cogdirect)
SBprof_data["totflux_direct" if fluxunits == "intensity" else "totmag_direct"] = list(cogdirect)

if "ap_iso_measurecoefs" in options and not options["ap_iso_measurecoefs"] is None:
whichcoefs = [0] + list(options["ap_iso_measurecoefs"])
Expand Down Expand Up @@ -363,9 +343,7 @@ def _Generate_Profile(IMG, results, R, parameters, options):
SBprof_data["C"] = list(p["C"] for p in parameters[:end_prof])

if "ap_doplot" in options and options["ap_doplot"]:
Plot_Phase_Profile(
np.array(SBprof_data["R"]), parameters[:end_prof], results, options
)
Plot_Phase_Profile(np.array(SBprof_data["R"]), parameters[:end_prof], results, options)
if fluxunits == "intensity":
Plot_I_Profile(
dat,
Expand Down Expand Up @@ -406,7 +384,7 @@ def Isophote_Extract_Forced(IMG, results, options):
ap_fluxunits : str, default "mag"
units for outputted photometry. Can either be "mag" for log
units, or "intensity" for linear units.

ap_isoband_start : float, default 2
The noise level at which to begin sampling a band of pixels to
compute SB instead of sampling a line of pixels near the
Expand Down Expand Up @@ -483,7 +461,7 @@ def Isophote_Extract_Forced(IMG, results, options):
Tuple with axes limits for the y-axis in the SB profile
diagnostic plot. Be careful when using intensity units
since this will change the ideal axis limits.

ap_plot_sbprof_xlim : tuple, default None
Tuple with axes limits for the x-axis in the SB profile
diagnostic plot.
Expand All @@ -492,7 +470,7 @@ def Isophote_Extract_Forced(IMG, results, options):
Float value by which to scale errorbars on the SB profile
this makes them more visible in cases where the statistical
errors are very small.

Notes
----------
:References:
Expand Down Expand Up @@ -522,13 +500,23 @@ def Isophote_Extract_Forced(IMG, results, options):
with open(options["ap_forcing_profile"], "r") as f:
raw = f.readlines()
for i, l in enumerate(raw):
if l[0] != "#":
readfrom = i
break
if len(l.strip()) == 0 or l[0] == "#":
continue
readfrom = i
break
header = list(h.strip() for h in raw[readfrom].split(","))
force = dict((h, []) for h in header)
for l in raw[readfrom + 2 :]:
for d, h in zip(l.split(","), header):
for l in raw[readfrom + 1 :]:
if len(l) > 0 and l[0] == "#":
continue # skip comments
D = list(l.split(","))
if len(D) != len(header):
continue # Skip missmatched rows with header
try:
float(D[0].strip())
except ValueError:
continue # Skip non-numeric rows
for d, h in zip(D, header):
force[h].append(float(d.strip()))

force["pa"] = PA_shift_convention(np.array(force["pa"]), deg=True) * np.pi / 180
Expand All @@ -538,11 +526,7 @@ def Isophote_Extract_Forced(IMG, results, options):
"ellip": force["ellip"][i],
"pa": (
force["pa"][i]
+ (
options["ap_forced_pa_shift"]
if "ap_forced_pa_shift" in options
else 0.0
)
+ (options["ap_forced_pa_shift"] if "ap_forced_pa_shift" in options else 0.0)
)
% np.pi,
}
Expand Down Expand Up @@ -584,7 +568,7 @@ def Isophote_Extract(IMG, results, options):
ap_fluxunits : str, default "mag"
units for outputted photometry. Can either be "mag" for log
units, or "intensity" for linear units.

ap_samplegeometricscale : float, default 0.1
growth scale for isophotes when sampling for the final output
profile. Used when sampling geometrically. By default, each
Expand Down Expand Up @@ -691,16 +675,16 @@ def Isophote_Extract(IMG, results, options):
Tuple with axes limits for the y-axis in the SB profile
diagnostic plot. Be careful when using intensity units
since this will change the ideal axis limits.

ap_plot_sbprof_xlim : tuple, default None
Tuple with axes limits for the x-axis in the SB profile
diagnostic plot.

ap_plot_sbprof_set_errscale : float, default None
Float value by which to scale errorbars on the SB profile
this makes them more visible in cases where the statistical
errors are very small.

Notes
----------
:References:
Expand Down Expand Up @@ -735,9 +719,11 @@ def Isophote_Extract(IMG, results, options):

# Radius values to evaluate isophotes
R = [
options["ap_sampleinitR"]
if "ap_sampleinitR" in options
else min(1.0, results["psf fwhm"] / 2)
(
options["ap_sampleinitR"]
if "ap_sampleinitR" in options
else min(1.0, results["psf fwhm"] / 2)
)
]
while (
(
Expand All @@ -746,10 +732,7 @@ def Isophote_Extract(IMG, results, options):
)
or (options["ap_extractfull"] if "ap_extractfull" in options else False)
) and R[-1] < max(IMG.shape) / np.sqrt(2):
if (
"ap_samplestyle" in options
and options["ap_samplestyle"] == "geometric-linear"
):
if "ap_samplestyle" in options and options["ap_samplestyle"] == "geometric-linear":
if len(R) > 1 and abs(R[-1] - R[-2]) >= (
options["ap_samplelinearscale"]
if "ap_samplelinearscale" in options
Expand Down Expand Up @@ -797,24 +780,16 @@ def Isophote_Extract(IMG, results, options):
)
)
R = np.array(R)
logging.info(
"%s: R complete in range [%.1f,%.1f]" % (options["ap_name"], R[0], R[-1])
)
logging.info("%s: R complete in range [%.1f,%.1f]" % (options["ap_name"], R[0], R[-1]))

# Interpolate profile values, when extrapolating just take last point
tmp_pa_s = UnivariateSpline(
results["fit R"], np.sin(2 * results["fit pa"]), ext=3, s=0
)(R)
tmp_pa_c = UnivariateSpline(
results["fit R"], np.cos(2 * results["fit pa"]), ext=3, s=0
)(R)
tmp_pa_s = UnivariateSpline(results["fit R"], np.sin(2 * results["fit pa"]), ext=3, s=0)(R)
tmp_pa_c = UnivariateSpline(results["fit R"], np.cos(2 * results["fit pa"]), ext=3, s=0)(R)
E = _x_to_eps(
UnivariateSpline(
results["fit R"], _inv_x_to_eps(results["fit ellip"]), ext=3, s=0
)(R)
UnivariateSpline(results["fit R"], _inv_x_to_eps(results["fit ellip"]), ext=3, s=0)(R)
)
# np.arctan(tmp_pa_s / tmp_pa_c) + (np.pi * (tmp_pa_c < 0))
PA = _x_to_pa(((np.arctan2(tmp_pa_s, tmp_pa_c)) % (2 * np.pi)) / 2)
PA = _x_to_pa(((np.arctan2(tmp_pa_s, tmp_pa_c)) % (2 * np.pi)) / 2)
parameters = list({"ellip": E[i], "pa": PA[i]} for i in range(len(R)))

if "fit Fmodes" in results:
Expand Down Expand Up @@ -857,16 +832,12 @@ def Isophote_Extract(IMG, results, options):
and (not results["fit pa_err"] is None)
):
parameters[i]["ellip err"] = np.clip(
UnivariateSpline(
results["fit R"], results["fit ellip_err"], ext=3, s=0
)(R[i]),
UnivariateSpline(results["fit R"], results["fit ellip_err"], ext=3, s=0)(R[i]),
a_min=1e-3,
a_max=None,
)
parameters[i]["pa err"] = np.clip(
UnivariateSpline(results["fit R"], results["fit pa_err"], ext=3, s=0)(
R[i]
),
UnivariateSpline(results["fit R"], results["fit pa_err"], ext=3, s=0)(R[i]),
a_min=1e-3,
a_max=None,
)
Expand Down Expand Up @@ -895,21 +866,21 @@ def Isophote_Extract_Photutils(IMG, results, options):
ap_fluxunits : str, default "mag"
units for outputted photometry. Can either be "mag" for log
units, or "intensity" for linear units.

ap_plot_sbprof_ylim : tuple, default None
Tuple with axes limits for the y-axis in the SB profile
diagnostic plot. Be careful when using intensity units
since this will change the ideal axis limits.

ap_plot_sbprof_xlim : tuple, default None
Tuple with axes limits for the x-axis in the SB profile
diagnostic plot.

ap_plot_sbprof_set_errscale : float, default None
Float value by which to scale errorbars on the SB profile
this makes them more visible in cases where the statistical
errors are very small.

Notes
----------
:References:
Expand Down Expand Up @@ -1026,9 +997,7 @@ def Isophote_Extract_Photutils(IMG, results, options):
res = {}
dat = IMG - results["background"]
if not "fit R" in results and not "fit photutils isolist" in results:
logging.info(
"%s: photutils fitting and extracting image data" % options["ap_name"]
)
logging.info("%s: photutils fitting and extracting image data" % options["ap_name"])
geo = EllipseGeometry(
x0=results["center"]["x"],
y0=results["center"]["y"],
Expand All @@ -1042,8 +1011,7 @@ def Isophote_Extract_Photutils(IMG, results, options):
res.update(
{
"fit photutils isolist": isolist,
"auxfile fitlimit": "fit limit semi-major axis: %.2f pix"
% isolist.sma[-1],
"auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % isolist.sma[-1],
}
)
elif not "fit photutils isolist" in results:
Expand All @@ -1069,8 +1037,7 @@ def Isophote_Extract_Photutils(IMG, results, options):
res.update(
{
"fit photutils isolist": isolist,
"auxfile fitlimit": "fit limit semi-major axis: %.2f pix"
% isolist.sma[-1],
"auxfile fitlimit": "fit limit semi-major axis: %.2f pix" % isolist.sma[-1],
}
)
else:
Expand All @@ -1093,9 +1060,7 @@ def Isophote_Extract_Photutils(IMG, results, options):
zeropoint,
)
)
SBprof_data["SB_e"].append(
2.5 * isolist.int_err[i] / (isolist.intens[i] * np.log(10))
)
SBprof_data["SB_e"].append(2.5 * isolist.int_err[i] / (isolist.intens[i] * np.log(10)))
SBprof_data["totmag"].append(flux_to_mag(isolist.tflux_e[i], zeropoint))
SBprof_data["totmag_e"].append(
2.5
Expand All @@ -1117,9 +1082,7 @@ def Isophote_Extract_Photutils(IMG, results, options):
for k in SBprof_data.keys():
if not np.isfinite(SBprof_data[k][-1]):
SBprof_data[k][-1] = 99.999
res.update(
{"prof header": params, "prof units": SBprof_units, "prof data": SBprof_data}
)
res.update({"prof header": params, "prof units": SBprof_units, "prof data": SBprof_data})

if "ap_doplot" in options and options["ap_doplot"]:
if fluxunits == "intensity":
Expand Down