maps

scripts/03_grxe_map.py

@@ -0,0 +1,155 @@
+#!/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="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)])
+
+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 = np.array(df0['GRXE'])
+
+        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])
+        
+        #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]
+
+
+# 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(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)