IKI

scripts/03_grxe_galprof.py

@@ -0,0 +1,190 @@
+#!/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]
+lon_nbin = int(sys.argv[2])
+key="ALL"
+fn='detcnts.{}.{}.resid.fits'.format(enkey,key)
+
+print("Reading {}".format(proddir+fn))
+dat = Table.read(proddir+fn, unit_parse_strict='silent')
+df = dat.to_pandas()
+
+#df = df.query('abs(LAT) < {} & abs(LON) < {}'.format(5,5))
+
+n_bins = 80
+minmax=[-300,800]
+sigma=3
+maxiters=10
+
+with open(ignored_rev_file, 'rb') as fp:
+    ignored_rev = pickle.load(fp)
+    print(ignored_rev)
+    print("{} orbits ignored".format(len(ignored_rev)))
+
+ign=ignored_rev.tolist()
+
+
+glon, step = np.linspace(-lon_max, lon_max, num=lon_nbin, endpoint=False,retstep=True)
+print(glon,step)
+
+
+#
+# initiate 2d arrays
+#
+mean_map = np.array([0.0 for j in range(lon_nbin)])
+sign_map = np.array([0.0 for j in range(lon_nbin)])
+sem_map = np.array([0.0 for j in range(lon_nbin)])
+cnt_map = np.array([0 for j in range(lon_nbin)])
+
+mean_sim_map = np.array([0.0 for j in range(lon_nbin)])
+error_sim_map = np.array([0.0 for j in range(lon_nbin)])
+
+mean_sim={}
+for i in range(lon_nbin):
+    dkey="{:04d}".format(i)
+    mean_sim[dkey] = []
+
+
+obsid_map = {}
+grxe_map = {}
+grxe_err_map = {}
+
+# redefine simfrac for low number of ScWs in pixel
+#simfrac=2
+#nsim=10
+
+lon_arr=[]
+lat_arr=[]
+
+for i in range(lon_nbin):
+    dkey="{:04d}".format(i)
+
+    
+    df0 = df.query('LON > {} & LON <= {} & LAT > {} & LAT < {} & REV != @ign'.format(glon[i],glon[i]+step,-lat_max,lat_max))
+
+    if  (df0.shape[0] < nscw_min):
+        continue
+
+    sg_mean,sg_sem = get_spec(df0, sigma=sigma,  grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
+        
+    nsel = int(df0.shape[0]/simfrac)
+    print("lon {:.2f} ".format(glon[i]),"nsel=",nsel,df0.shape[0])
+    for n in range(nsim):
+        df1=df0.sample(nsel)
+        sg_mean1,sg_sem1 = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
+        mean_sim[dkey].append(sg_mean1)   
+
+        
+    #print('sg_sem',sg_sem)
+    mean_map[i] = sg_mean
+    sem_map[i]  = sg_sem
+    cnt_map[i]  = df0.shape[0]
+
+    for index, row in df0.iterrows():
+        lon_arr.append(row['LON'])
+        lat_arr.append(row['LAT'])
+
+
+
+        
+if not os.path.exists(profdir):
+    os.makedirs(profdir)
+
+print("saving simulations")
+for i in range(lon_nbin):
+    dkey="{:04d}".format(i)
+    data=mean_sim[dkey]
+    #print("{} size {}".format(dkey,len(data)))
+    #if(len(data)>10):
+    (mu, sg) = norm.fit(data)
+    mean_sim_map[i] = mu
+    error_sim_map[i] = sg
+    
+
+coldefs = fits.ColDefs([
+    fits.Column(name='LON1', format='D', array=glon), 
+    fits.Column(name='LON2', format='D', array=glon+step), 
+    fits.Column(name='GRXE_FLUX', format='D', array=mean_map), 
+    fits.Column(name='GRXE_ERROR', format='D', array=sem_map), 
+    fits.Column(name='GRXE_SIM_FLUX', format='D', array=mean_sim_map), 
+    fits.Column(name='GRXE_SIM_ERROR', format='D', array=error_sim_map), 
+    fits.Column(name='nScW', format='J', array=cnt_map), 
+])
+
+fout = fn.replace(".fits",".galprof.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(profdir+fout, overwrite=True)
+
+with fits.open(profdir+fout, mode='update') as hdus:
+    hdus[1].add_checksum() 
+
+    
+# to check, which ScWs were used
+coldefs = fits.ColDefs([
+    fits.Column(name='LON', format='D', array=lon_arr), 
+    fits.Column(name='LAT', format='D', array=lat_arr), 
+])
+
+
+fout = fn.replace(".fits",".gal.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(profdir+fout, overwrite=True)
+
+with fits.open(profdir+fout, mode='update') as hdus:
+    hdus[1].add_checksum() 
+
+

scripts/03_grxe_galprof.sh

@@ -0,0 +1,4 @@
+./03_grxe_galprof.py E01 71
+./03_grxe_galprof.py E13 17
+./03_grxe_galprof.py E14 21
+

scripts/03_grxe_galprof_plot.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+
+__author__      = "Roman Krivonos"
+__copyright__   = "Space Research Institute (IKI)"
+
+from astropy.table import Table, Column
+
+import matplotlib.pyplot as plt
+import sys
+
+from ridge.utils import *
+from ridge.config import *
+
+
+scale = 100.0
+
+fn="detcnts.E01.ALL.resid.galprof.fits"
+dat = Table.read(profdir+fn, unit_parse_strict='silent')
+df1 = dat.to_pandas().sort_values(by=['LON1'])
+
+fn="detcnts.E14.ALL.resid.galprof.fits"
+dat = Table.read(profdir+fn, unit_parse_strict='silent')
+df2 = dat.to_pandas().sort_values(by=['LON1'])
+
+fn="detcnts.E13.ALL.resid.galprof.fits"
+dat = Table.read(profdir+fn, unit_parse_strict='silent')
+df3 = dat.to_pandas().sort_values(by=['LON1'])
+
+
+fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, figsize=(9, 7), dpi=100)
+#fig.suptitle('Vertically stacked subplots')
+#plt.figure(figsize=(8, 6), dpi=80)
+
+for axis in ['top','bottom','left','right']:
+    ax1.spines[axis].set_linewidth(1)
+    ax2.spines[axis].set_linewidth(1)
+    ax3.spines[axis].set_linewidth(1)
+ax1.tick_params(axis="both", width=1, labelsize=14)
+ax2.tick_params(axis="both", width=1, labelsize=14)
+ax3.tick_params(axis="both", width=1, labelsize=14)
+
+ax1.set_title("25-60 keV")
+ax2.set_title("60-80 keV")
+ax3.set_title("80-200 keV")
+
+ax1.set_ylim(-0.3,1.8)
+ax2.set_ylim(-0.3,1.8)
+ax3.set_ylim(-0.3,1.8)
+
+ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ax3.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+lon=(df1['LON1']+df1['LON2'])/2
+ax1.errorbar(lon, df1['GRXE_SIM_FLUX']/scale,
+             yerr=df1['GRXE_SIM_ERROR']/scale,
+             xerr=(df1['LON2']-df1['LON1'])/2,
+             fmt='o' )
+ax1.grid(visible=True)
+
+lon=(df2['LON1']+df2['LON2'])/2
+ax2.errorbar(lon, df2['GRXE_SIM_FLUX']/scale,
+             yerr=df2['GRXE_SIM_ERROR']/scale,
+             xerr=(df2['LON2']-df2['LON1'])/2,fmt='o' )
+ax2.grid(visible=True)
+
+lon=(df3['LON1']+df3['LON2'])/2
+ax3.errorbar(lon, df3['GRXE_SIM_FLUX']/scale,
+             yerr=df3['GRXE_SIM_ERROR']/scale,
+             xerr=(df3['LON2']-df3['LON1'])/2, fmt='o' )
+ax3.grid(visible=True)
+
+#x = np.arange(-10, 10, 0.001)
+#plot normal distribution with mean 0 and standard deviation 1
+#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
+	    
+plt.xlabel('Galactic Longitude',fontsize=14, fontweight='normal')
+ax2.set_ylabel('GRXE flux, x100 mCrab',fontsize=14, fontweight='normal')
+
+#plt.xscale('linear')
+#plt.yscale('linear')
+
+plt.savefig(profdir+'galprof.png', bbox_inches='tight')
+plt.close(fig)

scripts/03_grxe_map.py

@@ -72,6 +72,8 @@ sx=int(hdulist[0].header['NAXIS1'])
 sy=int(hdulist[0].header['NAXIS2'])
 
 # fill AITOF map indexes
+# Already done in 02_grxe_resid.py
+"""
 ds9x=[]
 ds9y=[]
 for i,row in df.iterrows():
@@ -88,6 +90,7 @@ for i,row in df.iterrows():
     #print(x,y,smap[y-1,x-1])
 df['DS9Y']=ds9x
 df['DS9X']=ds9y
+"""
 
 #
 # initiate 2d arrays
@@ -128,28 +131,31 @@ for i in range(sx):
         
         ds9i=i+1
         ds9j=j+1
+        #df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j))
         df0 = df.query('DS9X == {} & DS9Y == {} & REV != @ign'.format(ds9i,ds9j))
-        if (df0.shape[0] <= nscw_min):
+        if  (df0.shape[0] < nscw_min):
             continue
-        print("*** *** SUM *** ***")
-        print(np.sum(df0["GRXE"]))
+
+        #print("*** *** REV *** ***")
+        #print(df0["REV"])
+        
         # 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'])
         
 
-        sg_mean,sg_sem = get_spec(df0, sigma=sigma,  grxe_err_cut=grxe_err_cut, enkey=enkey)
+        sg_mean,sg_sem = get_spec(df0, sigma=sigma,  grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
         
         nsel = int(df0.shape[0]/simfrac)
-        print("nsel=",nsel,df0.shape[0],len(df0['GRXE']))
+        #print("nsel=",nsel,df0.shape[0])
         for n in range(nsim):
             df1=df0.sample(nsel)
-            sg_mean1,sg_sem1 = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey)
+            sg_mean1,sg_sem1 = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
             mean_sim[dkey].append(sg_mean1)   
 
         
-        print('sg_sem',sg_sem)
+        #print('sg_sem',sg_sem)
         mean_map[j][i] = sg_mean
         sem_map[j][i]  = sg_sem
         sign_map[j][i] = sg_mean/sg_sem
@@ -166,16 +172,22 @@ for i in range(sx):
         
 """ 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)
 
+
+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)
+print("index size {}".format(len(idx)))
 mean_map[idx]=0.0
 sem_map[idx]=0.0
 cnt_map[idx]=0
 sign_map[idx]=0.0
 
+#mean_sim_map[idx]=0.0
+#error_sim_map[idx]=0.0
+
+
 if not os.path.exists(mapsdir):
     os.makedirs(mapsdir)
 
@@ -191,15 +203,24 @@ 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)
 
-
+print("saving simulations")
 for i in range(sx):
     for j in range(sy):
         dkey="{:04d}{:04d}".format(j,i)
         data=mean_sim[dkey]
-        (mu, sg) = norm.fit(data)
-        mean_sim_map[j][i] = mu
-        error_sim_map[j][i] = sg
+        #print("{} size {}".format(dkey,len(data)))
+        if(len(data)>10):
+            (mu, sg) = norm.fit(data)
+            mean_sim_map[j][i] = mu
+            error_sim_map[j][i] = sg
     
+perc = np.percentile(error_sim_map, sem_cut, axis=(0, 1), keepdims=False)
+print("{} {}: {}% cut of SEM map: {:.2f} mCrab".format(key,enkey,sem_cut,perc))
+
+idx=np.where(error_sim_map > perc)
+print("index size {}".format(len(idx)))
+mean_sim_map[idx]=0.0
+error_sim_map[idx]=0.0
 
 hdulist[0].data=mean_sim_map
 hdulist.writeto(mapsdir+fn.replace(".fits",".sim.mean.fits"),overwrite=True)

scripts/03_grxe_map.sh

@@ -1,14 +1,15 @@
-./03_grxe_map.py E02
-./03_grxe_map.py E03
-./03_grxe_map.py E04
-./03_grxe_map.py E05
-./03_grxe_map.py E06
-./03_grxe_map.py E07
-./03_grxe_map.py E08
-./03_grxe_map.py E09
-./03_grxe_map.py E10
-./03_grxe_map.py E11
-./03_grxe_map.py E12
+./03_grxe_map.py E01
+#./03_grxe_map.py E02
+#./03_grxe_map.py E03
+#./03_grxe_map.py E04
+#./03_grxe_map.py E05
+#./03_grxe_map.py E06
+#./03_grxe_map.py E07
+#./03_grxe_map.py E08
+#./03_grxe_map.py E09
+#./03_grxe_map.py E10
+./03_grxe_map.py E13
+./03_grxe_map.py E14