In [1]:
from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
import aplpy
from astropy.modeling import models, fitting
from astropy import wcs
import wcsaxes
%matplotlib inline
#import ROOT
from header_utils import verify_header
import aplpy.regions
import pyregion
#from astropy import pyregion
from os import path, listdir
import warnings
In [3]:
#warnings.filterwarnings('module')
home = path.expanduser("~")
gc_dir = home + "/Dropbox/GalacticCenter"
fitsdir = gc_dir +"/fits/"
filename = fitsdir + "/stage6_wobble_disp5t_skymap_subtraction_center_SgrA_maps.fits"
update_file = fits.open(filename, mode='update')
#update_header = fits.getheader(filename)
data_up, update_header = fits.getdata(filename, header=True)
print(type(update_header))
#verify_header(update_header) # could not find RADECSYS
print(update_header)
update_file.close()
#fits.writeto(filename, data_up, update_header, clobber=False)
In [4]:
data_fixed, header_fixed = fits.getdata(fitsdir + "release_galactic_skymap_fixed_head.fits", header=True)
#print(header_fixed)
In [5]:
fig = aplpy.FITSFigure(fitsdir + "release_skymap.fits")
fig.show_colorscale(vmin=10, vmax=35, cmap='hot')
fig.show_colorbar()
fig.set_tick_labels_xformat('dd.dd')
fig.set_tick_labels_yformat('dd.dd')
#fig.recenter(0,0,width=1,height=1)
fig.recenter(266.415, -29.006, width=0.4, height=0.4)
#fig.show_circles(266.5, -29, 0.1)
#fig.add_label(266.5, -29, "test", color='white', fontsize='24')
In [6]:
fil = fitsdir + "release_galactic_skymap.fits"
wcsGC = wcs.WCS(fil)
data = fits.getdata(fil)
data[120,120]
Out[6]:
In [7]:
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(1, 1, 1, projection = wcsGC)
plt.imshow(data, origin='lower')
ax.set_xlabel("Galactic Longitude")
ax.set_ylabel("Galactic Latitude")
#ax.cla()
#ax.set_xlim(110,128)
#ax.set_ylim(112,127)
#fig.add_subplot?
In [8]:
HDUlist = fits.open(fitsdir + "release_skymap.fits")
HDUlist.info()
hdu_2 = HDUlist[2]
excess_data_x = HDUlist[2].data
array1 = excess_data_x.copy()
print(array1)
print(array1[120,120])
print(type(HDUlist))
print(type(HDUlist[2]))
In [9]:
plt.imshow(excess_data_x)
#matplotlib.pyplot.colorbar()
Out[9]:
In [10]:
print("Max: ",np.std(excess_data_x))
In [11]:
print(models.Gaussian2D)
symGauss = models.Gaussian2D(amplitude=30,x_mean=266,y_mean=-28,x_stddev=.17,y_stddev=.17)
In [12]:
pinit2 = models.Gaussian2D(x_mean=5, y_mean=5, x_stddev=2, y_stddev=2, theta=0.)
pinit = models.Gaussian2D(amplitude=30, x_mean=2.5, y_mean=2.5, x_stddev=0.1, y_stddev=0.1, theta=0.)
pinit3 = pinit + pinit2
fit_p = fitting.LevMarLSQFitter()
In [13]:
#fit_p?
#fit_p.supported_constraints?
#newmodel = fit_p(pinit, excess_data)
In [14]:
pinit.theta.fixed = True
pinit.amplitude.min = 0.
pinit.bounds
Out[14]:
In [15]:
data, header = fits.getdata(fitsdir+"release_skymap.fits", header=True)
np.nanmax(data)
np.where(data == np.nanmax(data))
data2 = np.nan_to_num(data)
In [16]:
plt.imshow(data2)
Out[16]:
In [17]:
#np.floor(fig.world2pixel(266.5, -29))
#w.wcs_pix2world(100,100)
#w.printwcs()
In [18]:
figure_SgrA_excess = aplpy.FITSFigure(hdu_2)
figure_SgrA_excess.show_colorscale(cmap='hot')
figure_SgrA_excess.show_colorbar()
#aplpy.FITSFigure?
HDUlist.close()
In [19]:
skymap_data, skymap_header = fits.getdata(fitsdir+"simpleFITStestRalph.fits", header = True)
excess_data_N = np.nan_to_num(skymap_data)
# Sgr A* middle bins
print(excess_data_N[119:123,119:122])
print(skymap_header)
# G0.9+0.1
print(skymap_data[148:151,159:161])
In [20]:
def setup_sig_fits(fig):
"""set standard options for sig map"""
fig.show_colorscale(smooth=1, interpolation='spline16')
#plasma
fig.show_colorbar()
cbar = fig_apl.colorbar
cbar.set_axis_label_text('Significance')
fig_apl.tick_labels.set_xformat('dd.d')
fig_apl.tick_labels.set_yformat('dd.d')
# end set_sig_fits
In [21]:
sig_filename_fits = fitsdir+"stage6_wobble_disp5t_skymap_subtraction_center_SgrA_maps.fits"
fig_apl = aplpy.FITSFigure(sig_filename_fits)
fig_apl.set_title("J1746-285 Significance Skymap")
setup_sig_fits(fig_apl)
fig_apl.show_regions(fitsdir+"radio_arc_exclusions.reg")
#fig_apl.show_regions(fitsdir+"stage6_wobble_disp5t_skymap_subtraction_center_SgrA_maps_ds9_full.reg")
fig_apl.recenter(0.0625, -0.1, width=0.65, height=0.5) # 1746
#fig_apl.recenter(-0.25, -0.25, width=3.5, height=1.5) # center
#fig_apl.recenter(0.7, -0.05, width=0.5, height=0.5) # B2
#plt.savefig(gc_dir+"/plots/skymaps/J1746-285_radioArc_subtracted.png")
#fig_apl.save('SgrB2_sigmap.pdf')
#fig_apl.show_contour(data)
In [26]:
sig_filename_fits = fitsdir + "SgrA_disp5t_excludeMore_skymap_4tels_s6.fits"
fig_apl = aplpy.FITSFigure(sig_filename_fits)
setup_sig_fits(fig_apl)
fig_apl.show_colorscale(cmap='hot', vmax=10., smooth=1, interpolation='spline16')
fig_apl.recenter(359.95, -0.05, width=0.75, height=0.75)
fig_apl.set_title("SgrA* Significance")
fig_apl.show_regions(fitsdir+"SgrA_position_uncertainty.reg")
fig_apl.show_markers(359.94, -0.04, marker='.') # radio position
#, fillstyles='full')
plt.savefig(home+"/Downloads/GC_SgrA_closeup.png")
In [23]:
# look at G1.9
warnings.filterwarnings('ignore')
sig_filename_fits = fitsdir + "SgrA_disp5t_excludeMore_skymap_4tels_s6.fits"
fig_apl = aplpy.FITSFigure(sig_filename_fits)
setup_sig_fits(fig_apl)
fig_apl.show_colorscale(vmax=10)
#fig_apl.recenter(359.95, -0.05, width=3., height=3.)
#fig_apl.set_title("Galactic Center Significance")
#fig_apl.show_regions(fitsdir+"G1.9+0.3_and_wobbles_ds9.reg")
#plt.savefig(gc_dir+"/plots/skymaps/GC_wobbles_G1.9+0.3.png")
#fig_apl.show_circles(100, 100, 25)
In [ ]:
warnings.filterwarnings('ignore')
#sig_filename_fits = fitsdir+"release_galactic_skymap_fixed_head.fits"
sig_filename_fits = fitsdir + "SgrA_disp5t_excludeMore_skymap_4tels_s6.fits"
fig_apl = aplpy.FITSFigure(sig_filename_fits)
fig_apl.set_title("Galactic Center Annulus Region")
setup_sig_fits(fig_apl)
#fig_apl.show_colorscale(smooth=1)
#cbar = fig_apl.show_colorbar()
#fig_apl.show_regions(fitsdir+"GC_annulus_region.reg")
fig_apl.recenter(359.95, -0.05, width=1.0, height=1.0) # center
In [ ]:
fig_apl = aplpy.FITSFigure(fitsdir+"/SgrA_disp5t_E5000_4tels_skymap.fits")
setup_sig_fits(fig_apl)
fig_apl.set_title("Significance skymap (E > 5TeV)")
#fig_apl.show_colorscale(vmax=7.)
fig_apl.recenter(359.7, -0.05, width=3., height=1.5)
fig_apl.show_regions(fitsdir+"/stage6_wobble_disp5t_skymap_subtraction_center_SgrA_maps_ds9_full.reg")
#plt.savefig(gc_dir+"/plots/skymaps/SgrA_sigmap_E5000_fits.png")
In [ ]:
filebase = home+"/Documents/GC_material/NRO_radio_fits_data/GC.CS10.10.FITS.files/A10MAP"
data, header = fits.getdata(filebase+".1.FITS", header=True)
for i in range(2, 41):
data = data + fits.getdata(filebase+'.'+str(i)+".FITS")
filename = filebase + "_sum.fits"
#fits.writeto(filename, data, header, clobber=False)
fig_apl = aplpy.FITSFigure(filename)
setup_sig_fits(fig_apl)
fig_apl.show_contour(data)
In [ ]:
char = 'V'
filedir = home + "/Documents/GC_material/NRO_radio_fits_data/GC.CS10."+char+"V.FITS.files/"
filelist = listdir(filedir)
#data, header = fits.getdata(filedir+filelist[0], header=True)
#for file in filelist[1:]:
# data = data + fits.getdata(filedir+file)
filename = filedir + char + "VMAP_sum.fits"
#fits.writeto(filename, data, header, clobber=False)
fig_apl = aplpy.FITSFigure(filename)
#setup_sig_fits(fig_apl)
fig_apl.show_colorscale(smooth=1, interpolation='spline16')
fig_apl.show_colorbar()
#fig_apl.show_contour(data)
In [ ]:
fig_apl.show_colorscale(smooth=1, interpolation='spline16')
In [ ]:
fig.show_colorscale(smooth=1, interpolation='spline16')
#plasma
fig.show_colorbar()
cbar = fig_apl.colorbar
cbar.set_axis_label_text('Significance')
fig_apl.tick_labels.set_xformat('dd.d')
fig_apl.tick_labels.set_yformat('dd.d')
In [ ]:
fig_apl = aplpy.FITSFigure(home+"/Documents/GC_material/NRO_radio_fits_data/GCCS.fits", slices=[87])
setup_sig_fits(fig_apl)