generated from erosita/uds
229 lines
6.6 KiB
Python
Executable File
229 lines
6.6 KiB
Python
Executable File
#!/usr/bin/env python
|
|
|
|
__author__ = "Roman Krivonos"
|
|
__copyright__ = "Space Research Institute (IKI)"
|
|
|
|
import numpy as np
|
|
import numpy.ma as ma
|
|
import pandas as pd
|
|
from astropy.wcs import WCS
|
|
from astropy import wcs
|
|
from astropy.io import fits
|
|
from astropy.table import Table, Column
|
|
import matplotlib.pyplot as plt
|
|
import math, sys, os
|
|
import pickle
|
|
|
|
from astropy.coordinates import SkyCoord # High-level coordinates
|
|
from astropy.coordinates import ICRS, Galactic, FK4, FK5 # Low-level frames
|
|
from astropy.coordinates import Angle, Latitude, Longitude # Angles
|
|
import astropy.units as u
|
|
|
|
from sklearn.linear_model import LinearRegression
|
|
from sklearn.linear_model import HuberRegressor
|
|
from sklearn.linear_model import RANSACRegressor
|
|
from sklearn.linear_model import TheilSenRegressor
|
|
from sklearn.model_selection import cross_val_score
|
|
from sklearn.model_selection import RepeatedKFold
|
|
|
|
from astropy.stats import sigma_clip
|
|
from astropy.stats import sigma_clipped_stats
|
|
|
|
from scipy.stats import norm
|
|
from scipy.stats import describe
|
|
from scipy.stats import sem
|
|
|
|
|
|
from numpy import absolute
|
|
from numpy import arange
|
|
|
|
from ridge.utils import *
|
|
from ridge.config import *
|
|
|
|
enkey = sys.argv[1]
|
|
key = sys.argv[2]
|
|
#key="ALL"
|
|
fn='detcnts.{}.{}.resid.fits'.format(enkey,key)
|
|
|
|
d = fits.getdata(proddir+fn)
|
|
df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
|
|
#print(df.columns)
|
|
|
|
#df = df.query('abs(LAT) < {} & abs(LON) < {}'.format(5,5))
|
|
|
|
n_bins = 80
|
|
minmax=[-300,800]
|
|
sigma=3
|
|
maxiters=10
|
|
|
|
modelrxte="allsky/modelrxte_ait_3to20keV_flux_2deg.fits"
|
|
hdulist = fits.open(datadir+modelrxte)
|
|
w = wcs.WCS(hdulist[0].header)
|
|
smap =hdulist[0].data
|
|
|
|
sx=int(hdulist[0].header['NAXIS1'])
|
|
sy=int(hdulist[0].header['NAXIS2'])
|
|
|
|
# fill AITOF map indexes
|
|
ds9x=[]
|
|
ds9y=[]
|
|
for i,row in df.iterrows():
|
|
lon=row['LON']
|
|
lat=row['LAT']
|
|
world = SkyCoord(lon,lat, frame=Galactic, unit="deg")
|
|
ra=world.fk5.ra.deg
|
|
dec=world.fk5.dec.deg
|
|
pixcrd = w.wcs_world2pix([(lon,lat)], 1)
|
|
y=int(pixcrd[0][0])
|
|
x=int(pixcrd[0][1])
|
|
ds9x.append(x)
|
|
ds9y.append(y)
|
|
#print(x,y,smap[y-1,x-1])
|
|
df['DS9Y']=ds9x
|
|
df['DS9X']=ds9y
|
|
|
|
mean_map = np.array([[0.0 for i in range(sx)] for j in range(sy)])
|
|
sign_map = np.array([[0.0 for i in range(sx)] for j in range(sy)])
|
|
sem_map = np.array([[0.0 for i in range(sx)] for j in range(sy)])
|
|
cnt_map = np.array([[0 for i in range(sx)] for j in range(sy)])
|
|
|
|
obsid_map = {}
|
|
grxe_map = {}
|
|
grxe_err_map = {}
|
|
|
|
for i in range(sx):
|
|
for j in range(sy):
|
|
world = w.wcs_pix2world([(i+1,j+1)], 1)
|
|
lon = world[0][0]
|
|
lat = world[0][1]
|
|
if(np.isnan(lon) or np.isnan(lat)):
|
|
continue
|
|
|
|
|
|
|
|
ds9i=i+1
|
|
ds9j=j+1
|
|
df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j))
|
|
if (len(df0) <= nscw_min):
|
|
continue
|
|
|
|
# check coordinates
|
|
#print("***",i+1,j+1,lon,lat,smap[j][i])
|
|
#for i0,row0 in df0.iterrows():
|
|
# print(row0['LON'],row0['LAT'],row0['GRXE'])
|
|
grxe_err = np.array(df0['GRXE_ERR'])
|
|
perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False)
|
|
print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(key,enkey,grxe_err_cut,perc))
|
|
|
|
df0 = df.query('DS9X == {} & DS9Y == {} & GRXE_ERR < {}'.format(ds9i,ds9j,perc))
|
|
if (len(df0) <= nscw_min):
|
|
continue
|
|
|
|
obsid = np.array(df0['OBSID'])
|
|
grxe = np.array(df0['GRXE'])
|
|
grxe_err = np.array(df0['GRXE_ERR'])
|
|
|
|
sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=maxiters)
|
|
filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=maxiters, return_bounds=True)
|
|
filtered_grxe = filtered_data[0]
|
|
filtered_min = filtered_data[1]
|
|
filtered_max = filtered_data[2]
|
|
|
|
# final error on flux measurement ~RMS/sqrt(n)
|
|
sg_sem = sem(grxe[filtered_grxe.mask==False])
|
|
obsid=obsid[filtered_grxe.mask==False]
|
|
|
|
#print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem))
|
|
|
|
#plt.hist(grxe, bins=n_bins, range=minmax)
|
|
#plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=minmax)
|
|
#plt.show()
|
|
|
|
mean_map[j][i] = sg_mean
|
|
sem_map[j][i] = sg_sem
|
|
sign_map[j][i] = sg_mean/sg_sem
|
|
cnt_map[j][i] = df0.shape[0]
|
|
|
|
dkey="{:04d}{:04d}".format(j,i)
|
|
obsid_map[dkey] = obsid
|
|
grxe_map[dkey] = grxe
|
|
grxe_err_map[dkey] = grxe_err
|
|
|
|
|
|
|
|
""" Filter by error map """
|
|
# Calculate the percentiles across the x and y dimension
|
|
perc = np.percentile(sem_map, sem_cut, axis=(0, 1), keepdims=False)
|
|
|
|
print("{} {}: {}% cut of SEM map: {:.2f} mCrab".format(key,enkey,sem_cut,perc))
|
|
|
|
idx=np.where(sem_map > perc)
|
|
mean_map[idx]=0.0
|
|
sem_map[idx]=0.0
|
|
cnt_map[idx]=0
|
|
sign_map[idx]=0.0
|
|
|
|
if not os.path.exists(mapsdir):
|
|
os.makedirs(mapsdir)
|
|
|
|
hdulist[0].data=mean_map
|
|
hdulist.writeto(mapsdir+fn.replace(".fits",".mean.fits"),overwrite=True)
|
|
|
|
hdulist[0].data=sem_map
|
|
hdulist.writeto(mapsdir+fn.replace(".fits",".error.fits"),overwrite=True)
|
|
|
|
hdulist[0].data=cnt_map
|
|
hdulist.writeto(mapsdir+fn.replace(".fits",".cnt.fits"),overwrite=True)
|
|
|
|
hdulist[0].data=sign_map
|
|
hdulist.writeto(mapsdir+fn.replace(".fits",".sign.fits"),overwrite=True)
|
|
|
|
sys.exit()
|
|
print("Prepare data for fine map")
|
|
|
|
obsid_list=[]
|
|
grxe_list=[]
|
|
grxe_err_list=[]
|
|
for i in range(sx):
|
|
for j in range(sy):
|
|
""" Use mask """
|
|
if not (cnt_map[j][i]>0):
|
|
continue
|
|
world = w.wcs_pix2world([(i+1,j+1)], 1)
|
|
lon = world[0][0]
|
|
lat = world[0][1]
|
|
if(np.isnan(lon) or np.isnan(lat)):
|
|
continue
|
|
ds9i=i+1
|
|
ds9j=j+1
|
|
df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j))
|
|
if (len(df0) <= nscw_min):
|
|
continue
|
|
|
|
dkey="{:04d}{:04d}".format(j,i)
|
|
for scw in obsid_map[dkey]:
|
|
obsid_list.append(scw.decode("UTF-8"))
|
|
for grxe in grxe_map[dkey]:
|
|
grxe_list.append(grxe)
|
|
for grxe in grxe_err_map[dkey]:
|
|
grxe_err_list.append(grxe)
|
|
|
|
|
|
coldefs = fits.ColDefs([
|
|
fits.Column(name='OBSID', format='11A', array=obsid_list),
|
|
])
|
|
|
|
fout = fn.replace(".fits",".grxe.fits")
|
|
hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE')
|
|
hdu.header['MISSION'] = ('INTEGRAL', '')
|
|
#hdu.header['TELESCOP'] = (outkey, '')
|
|
hdu.header['INSTITUT'] = ('IKI', 'Affiliation')
|
|
hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person')
|
|
hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail')
|
|
#hdu.add_checksum()
|
|
print(hdu.columns)
|
|
hdu.writeto(proddir+fout, overwrite=True)
|
|
|
|
with fits.open(proddir+fout, mode='update') as hdus:
|
|
hdus[1].add_checksum()
|