Skip to content

Simu pa #14

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

Open
PAduverne wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Simu pa #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
268a011
add changes for potonetric conversion and possibility to insert objec…
PAduverne Oct 19, 2020
9bf314f
resolved conflicts
PAduverne Oct 21, 2020
0f35fb1
Merge branch 'master' of https://github.com/dcorre/gmadet into simu_pa
PAduverne Oct 23, 2020
ed6e2ab
rebase sim.py with current version
PAduverne Oct 23, 2020
120ef5a
minor changes and bug fix
PAduverne Oct 23, 2020
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
89 changes: 56 additions & 33 deletions gmadet/cnn/sim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,26 +24,27 @@
from gmadet.utils import rm_p, mkdir_p


def sim(datapath, filenames, Ntrans=50, size=48,
magrange=[14, 22], gain=None, magzp=30):
def sim(datapath, telescope, filenames, Ntrans=50, size=48, magrange=[14, 22],
gain=1, magzp=30, radec=None, magnitude=None):
"""Insert point sources in real images """

filenames = np.atleast_1d(filenames)

#simdir = os.path.join(datapath, "simulation")
# simdir = os.path.join(datapath, "simulation")
simdir = datapath
mkdir_p(simdir)

cutsize = np.array([size, size], dtype=np.int32)

hcutsize = cutsize // 2

#  List to store position of simulated transients
# List to store position of simulated transients
trans_pix = []
trans_wcs = []
filelist = []
maglist = []
filterlist = []
trans_count = []
cpsf1 = np.zeros((cutsize[1], cutsize[0]))

counter = 0
Expand Down Expand Up @@ -89,30 +90,45 @@ def sim(datapath, filenames, Ntrans=50, size=48,
filelist.append(os.path.abspath(newfile))

filterlist.append(band)

pos[j] = np.random.random_sample(
2) * (imsize - cutsize) + cutsize / 2.0
# store positions
ra, dec = w.wcs_pix2world(pos[j][1], pos[j][0], 1)
trans_pix.append(pos[j])
trans_wcs.append([ra, dec])

# get pixels indexes in the image
ipos = pos[j].astype(int)
# extract subimage centered on position of the object and
# store in cima1
# same for weight maps
iposrange = np.s_[
ipos[0] - hcutsize[0]: ipos[0] + hcutsize[0],
ipos[1] - hcutsize[1]: ipos[1] + hcutsize[1],
]

if radec is None:
pos[j] = np.random.random_sample(
2) * (imsize - cutsize) + cutsize / 2.0
# store positions
ra, dec = w.wcs_pix2world(pos[j][1], pos[j][0], 1)
trans_pix.append(pos[j])
trans_wcs.append([ra, dec])

# get pixels indexes in the image
ipos = pos[j].astype(int)
# extract subimage centered on position of the object and
# store in cima1
# same for weight maps
iposrange = np.s_[
ipos[0] - hcutsize[0]: ipos[0] + hcutsize[0],
ipos[1] - hcutsize[1]: ipos[1] + hcutsize[1],
]
else:
ra, dec = radec[j][0], radec[j][1]
# print('ra, dec = ', ra, dec)
pos[j] = w.all_world2pix(ra, dec, 0)
# print(pos[j])
trans_pix.append(pos[j])
trans_wcs.append([ra, dec])
# get pixels indexes in the image
ipos = pos[j].astype(int)
# print(ipos)
# extract subimage centered on position of the object
# and store in cima1
# same for weight maps
iposrange = np.s_[
ipos[1] - hcutsize[1]: ipos[1] + hcutsize[1],
ipos[0] - hcutsize[0]: ipos[0] + hcutsize[0]]
# Select the PSF corresponding to the image area,
# in case there are more than one PSF estimated per axis
if nb_psf_snaps == 1:
psf1 = psfs1
else:
#  get position with respect to number of PSF snapshots
# get position with respect to number of PSF snapshots
ppos = (pos[j] * posfac).astype(int)
psf1 = psfs1[ppos[0], ppos[1]]
# step1 is psf_samp parameter from psfex, used in mat1 and 2
Expand All @@ -132,13 +148,17 @@ def sim(datapath, filenames, Ntrans=50, size=48,

# define the object magnitude using this random number and
# predefined ranges
mag = np.random.uniform(
low=magrange[0], high=magrange[1], size=(1,))
# convert the magnitude in ADU using the zeropoint magnitude.
# Note that the zeropoint magnitude is define as 30,
# so did not care of the exact value.
# simply needed to draw random magnitudes. We could estimate
# the proper one for our telescopes
if magnitude is None:
mag = np.random.uniform(
low=magrange[0], high=magrange[1], size=(1,))
# convert the magnitude in ADU using the zeropoint
# magnitude.
# Note that the zeropoint magnitude is define as 30,
# so did not care of the exact value.
# simply needed to draw random magnitudes.
# We could estimate the proper one for our telescopes
else:
mag = np.array([magnitude[j]])
maglist.append(mag[0])
amp1 = np.exp(0.921034 * (magzp - mag))

Expand All @@ -150,9 +170,11 @@ def sim(datapath, filenames, Ntrans=50, size=48,
noisy_object1 = noisy_object1 / gain

ima1[iposrange] += noisy_object1
count = np.sum(noisy_object1)
trans_count.append(count)

hdusi1[0].data = ima1
#  Write new fits file
# Write new fits file
hdusi1.writeto(newfile, overwrite=True)

hdusi1.close()
Expand All @@ -171,6 +193,7 @@ def sim(datapath, filenames, Ntrans=50, size=48,
wcspos[:, 1],
maglist,
filterlist,
trans_count
],
names=[
"idx",
Expand All @@ -180,12 +203,12 @@ def sim(datapath, filenames, Ntrans=50, size=48,
"RA",
"Dec",
"mag",
"filter"],
"filter",
"count_ADU"],
)
table.write(
os.path.join(simdir, "simulated_objects.list"),
format="ascii.commented_header",
overwrite=True,
)
return table

193 changes: 68 additions & 125 deletions gmadet/phot_conversion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,76 +187,53 @@ def SDSS2Johnson(band, SDSS_Table):
R_band = ["R"] # found in header
I_band = ["I"]
B_band = ["B"]

if band in V_band:
coefficients = [-0.016, -0.573]
SDSS_Table["g-r"] = SDSS_Table["gmag"] - SDSS_Table["rmag"]
SDSS_Table["VMag"] = SDSS_Table["gmag"] + \
poly(SDSS_Table["g-r"], coefficients)
SDSS_Table["calib_err"] = (
np.sqrt(
(0.002 / 0.573) ** 2
+ (0.002 / 0.016) ** 2
+ (SDSS_Table["e_gmag"] / SDSS_Table["gmag"]) ** 2
+ (SDSS_Table["e_rmag"] / SDSS_Table["rmag"]) ** 2
)
* SDSS_Table["VMag"]
)

if band in R_band:

if filter in V_band:
coefficients = [- 0.016, -0.573]
SDSS_Table['g-r'] = SDSS_Table['gmag'] - SDSS_Table['rmag']
SDSS_Table["VMag"] = SDSS_Table["gmag"] + poly(SDSS_Table['g-r'],
coefficients)
SDSS_Table["calib_err"] = np.sqrt((0.002 * SDSS_Table['g-r'])**2 +
0.573**2 * SDSS_Table['e_rmag']**2 +
0.427**2 * SDSS_Table['e_gmag']**2 +
0.002**2)

if filter in R_band:
coefficients = [0.152, -0.257]
SDSS_Table["r-i"] = SDSS_Table["rmag"] - SDSS_Table["imag"]
SDSS_Table["RMag"] = SDSS_Table["rmag"] + \
poly(SDSS_Table["r-i"], coefficients)
SDSS_Table["calib_err"] = (
np.sqrt(
(0.004 / 0.257) ** 2
+ (0.002 / 0.152) ** 2
+ (SDSS_Table["e_rmag"] / SDSS_Table["rmag"]) ** 2
+ (SDSS_Table["e_imag"] / SDSS_Table["imag"]) ** 2
)
* SDSS_Table["RMag"]
)

if band in I_band:
SDSS_Table['r-i'] = SDSS_Table['rmag'] - SDSS_Table['imag']
SDSS_Table["RMag"] = SDSS_Table["rmag"] + poly(SDSS_Table['r-i'],
coefficients)
SDSS_Table["calib_err"] = np.sqrt(0.743**2 * SDSS_Table['e_rmag']**2 +
0.257**2 * SDSS_Table['e_imag']**2 +
0.002**2 + SDSS_Table['r-i']**2 * 0.004**2)


if filter in I_band:
coefficients = [-0.394, -0.409]
SDSS_Table["i-z"] = SDSS_Table["imag"] - SDSS_Table["zmag"]
SDSS_Table["IMag"] = SDSS_Table["imag"] + \
poly(SDSS_Table["i-z"], coefficients)
SDSS_Table["calib_err"] = (
np.sqrt(
(0.006 / 0.409) ** 2
+ (0.002 / 0.394) ** 2
+ (SDSS_Table["e_zmag"] / SDSS_Table["zmag"]) ** 2
+ (SDSS_Table["e_imag"] / SDSS_Table["imag"]) ** 2
)
* SDSS_Table["IMag"]
)

if band in B_band:
SDSS_Table['i-z'] = SDSS_Table['imag'] - SDSS_Table['zmag']
SDSS_Table["IMag"] = SDSS_Table["imag"] + poly(SDSS_Table['i-z'],
coefficients)
SDSS_Table["calib_err"] = np.sqrt(0.591**2 * SDSS_Table['e_imag']**2 +
0.409**2 * SDSS_Table['e_zmag']**2 +
0.002**2 + SDSS_Table['i-z']**2 * 0.006**2)

if filter in B_band:
coefficients = [0.219, 0.312]
SDSS_Table["g-r"] = SDSS_Table["gmag"] - SDSS_Table["rmag"]
SDSS_Table["BMag"] = SDSS_Table["gmag"] + \
poly(SDSS_Table["g-r"], coefficients)
SDSS_Table["calib_err"] = (
np.sqrt(
(0.002 / 0.219) ** 2
+ (0.003 / 0.312) ** 2
+ (SDSS_Table["e_gmag"] / SDSS_Table["gmag"]) ** 2
+ (SDSS_Table["e_rmag"] / SDSS_Table["rmag"]) ** 2
)
* SDSS_Table["BMag"]
)

SDSS_Table['g-r'] = SDSS_Table['gmag'] - SDSS_Table['rmag']
SDSS_Table["BMag"] = SDSS_Table["gmag"] + poly(SDSS_Table['g-r'],
coefficients)
SDSS_Table["calib_err"] = np.sqrt(1.312**2 * SDSS_Table['e_gmag']**2 +
0.312**2 * SDSS_Table['e_rmag']**2 + 0.002**2 +
SDSS_Table['g-r']**2 * 0.003**2)

return SDSS_Table


def PS2Johnson(band, PS_Table):
"""
Give the transformation laws to go from Pan-STARRS photometric system
to Johnson photometric system
Transformation given by http://www.sdss3.org/dr8/algorithms/
sdssUBVRITransform.php
Transformation given by https://arxiv.org/pdf/1706.06147.pdf

parameters: band: filter used to get the image
PS_Table: astropy.table with information from PS for
Expand All @@ -265,69 +242,35 @@ def PS2Johnson(band, PS_Table):
returns: astropy.table with informations from PS and the computation
for the filter of the image
"""
V_band = ["V"] # Different names that can be
R_band = ["R"] # found in header
I_band = ["I"]
B_band = ["B"]

if band in V_band:
coefficients = [-0.016, -0.573]
PS_Table["g-r"] = PS_Table["gmag"] - PS_Table["rmag"]
PS_Table["VMag"] = PS_Table["gmag"] + \
poly(PS_Table["g-r"], coefficients)
PS_Table["calib_err"] = (
np.sqrt(
(0.002 / 0.573) ** 2
+ (0.002 / 0.016) ** 2
+ (PS_Table["e_gmag"] / PS_Table["gmag"]) ** 2
+ (PS_Table["e_rmag"] / PS_Table["rmag"]) ** 2
)
* PS_Table["VMag"]
)

if band in R_band:
coefficients = [0.152, -0.257]
PS_Table["r-i"] = PS_Table["rmag"] - PS_Table["imag"]
PS_Table["RMag"] = PS_Table["rmag"] + \
poly(PS_Table["r-i"], coefficients)
PS_Table["calib_err"] = (
np.sqrt(
(0.004 / 0.257) ** 2
+ (0.002 / 0.152) ** 2
+ (PS_Table["e_rmag"] / PS_Table["rmag"]) ** 2
+ (PS_Table["e_imag"] / PS_Table["imag"]) ** 2
)
* PS_Table["RMag"]
)

if band in I_band:
coefficients = [-0.394, -0.409]
PS_Table["i-z"] = PS_Table["imag"] - PS_Table["zmag"]
PS_Table["IMag"] = PS_Table["imag"] + \
poly(PS_Table["i-z"], coefficients)
PS_Table["calib_err"] = (
np.sqrt(
(0.006 / 0.409) ** 2
+ (0.002 / 0.394) ** 2
+ (PS_Table["e_zmag"] / PS_Table["zmag"]) ** 2
+ (PS_Table["e_imag"] / PS_Table["imag"]) ** 2
)
* PS_Table["IMag"]
)

if band in B_band:
coefficients = [0.219, 0.312]
PS_Table["g-r"] = PS_Table["gmag"] - PS_Table["rmag"]
PS_Table["BMag"] = PS_Table["gmag"] + \
poly(PS_Table["g-r"], coefficients)
PS_Table["calib_err"] = (
np.sqrt(
(0.002 / 0.219) ** 2
+ (0.003 / 0.312) ** 2
+ (PS_Table["e_gmag"] / PS_Table["gmag"]) ** 2
+ (PS_Table["e_rmag"] / PS_Table["rmag"]) ** 2
)
* PS_Table["BMag"]
)


V_band = ['V'] #Different names that can be
R_band = ['R'] #found in header
I_band = ['I']
B_band = ['B']

if filter in R_band:
coefficients = [-0.163, -0.086, -0.061]
PS_Table['g-r'] = PS_Table['gmag'] - PS_Table['rmag']
PS_Table["RMag"] = PS_Table["rmag"] + poly(PS_Table['g-r'], coefficients)
PS_Table["calib_err"] = 0.041

if filter in V_band:
coefficients = [-0.020, -0.498, -0.008]
PS_Table['g-r'] = PS_Table['gmag'] - PS_Table['rmag']
PS_Table["VMag"] = PS_Table["gmag"] + poly(PS_Table['g-r'], coefficients)
PS_Table["calib_err"] = 0.032

if filter in B_band:
coefficients = [0.199, 0.540, 0.016]
PS_Table['g-r'] = PS_Table['gmag'] - PS_Table['rmag']
PS_Table["BMag"] = PS_Table["gmag"] + poly(PS_Table['g-r'], coefficients)
PS_Table["calib_err"] = 0.056

if filter in I_band:
coefficients = [-0.387, -0.123, -0.03]
PS_Table['g-r'] = PS_Table['gmag'] - PS_Table['rmag']
PS_Table["IMag"] = PS_Table["imag"] + poly(PS_Table['g-r'], coefficients)
PS_Table["calib_err"] = 0.054

return PS_Table