Commit a6ea0ab6 authored by Peter Lustig's avatar Peter Lustig

Added fits output.

parent 30cfdcf7
......@@ -25,38 +25,21 @@ ipython = get_ipython()
import myfunctions as mfct
from time import clock
from astropy.io import fits
if '__IPYTHON__' in globals():
ipython.magic('load_ext autoreload')
ipython.magic('autoreload 2')
ipython.magic('matplotlib tk')
#%load_ext autoreload
#%autoreload 2
#%matplotlib tk
comm = MPI.COMM_WORLD
rank, size = comm.Get_rank(), comm.Get_size()
plt.ion()
def Wait():
while True:
a=input("[c]ontinue? ")
if a=="c":
break
# DATA_DIR = "/data/NIKA/Reduced/G2_kidpar_20171025s41_v2_LP_md_recal_calUranus/v_1/"
DATA_DIR = "/home/peter/Dokumente/Uni/Paris/Stage/N2CLS/ForMarseille/HLS091828_common_mode_one_block/"
# def scanfile_to_filenames(scan_file):
# scan_list = ascii.read(scan_file, format="fast_no_header", names=['SCAN'])
# filenames = [os.path.join(DATA_DIR, scan, 'map.fits') for scan in scan_list['SCAN']]
# return filenames
def add_axis(name, range, bins, unit=None, i_axe=3, log=False):
"""Define a dictionnary for additionnal wcs axes (linear or log)"""
......@@ -172,13 +155,7 @@ def completness_worker(shape, wcs, sources, fake_sources, min_threshold=2, max_t
# covering ALL possible SNR, otherwise loose flux, or cap the thresholds...
fake_snr = np.ma.array(sources[fake_sources['find_peak']]['SNR'],
mask=fake_sources['find_peak'].mask)
#print('-----------------------')
#print(fake_sources['find_peak'].mask)
#print('-----------------------')
#print(fake_sources['find_peak'])
#print('-----------------------')
#print(fake_snr)
#print('-----------------------')
# As we are interested by the cumulative numbers, keep all inside the upper pixel
fake_snr[fake_snr > max_threshold] = max_threshold
......@@ -186,14 +163,6 @@ def completness_worker(shape, wcs, sources, fake_sources, min_threshold=2, max_t
# TODO: Consider keeping all pixels information in fake_source and source...
# This would imply to do only a simple wcs_threshold here...
xx, yy, zz = wcs.wcs_world2pix(fake_sources['ra'], fake_sources['dec'], fake_snr.filled(min_threshold), 0)
##
## Warum werden jetzt ?nicht detektierte fake quellen? mit niedrigster schwelle aufgefüllt?
## Ich muss unbedingt rausbekommen wo diese maskierte spalte herkommt und was sie bedeutet
#print('-----------------------')
#print(xx, yy, zz)
#print('-----------------------')
#print(shape)
#print('-----------------------')
# Number of fake sources recovered
_completness, _ = np.histogramdd(np.asarray([xx, yy, zz]).T + 0.5, bins=np.asarray(shape),
......@@ -296,7 +265,7 @@ def completness_purity(flux, wcs=None, shape=None, nsources=8**2, within=(1 / 4,
def completness_purity_2(flux, nsim=2, wcs=None, shape=None, nsources=8**2, within=(1 / 4, 3 / 4), jk_map=None):
"""Compute completness map for a given flux"""
print(flux)
#print(flux)
# wcs_celestial = wcs.celestial
# Lower and upper edges ... take the center of the pixel for the upper edge
......@@ -354,6 +323,34 @@ def Plot_CompPur(completness, purity, threshold,nsim=None):
plt.show(block=True)
def AddParametersToHeader(header, threshold, nsim, method=None):
header['thresh'] = threshold
header['nsim'] = nsim
if method is not None:
header['method'] = method
def WriteResult(filename, threshold, completness, purity, nsim, wcs4D):
allhducomp, allhdupur = [],[]
for i, thresh in enumerate(threshold):
hducomp = fits.ImageHDU(data=completness[i], name='completness_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hducomp.header, thresh, nsim)
allhducomp.append(hducomp)
hdupur=fits.ImageHDU(data=purity[i], name='purity_thresh_{:.1f}'.format(float(threshold[i])))
AddParametersToHeader(hdupur.header, thresh, nsim)
allhdupur.append(hdupur)
wcs2D = wcs4D.sub([1,2])
hdu = fits.PrimaryHDU(header=wcs2D.to_header())
hdul = [hdu]+allhducomp+allhdupur
hdul = fits.HDUList(hdul)
hdul.writeto(filename, overwrite=True)
plt.close('all')
filenames = list(Path(DATA_DIR).glob('../*/map.fits'))
......@@ -397,12 +394,16 @@ assert np.all(np.abs(wcs_flux.all_pix2world(np.arange(flux_bins+1)-0.5, 0) - flu
FAKE_SOURCES = []
DETECTED_SOURCES = []
#wcsheader = wcs_4D.to_header()
#hdu = fits.PrimaryHDU(header=wcsheader)
#hdu.writeto('wcs.fits')
#print('written')
# DEBUG :
# flux, nsources, within, wcs, shape, nsim, jk_filenames = 10*u.mJy, 8**2, (0, 1), wcs_4D.sub([1, 2, 3]), shape_4D[0:3], np.multiply(*shape_4D[0:2]) * 100, filenames
# This is a single run check for a single flux
nsim = 5
nsim = 4
#ncores = 3
allflux = [10*u.mJy]
......@@ -446,14 +447,12 @@ for irun, flux in enumerate(allflux):
all_comp_n = np.sum(all_comp_n, axis=0)
all_pur = np.sum(all_pur, axis=0)
all_pur_n = np.sum(all_pur_n, axis=0)
print(all_comp_n.shape)
print('------------')
with np.errstate(divide='ignore', invalid='ignore'):
completness = all_comp / all_comp_n[..., np.newaxis]
purity = all_pur / all_pur_n
WriteResult('test.fits', threshold, completness, purity, nsim, wcs_4D)
Plot_CompPur(completness, purity, threshold, nsim=nsim)
# Save result !!
......
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