Commit 6b769b02 authored by LUSTIG Peter's avatar LUSTIG Peter

added some files...

parent fabef1f8
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from pathlib import Path
from matplotlib.figure import Figure
class myfigure(Figure):
def save(self, filename, outdir='/home/peter/theoryplots'):
outpath = Path(outdir)
formats = ['.png', '.svg', '.pdf']
for format in formats:
_filename = outpath / (filename + format)
self.savefig(str(_filename))
# %matplotlib tk
shape = (501, 501)
fwhm_pix = 6.25
circle_radius = 3 * fwhm_pix
idy, idx = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
mean = (250, 251)
pixdist2 = np.square(idy-mean[0]) + np.square(idx-mean[1])
mask = pixdist2 < (circle_radius ** 2)
phdu = fits.PrimaryHDU()
ihdu = fits.ImageHDU(data=mask.astype(int), name='HLSmask')
ihdu.header['radius'] = circle_radius
ihdu.header['xcenter'] = mean[0]
ihdu.header['ycenter'] = mean[1]
# Load data for plot
maskdir = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify')
with fits.open(maskdir / 'snrmasks.fits') as hdul:
areamask = (hdul['SNRMASK0'].data).astype(bool)
wcs = WCS(hdul['SNRMASK0'].header)
newfile = '/home/peter/Dokumente/Uni/Paris/Stage/NicePlots/newimage.fits'
SNRimg = fits.getdata(newfile, 'SNR')
# %%
if 1:
bandadd = ' 1mm'
plt.close('all')
import matplotlib.pylab as pl
from matplotlib.colors import ListedColormap
cmap = pl.cm.get_cmap('binary')
my_cmap = cmap(np.arange(cmap.N))
my_cmap[:,-1] = np.linspace(0, 1, cmap.N)
my_cmap = ListedColormap(my_cmap)
globmask = np.array(areamask | mask, dtype=bool)
labelsize = 23
titlesize = 28
ticklabelsize = 18
plotwidth = .35
plotheight = .8
# %matplotlib tk
# Hitmap, Variancemap with colorbar try
fig = plt.figure(figsize=(14, 6))
# ax1 = fig.add_subplot(121, projection=nm.wcs)
ax1 = fig.add_axes([.07, .1, plotwidth, plotheight], projection=nm.wcs)
imhits = ax1.imshow(SNRimg, origin='lower')
ra = ax1.coords[0]
dec = ax1.coords[1]
# ax1.set_title('Nhits{}'.format(bandadd), fontsize=25, y=1.02)
ax1.set_title('HLS091828{}'.format(bandadd), fontsize=titlesize, y=1.02)
ra.set_axislabel('Ra', minpad=0.4, fontsize=labelsize)
ra.set_ticklabel(fontsize=ticklabelsize)
dec.set_axislabel('Dec', minpad=-0.5, fontsize=labelsize)
dec.set_ticklabel(fontsize=ticklabelsize)
# %matplotlib tk
# fig.colorbar(imhits, shrink=0.7)
# (ax1x0, ax1y0), (ax1x0, ax1y0) =
# ax2 = fig.add_subplot(122, projection=nm_mf.wcs)
ax2 = fig.add_axes([0.93-plotwidth, 0.1, plotwidth, plotheight], projection=nm.wcs)
snrvmin, snrvmax = -3, 5
imstddev = ax2.imshow(SNRimg, vmin=snrvmin, vmax=snrvmax, origin='lower')
ax2.imshow(globmask, origin='lower', alpha=0.6, cmap=my_cmap)
ra = ax2.coords[0]
dec = ax2.coords[1]
# ax2.set_title('1 / Var(pix)', fontsize=25, y=1.02)
ax2.set_title('SNR', fontsize=titlesize, y=1.02)
ra.set_axislabel('Ra', minpad=0.4, fontsize=labelsize)
ra.set_ticklabel(fontsize=ticklabelsize)
dec.set_ticklabel_position('r')
dec.set_ticklabel(fontsize=ticklabelsize)
dec.set_ticklabel_visible(True)
pos1 = ax1.get_position()
pos2 = ax2.get_position()
caxwidth = 0.02
caxmean = (pos1.x1 + pos2.x0) / 2
cax = fig.add_axes([caxmean-.5*caxwidth, pos1.y0, caxwidth,
(pos1.y1-pos1.y0)])
pos1.height
cbar = fig.colorbar(imhits, cax=cax)
# cbar.set_label('', fontsize=20)
cax.yaxis.set_label_position('left')
cax.yaxis.set_tick_params(labelsize=ticklabelsize)
cax.yaxis.set_ticks_position('left')
cax.set_ylabel('Flux [mJy]', fontsize=labelsize, rotation=270,
va='center', labelpad=14.7)
cax2 = cax.twinx()
cax2labelpos = cax2.get_yticks()
nlabels = 9
cax2labelpos = np.linspace(0, 1, nlabels)
cax2labels = np.linspace(snrvmin, snrvmax, nlabels)
cax2.set_yticks(cax2labelpos)
cax2formattedlabels = ['{:.0f}'.format(label) for label in cax2labels]
cax2.set_yticklabels(cax2formattedlabels)
cax2.yaxis.set_tick_params(labelsize=ticklabelsize)
cax2.set_ylabel('SNR', fontsize=labelsize, labelpad=15)
# %%
#%matplotlib tk
# SaveFigure(fig, )
'''
fig = plt.figure(figsize=(12, 6))
ax.set_title('Detection Mask', fontsize=30, y=1.02)
ax.imshow(globmask, origin='lower', alpha=0.5, cmap=my_cmap)
ra = ax.coords[0]
dec = ax.coords[1]
print(type(ra))
ra.set_axislabel('Ra', minpad=0.4, fontsize=25)
ra.set_ticklabel(fontsize=20)
dec.set_axislabel('Dec', minpad=0.17, fontsize=25)
dec.set_ticklabel(fontsize=20)
plt.subplots_adjust(left=0.16)
plt.show()
'''
if 0:
import matplotlib.pylab as pl
from matplotlib.colors import ListedColormap
cmap = pl.cm.get_cmap('binary')
my_cmap = cmap(np.arange(cmap.N))
my_cmap[:,-1] = np.linspace(0, 1, cmap.N)
my_cmap = ListedColormap(my_cmap)
#sys.exit()
maskfile = 'shit...'
fname = 'jgfsjdg'
hdul = fits.open(fname)
areamask = hdul['MATCHMASK'].data != 0
globmask = np.array(areamask | mask, dtype=bool)
realfile = '/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits'
from nikamap import NikaMap
nm = NikaMap.read(realfile)
mf_nm = nm.match_filter(nm.beam)
std = mf_nm.check_SNR()
mf_nm.uncertainty.array *= std
# mf_nm.mask = np.array(mf_nm.mask | mask, dtype=bool)
# plt.figure()
im = mf_nm.plot_SNR()
ax = plt.gca()
ax.set_title('Detection Mask', fontsize=30, y=1.02)
ax.imshow(globmask, origin='lower', alpha=0.5, cmap=my_cmap)
ra = ax.coords[0]
dec = ax.coords[1]
print(type(ra))
ra.set_axislabel('Ra', minpad=0.4, fontsize=25)
ra.set_ticklabel(fontsize=20)
dec.set_axislabel('Dec', minpad=0.17, fontsize=25)
dec.set_ticklabel(fontsize=20)
plt.subplots_adjust(left=0.16)
plt.show()
# fits.HDUList([phdu, ihdu]).writeto('hlsmask3.fits', overwrite=False)
......@@ -159,7 +159,8 @@ directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/monte
# 'PhotometryNewrealisation3.5_MultTables.fits')
directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/realdetection'
outname = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/AllRealdetections.fits'
directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/fluxboost'
outname = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/allphotometry_3.3.fits'
key = '*'
# Combine(directory, key, outname)
MergeToOneTable(directory, key, outname, overwrite=False)
......
......@@ -542,6 +542,7 @@ class FluxArea(CompletenessArea):
center = np.zeros(len(bins)-1)
abs_mean = np.zeros(len(bins)-1)
abs_mean_var = np.zeros(len(bins)-1)
abs_var = np.zeros(len(bins)-1)
# fluxcenter = np.zeros(len(bins)-1)
for i in range(0, len(bins)-1):
......@@ -561,6 +562,7 @@ class FluxArea(CompletenessArea):
mean_var[i] = var[i] / len(tmp)
abs_mean[i] = np.nanmean(_det_flux)
abs_mean_var[i] = np.nanvar(_det_flux) / len(tmp)
abs_var[i] = np.nanvar(_det_flux)
else:
perc[i] = np.repeat(np.nan, len(percentiles))
mean[i] = np.nan
......@@ -568,6 +570,7 @@ class FluxArea(CompletenessArea):
mean_var[i] = np.nan
abs_mean[i] = np.nan
abs_mean_var[i] = np.nan
abs_var[i] = np.nan
self.snrbins = snrbins
self.fluxpercentiles = perc
......@@ -575,6 +578,7 @@ class FluxArea(CompletenessArea):
self.fluxmean_unc = np.sqrt(mean_var)
self.fluxmean_absolute = abs_mean
self.fluxmean_absolute_unc = np.sqrt(abs_mean_var)
self.flux_absolute_unc = np.sqrt(abs_var)
self.flux_unc = np.sqrt(var)
self.percentile_center = center
......
......@@ -21,7 +21,7 @@ with fits.open(maskfilename) as hdul:
# Load SNR masks
anothermaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify/'
'snrmasks_fine_new.fits')
'snrmasks_fine__newrealisation.fits')
snrmasks = []
with fits.open(anothermaskfile) as hdul:
i = 0
......
......@@ -24,7 +24,7 @@ with fits.open(maskfilename) as hdul:
# Load SNR masks
anothermaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify/'
'snrmasks_fine_new.fits')
'snrmasks_fine_newrealisation.fits')
snrmasks = []
with fits.open(anothermaskfile) as hdul:
i = 0
......@@ -111,8 +111,9 @@ if __name__ == '__main__':
len(comptbl)
comptbl
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
phottble.write(here / '1mm_photfit_results.fits', overwrite=True)
comptbl.write(here / '1mm_compfit_results.fits', overwrite=True)
phottble.write(here / '1mm_photfit_results_newres.fits', overwrite=True)
comptbl.write(here / '1mm_compfit_results_newres.fits', overwrite=True)
# phottble
if 0:
......
......@@ -23,7 +23,7 @@ with fits.open(maskfilename) as hdul:
# Load SNR masks
anothermaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify/'
'snrmasks_fine_new.fits')
'snrmasks_fine_newrealisation.fits')
snrmasks = []
with fits.open(anothermaskfile) as hdul:
i = 0
......@@ -98,29 +98,31 @@ rc = CountsAreaEvaluation(
abs_counts = rc.Counts()
plt.imshow(snrmasks[8] | hlsmask)
# matplotlib tk
plt.imshow(~snrmasks[8]+~snrmasks[7]*2)
np.sum(~(snrmasks[7] | hlsmask))
plt.imshow()
im = np.zeros(snrmasks[0].shape)
plt.close('all')
# fig = plt.figure()
fig, axes = plt.subplots(1, len(snrmasks))
for i, mask in enumerate(snrmasks):
_mask = (~mask).astype(int)
axes[i].imshow(_mask)
_mask *= (i+1)
print(np.amax(_mask))
im += _mask
print(np.sum(~(mask | hlsmask)))
# %%
if 0:
plt.imshow(snrmasks[8] | hlsmask)
# matplotlib tk
plt.imshow(~snrmasks[8]+~snrmasks[7]*2)
np.sum(~(snrmasks[7] | hlsmask))
# plt.imshow()
im = np.zeros(snrmasks[0].shape)
plt.close('all')
# fig = plt.figure()
fig, axes = plt.subplots(1, len(snrmasks))
for i, mask in enumerate(snrmasks):
_mask = (~mask).astype(int)
axes[i].imshow(_mask)
_mask *= (i+1)
print(np.amax(_mask))
im += _mask
print(np.sum(~(mask | hlsmask)))
# %matplotlib tk
plt.figure()
plt.imshow(im, origin='lower')
# plt.figure()
# plt.imshow(im, origin='lower')
pur = Purity(abs_counts, false_counts, threshbins[0:-1])
fc.Counts()
fc.threshold_edges
# fc.Counts()
# fc.threshold_edges
# %%
if __name__ == '__main__':
......@@ -150,7 +152,7 @@ if __name__ == '__main__':
t = vstack([t, tt])
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
t.write(here / '1mm_purity_table.fits', overwrite=True)
t.write(here / '1mm_purity_table_newres.fits', overwrite=True)
if 0:
......
import numpy as np
import matplotlib.pyplot as plt
from nikamap import NikaMap
from astropy.io import fits
from pathlib import Path
def plot(data, wcs=None, ax=None, steps=None):
if ax is None:
fig = plt.figure(figsize=(8, 7))
# ax = fig.add_subplot(111, projection=wcs)
ax = fig.add_axes([.15, .1, .6, .8], projection=wcs)
img = ax.imshow(data, origin='lower')
ra = ax.coords[0]
dec = ax.coords[1]
ax.set_title('Evaluation Masks'.format(), fontsize=25, y=1.02)
ra.set_axislabel('Ra', minpad=0.4, fontsize=20)
ra.set_ticklabel(fontsize=15)
dec.set_axislabel('Dec', minpad=-0.5, fontsize=20)
dec.set_ticklabel(fontsize=15)
if steps is not None:
# fig2 = plt.figure()
# ax2 = fig2.add_subplot(111)
ax2 = fig.add_axes([.8, .1, .1, .8])
# steps = np.repeat(steps, 2)
ax2.imshow((np.array([np.arange(len(steps)), np.arange(len(steps))])
.T[::-1]), aspect=2)
for i, step in enumerate(steps):
ax2.text(.5, i, '{}'.format(step), ha='center', va='center',
fontsize=20, color='r') # np.repeat(i/len(steps), 3))
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_xticklabels([])
ax2.set_ylabel('SNR Percentile', fontsize=20)
ax2.yaxis.set_label_position("right")
orimgfile = ('/home/peter/Dokumente/Uni/Paris/Stage/data/map.fits')
nm = NikaMap.read(orimgfile)
hlsmaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'masks.fits')
hlsmask = np.array(fits.getdata(hlsmaskfile, 'HLS2FWHM'), dtype=bool)
# %matplotlib tk
# plot(1/nm.uncertainty.array**2, wcs=nm.wcs)
stddev = nm.uncertainty.array
stddev = fits.getdata('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'newuncertainty.fits')
steps = np.array([10, 30, 50, 70, 100])
steps = np.arange(10, 101, 10)
steps
steps = np.insert(steps, 0, 5)
steps
# steps = np.array([10, 50, 70, 100])
thresholds = np.nanpercentile(stddev, steps)
# masks = stddev < np.moveaxis(np.atleast_3d(thresholds), 1, 0)
masks = stddev < thresholds[:, np.newaxis, np.newaxis]
masks = np.sum(masks, axis=0, dtype=float)
masks[nm.mask] = np.nan
masks[hlsmask] = np.nan
# plt.close('all')
plot(masks, steps=steps, wcs=nm.wcs)
plt.show(block=True)
# %matplotlib tk
masks = (~(masks == np.arange(len(thresholds)+1)[:, np.newaxis, np.newaxis]))[::-1]
masks.shape
# %matplotlib tk
perc_limits = np.insert(steps, 0, 0)
lowlims = perc_limits[0:-1]
highlims = perc_limits[1:]
hdul = fits.HDUList([fits.PrimaryHDU()])
header = nm.wcs.to_header()
for i, (mask, lowlim, highlim) in enumerate(zip(masks, lowlims, highlims)):
_header = header.copy()
_header['snrlow'] = lowlim
_header['snrhigh'] = highlim
hdul.append(fits.ImageHDU(data=mask.astype(int), header=_header,
name='snrmask{}'.format(i)))
# plt.figure()
# plt.title('{} - {}, {}'.format(lowlim, highlim, i), y=1.02, fontsize=25)
# plt.imshow(mask, origin='lower')
# plt.show()
# nmmask = NikaMap(data=np.ones(mask.shape), wcs=nm.wcs, mask=mask)
# print(nmmask.surface().to(u.arcmin**2))
# nmmask.plot()
outdir = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify')
# hdul.writeto(outdir / 'snrmasks_fine_newrealisation.fits')
'''
coverage = np.sum(~masks, axis=(1, 2))
conversion = (u.pix.to(u.arcmin, equivalencies=nm._pixel_scale)) ** 2
coverage = coverage * conversion * u.arcmin ** 2
print(coverage)
'''
# plt.figure()
# plt.plot(1/nm.uncertainty.array[250,:]**2)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment