#!/usr/bin/env python3
import os
import os.path
import subprocess as sp
import sys
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from galfitools.galin.galfit import Galfit
from galfitools.galin.std import GetAxis
from galfitools.galin.std import GetSize
from galfitools.galin.std import Ds9ell2Kronellv2
from galfitools.galin.std import GetInfoEllip
from galfitools.galin.std import GetPmax
from galfitools.galin.std import GetExpTime
from galfitools.galin.MakeMask import CheckFlag
from mgefit.mge_fit_sectors import mge_fit_sectors
from mgefit.mge_fit_sectors_regularized import mge_fit_sectors_regularized
from mgefit.mge_fit_sectors_twist import mge_fit_sectors_twist
from mgefit.sectors_photometry import sectors_photometry
from mgefit.sectors_photometry_twist import sectors_photometry_twist
from scipy.special import gammaincinv
[docs]
def mge2gal(
galfitFile,
regfile,
center,
psf,
twist,
gauss,
freeser,
freesky,
numgauss,
sky=0,
xypos=None,
ellip=None,
posang=None,
) -> str:
"""
Creates MGE initial parameters for GALFIT
Creates a Multi-Gaussian Expansion (MGE) model and
formats it into an initial parameter file for fitting by GALFIT.
galfitFile : str
GALFIT file from which the header information
will be extracted to create the MGE model.
regfile : str
DS9 ellipse region file, which must enclose the galaxy to be fitted.
center : Bool
if True it will take the geometric's center of the DS9 ellipse
as the center, otherwise it will take the pixel with the peak's value
psf : float
value of the PSF sigma
twist : bool
If True, the twist option for MGE is enabled, allowing the
angular positions of the Gaussians to differ from one another.
gauss : bool
if True, it uses the gaussian model instead of the Sersic model
with n = 0.5
freeser : bool
leaves the sersic index as a free parameter to fit
freesky : bool
leaves the sky parameter as a free parameter to fit
numgauss : int
maximum number of gaussians allowed to fit
sky : float
value of the sky background. Default = 0
xypos : list, optional
provides the (x y) position center of the object to fit
ellip : float
ellipticity of the object.
posang : float
position angle of object. Measured from Y-axis
Returns
-------
T2 : str
name of the output FITS
Notes
-----
The creation of initial parameters for MGE are created through
the routine of Cappellari, MNRAS 333,400-410 (2002)
"""
# reading options from galfit header
galfit = Galfit(galfitFile)
head = galfit.ReadHead()
# galsky = galfit.ReadSky()
image = head.inputimage
imageout = head.outimage
magzpt = head.mgzpt
maskfile = head.maskimage
scale = head.scale
psfile = head.psfimage
sigfile = head.sigimage
convbox = head.convx
convboxy = head.convy
consfile = head.constraints
xlo = head.xmin
ylo = head.ymin
xhi = head.xmax
yhi = head.ymax
#################
initgauss = 12
regu = None # regu option now available
mgeoutfile = "mgegas.txt"
# Updating variables for empty files ###########
if (maskfile == "None") or (maskfile == "none"):
maskfile = None
if (sigfile == "None") or (sigfile == "none"): # pragma: no cover
sigfile = None
if (psfile == "None") or (psfile == "none"): # pragma: no cover
psfile = None
if regu: # pragma: no cover
try:
from mgefit.mge_fit_sectors_twist_regularized import (
mge_fit_sectors_twist_regularized,
)
except Exception:
print(
"mge_fit_sectors_twist_regularized is not installed. Regu desactivated"
)
regu = False
##########################################
############################################
if regu: # pragma: no cover
print("regularization mode activated")
root_ext = os.path.splitext(image)
timg = root_ext[0]
namepng = timg + ".png"
sectorspng = "sectors.png"
magzpt = float(magzpt)
fit = 1
serfit = 0
skyfit = 0
Z = 0
###################################################
# Create GALFIT input file header to compute sky #
###################################################
# image = "ngc4342.fits"
root_ext = os.path.splitext(imageout)
TNAM = root_ext[0]
exptime = GetExpTime(image)
######################
# Mask file
#################
if maskfile:
errmsg = "file {} does not exist".format(maskfile)
assert os.path.isfile(maskfile), errmsg
hdu = fits.open(maskfile)
mask = hdu[0].data
maskb = np.array(mask, dtype=bool)
hdu.close()
else:
mask = np.array([])
######################
# sky=back
######################
hdu = fits.open(image)
img = hdu[0].data
img = img.astype(float)
minlevel = 0
# plt.clf()
(ncol, nrow) = GetAxis(image)
eps = 0
theta = 0
if regfile: # pragma: no cover
if os.path.isfile(regfile):
obj, xpos, ypos, rx, ry, angle = GetInfoEllip(regfile)
xx, yy, Rkron, theta, eps = Ds9ell2Kronellv2(xpos, ypos, rx, ry, angle)
if center:
print("center of ds9 ellipse region will be used")
xpeak, ypeak = xpos, ypos
else:
(xmin, xmax, ymin, ymax) = GetSize(
xx, yy, Rkron, theta + 90, eps, ncol, nrow
)
xpeak, ypeak = GetPmax(img, mask, xmin, xmax, ymin, ymax)
if xypos:
xpeak = xypos[0] - 1
ypeak = xypos[1] - 1
if ellip:
eps = ellip
if posang:
theta = posang
print("galaxy found at ", xpeak + 1, ypeak + 1)
print("Ellipticity, Angle = ", eps, theta)
print("Sky = ", sky)
img = img - sky # subtract sky
# I have to switch x and y coordinates, don't ask me why
xtemp = xpeak
xpeak = ypeak
ypeak = xtemp
# plt.clf()
if twist: # pragma: no cover
print("twist mode mge fit sectors is used")
# Perform galaxy photometry
if maskfile:
s = sectors_photometry_twist(
img, theta, xpeak, ypeak, minlevel=minlevel, badpixels=maskb, plot=False
)
else:
s = sectors_photometry_twist(
img, theta, xpeak, ypeak, minlevel=minlevel, plot=False
)
# plt.savefig(sectorspng)
# plt.pause(1) # Allow plot to appear on the screen
# plt.clf()
tolpsf = 0.001
if np.abs(psf) > tolpsf:
if regu:
m = mge_fit_sectors_twist_regularized(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
sigmapsf=psf,
scale=scale,
plot=False,
)
else:
m = mge_fit_sectors_twist(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
sigmapsf=psf,
scale=scale,
plot=False,
)
else:
print("No convolution")
if regu:
m = mge_fit_sectors_twist_regularized(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
scale=scale,
plot=False,
)
else:
m = mge_fit_sectors_twist(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
scale=scale,
plot=False,
)
# plt.pause(1) # Allow plot to appear on the screen
# plt.savefig(namepng)
else:
# elif twist == False:
if maskfile:
s = sectors_photometry(
img,
eps,
theta,
xpeak,
ypeak,
minlevel=minlevel,
badpixels=maskb,
plot=False,
)
else:
s = sectors_photometry(
img, eps, theta, xpeak, ypeak, minlevel=minlevel, plot=False
)
# plt.pause(1) # Allow plot to appear on the screen
# plt.savefig(sectorspng)
# plt.clf()
tolpsf = 0.001
if np.abs(psf) > tolpsf:
if regu: # pragma: no cover
m = mge_fit_sectors_regularized(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
sigmapsf=psf,
scale=scale,
plot=False,
bulge_disk=0,
linear=0,
)
else:
m = mge_fit_sectors(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
sigmapsf=psf,
scale=scale,
plot=False,
bulge_disk=0,
linear=0,
)
else:
print("No convolution")
if regu: # pragma: no cover
m = mge_fit_sectors_regularized(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
scale=scale,
plot=False,
bulge_disk=0,
linear=0,
)
else:
m = mge_fit_sectors(
s.radius,
s.angle,
s.counts,
eps,
ngauss=initgauss,
scale=scale,
plot=False,
bulge_disk=0,
linear=0,
)
# plt.pause(1) # Allow plot to appear on the screen
# plt.savefig(namepng)
if twist: # pragma: no cover
(counts, sigma, axisrat, pa) = m.sol
theta2 = 270 - theta
alpha1 = pa - theta2
alphaf = alpha1 - 90
else:
# elif twist == False:
(counts, sigma, axisrat) = m.sol
# anglegass = 90 - theta
anglegass = theta
####################
# print GALFIT files
#####################
if gauss: # pragma: no cover
parfile = "mgeGALFIT.txt"
else:
parfile = "mseGALFIT.txt"
outname = TNAM
rmsname = sigfile
if psfile: # pragma: no cover
psfname = psfile
else:
psfname = "None"
# switch back
xtemp = xpeak
xpeak = ypeak
ypeak = xtemp
T1 = "{}".format(image)
T2 = outname + "-mge.fits"
T3 = "{}".format(rmsname)
# now this is read from galfit file
# xlo=1
# ylo=1
# (xhi, yhi) = (ncol, nrow)
fout1 = open(parfile, "w")
fout2 = open(mgeoutfile, "w")
outline2 = "# Mag Sig(pixels) FWHM(pixels) q angle \n"
fout2.write(outline2)
# sersic K for n = 0.5
K = gammaincinv(1, 0.5)
if psfile: # pragma: no cover
(ncol, nrow) = GetAxis(psfname)
convbox = ncol + 1
convboxy = nrow + 1
PrintHeader(
fout1,
T1,
T2,
T3,
psfname,
1,
maskfile,
consfile,
xlo,
xhi,
ylo,
yhi,
convbox,
convboxy,
magzpt,
scale,
scale,
"regular",
0,
0,
)
if freeser: # pragma: no cover
serfit = 1
else:
serfit = 0
# removing the last components:
totGauss = len(counts)
if not (numgauss): # pragma: no cover
numgauss = totGauss
print("total number of gaussians by mge_fit_sectors: ", totGauss)
if numgauss > totGauss: # pragma: no cover
print("number of gaussians is greater than the ones fitted by mge_fit_sectors")
print("total number of gaussians of mge_fit_sectors will be used instead")
numgauss = totGauss
print("total number of gaussians used by GALFIT: ", numgauss)
index = 0
while index < numgauss:
TotCounts = counts[index]
SigPix = sigma[index]
qobs = axisrat[index]
TotCounts = float(TotCounts)
SigPix = float(SigPix)
qobs = float(qobs)
if twist: # pragma: no cover
anglegass = alphaf[index] # - 90
anglegass = float(anglegass)
if gauss: # pragma: no cover
C0 = TotCounts / (2 * np.pi * qobs * SigPix**2)
Ftot = 2 * np.pi * SigPix**2 * C0 * qobs
mgemag = magzpt + 0.1 + 2.5 * np.log10(exptime) - 2.5 * np.log10(Ftot)
FWHM = 2.35482 * SigPix
outlinea = "Mag: {:.2f} Sig: {:.2f} FWHM: {:.2f} ".format(
mgemag, SigPix, FWHM
)
outlineb = "q: {:.2f} angle: {:.2f} \n".format(qobs, anglegass)
outline = outlinea + outlineb
print(outline)
outline2 = "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} \n".format(
mgemag, SigPix, FWHM, qobs, anglegass
)
fout2.write(outline2)
PrintGauss(
fout1,
index + 1,
xpeak + 1,
ypeak + 1,
mgemag,
FWHM,
qobs,
anglegass,
Z,
fit,
)
else:
C0 = TotCounts / (2 * np.pi * qobs * SigPix**2)
Ftot = 2 * np.pi * SigPix**2 * C0 * qobs
mgemag = magzpt + 0.1 + 2.5 * np.log10(exptime) - 2.5 * np.log10(Ftot)
FWHM = 2.35482 * SigPix
h = np.sqrt(2) * SigPix
Re = (K ** (0.5)) * h
outlinea = "Mag: {:.2f} Sig: {:.2f} FWHM: {:.2f} ".format(
mgemag, SigPix, FWHM
)
outlineb = "Re: {:.2f} q: {:.2f} angle: {:.2f} \n".format(
Re, qobs, anglegass
)
outline = outlinea + outlineb
print(outline)
outline2 = "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} \n".format(
mgemag, SigPix, FWHM, Re, qobs, anglegass
)
fout2.write(outline2)
PrintSersic(
fout1,
index + 1,
xpeak + 1,
ypeak + 1,
mgemag,
Re,
0.5,
qobs,
anglegass,
Z,
fit,
serfit,
)
index += 1
if freesky: # pragma: no cover
skyfit = 1
else:
skyfit = 0
PrintSky(fout1, index + 1, sky, Z, skyfit)
fout1.close()
fout2.close()
# makeConstraints(consfile, len(counts))
makeConstraints(consfile, numgauss)
print(
"Done. Gaussians are stored in {}, and {} for galfit format ".format(
mgeoutfile, parfile
)
)
# returns output name
return T2
####################################################
####################################################
#####################################################
[docs]
def ReadMgePsf(psfile):
"""Reads the file created by mge2gal
returns the normalized psf.
Notes
-----
Not used in GALFITools and will be removed
in the future
"""
counts, mag, sigpix, fwhm, qobs, angle = np.genfromtxt(
psfile, skip_header=1, delimiter="", unpack=True
)
tot = counts.sum()
psfs = sigpix
normpsf = counts / tot
return psfs, normpsf
'''
def Sextractor(filimage, X, Y):
"""
Runs Sextractor and recover photometry from
the object indicated by coordinates X, Y
Notes
-----
Not used in the library and will be removed in the future
"""
sexfile = "default.sex"
CatSex = "sim.cat"
runcmd = "sextractor -c {} {} ".format(sexfile, filimage)
print(runcmd)
sp.run(
[runcmd], shell=True, stdout=sp.PIPE, stderr=sp.PIPE, universal_newlines=True
) # Run GALFIT
# KronScale = 1
# Read in sextractor sorted data
if os.stat(CatSex).st_size != 0:
(
Num,
RA,
Dec,
XPos,
YPos,
Mag,
Kron,
FluxRad,
IsoArea,
AIm,
E,
Theta,
Background,
Class,
Flag,
) = np.genfromtxt(
CatSex, delimiter="", unpack=True
) # sorted
# checking if they are numpy arrays:
# index = Mag.argsort()
# i=0
dist = 100
for idx, item in enumerate(Num):
# if(isinstance(Num,np.ndarray)):
dx = XPos[idx] - X
dy = YPos[idx] - Y
dt = np.sqrt(dx**2 + dy**2)
cflag = CheckFlag(4, Flag[idx])
if cflag is False and dt < dist:
dist = dt
sindex = idx
Num = Num[sindex]
RA = RA[sindex]
Dec = Dec[sindex]
XPos = XPos[sindex]
YPos = YPos[sindex]
Mag = Mag[sindex]
Kron = Kron[sindex]
FluxRad = FluxRad[sindex]
IsoArea = IsoArea[sindex]
AIm = AIm[sindex]
E = E[sindex]
Theta = Theta[sindex]
Background = Background[sindex]
Class = Class[sindex]
Flag = Flag[sindex]
Angle = np.abs(Theta)
# AR = 1 - E
# RKron = KronScale * AIm * Kron
# Sky = Background
XPos = round(XPos)
YPos = round(YPos)
XPos = int(XPos)
YPos = int(YPos)
else:
print("Sextractor cat file not found ")
print("Sextractor done")
return E, Angle, XPos, YPos, Background
'''
# GALFIT functions
[docs]
def PrintSky(hdl, ncomp, sky, Z, fit):
"""Print the GALFIT sky function parameters to a
file specified by the filehandle.
Parameters
----------
hdl : str
filehandler
ncomp : int
number of component
sky : float
value of sky background
Z : int
skip object. Z GALFIT component parameter
fit : int
leave parameter free/fixed during fit
Returns
-------
bool
"""
# k Check
line00 = "# Object number: {} \n".format(ncomp)
line01 = " 0) sky # Object type \n"
line02 = " 1) {} {} # sky background [ADU counts] \n".format(sky, fit)
line03 = " 2) 0.000 0 # dsky/dx (sky gradient in x) \n"
line04 = " 3) 0.000 0 # dsky/dy (sky gradient in y) \n"
line05 = " Z) {} # Skip this model in output image?\n".format(Z)
line06 = "\n"
line07a = "=========================================="
line07b = "======================================\n"
line07 = line07a + line07b
hdl.write(line00)
hdl.write(line01)
hdl.write(line02)
hdl.write(line03)
hdl.write(line04)
hdl.write(line05)
hdl.write(line06)
hdl.write(line07)
return True
[docs]
def PrintSersic(
hdl, ncomp, xpos, ypos, magser, reser, nser, axratser, angleser, Z, fit, serfit
):
"""Print the GALFIT Sersic function parameters to a
file specified by the filehandle.
Parameters
----------
hdl : str
filehandler
ncomp : int
number of component
xpos, ypos : int, int
pixel position of the component's center
magser : float
magnitud of Sersic component
reser : float
effective radius in pixels
nser : float
Sersic index
axratser : axis ratio of Sersic component
angleser : angular position of component measured from Y-axis
in degrees
Z : int
skip object. Z GALFIT component parameter
fit : int
leave parameter free/fixed during fit
Returns
-------
bool
"""
line00 = "# Object number: {} \n".format(ncomp)
line01 = " 0) sersic # Object type \n"
line02 = " 1) {:.2f} {:.2f} {} {} # position x, y [pixel] \n".format(
xpos, ypos, fit, fit
)
line03 = " 3) {:.2f} {} # total magnitude \n".format(magser, fit)
line04 = " 4) {:.2f} {} # R_e [Pixels] \n".format(reser, fit)
line05 = " 5) {} {} # Sersic exponent (deVauc=4, expdisk=1) \n".format(
nser, serfit
)
line06 = " 6) 0.0000 0 # ---------------- \n"
line07 = " 7) 0.0000 0 # ---------------- \n"
line08 = " 8) 0.0000 0 # ---------------- \n"
line09 = " 9) {:.2f} {} # axis ratio (b/a) \n".format(axratser, fit)
line10 = "10) {:.2f} {} # position angle (PA) \n".format(angleser, fit)
lineZ = " Z) {} # Skip this model in output image? \n".format(Z)
line11 = "\n"
hdl.write(line00)
hdl.write(line01)
hdl.write(line02)
hdl.write(line03)
hdl.write(line04)
hdl.write(line05)
hdl.write(line06)
hdl.write(line07)
hdl.write(line08)
hdl.write(line09)
hdl.write(line10)
hdl.write(lineZ)
hdl.write(line11)
return True
[docs]
def PrintGauss(hdl, ncomp, xpos, ypos, magass, fwhm, axratgass, anglegass, Z, fit):
"""Print the GALFIT Gauss function parameters to a
file specified by the filehandle.
Parameters
----------
hdl : str
filehandler
ncomp : int
number of component
xpos, ypos : int, int
pixel position of the component's center
magass : float
magnitud of Gaussian component
fwhm : float
Full Width Half Maximum in pixels
axratgass : axis ratio of Gauss component
anglegass : angular position of component measured from Y-axis
in degrees
Z : int
skip object. Z GALFIT component parameter
fit : int
leave parameter free/fixed during fit
Returns
-------
bool
# repeated
"""
line00 = "# Object number: {} \n".format(ncomp)
line01 = " 0) gaussian # Object type \n"
line02 = " 1) {:.2f} {:.2f} {} {} # position x, y [pixel] \n".format(
xpos, ypos, fit, fit
)
line03 = " 3) {:.2f} {} # total magnitude \n".format(magass, fit)
line04 = " 4) {:.2f} {} # FWHM [Pixels] \n".format(fwhm, fit)
line05 = " 9) {:.2f} {} # axis ratio (b/a) \n".format(axratgass, fit)
line06 = "10) {:.2f} {} # position angle (PA) [Degrees: Up=0, Left=90]\n".format(
anglegass, fit
)
line07 = " Z) {} # Skip this model in output image? (yes=1, no=0) \n".format(Z)
line08 = "\n"
hdl.write(line00)
hdl.write(line01)
hdl.write(line02)
hdl.write(line03)
hdl.write(line04)
hdl.write(line05)
hdl.write(line06)
hdl.write(line07)
hdl.write(line08)
return True
[docs]
def makeConstraints(consfile: str, numcomp: int) -> True:
"""Creates a contraints file for GALFIT with the number of components
Parameters
----------
consfile : str
the new constraints file
numcomp : int
the number of components to constraint
Returns
-------
bool
"""
fout = open(consfile, "w")
line00 = "# Component/ parameter constraint Comment \n"
line01 = "# operation (see below) range \n"
linempty = " \n"
comp = np.arange(2, numcomp + 1)
cad = "1"
for idx, item in enumerate(comp):
cad = cad + "_" + str(item)
line02 = " " + cad + " " + "x" + " " + "offset \n"
line03 = " " + cad + " " + "y" + " " + "offset \n"
fout.write(line00)
fout.write(line01)
fout.write(linempty)
fout.write(line02)
fout.write(line03)
fout.close()
return True
#############################################################################
# End of program ###################################
# ______________________________________________________________________
# /___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/___/_/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/|
# |___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|__/|
# |_|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|___|/
##############################################################################