# read in data and prepare for plotting
hdu = fits.open('../atomic/FIRAS_LINE_EMISSION_MAP_HIGH.fits')
# get the NII line flux at each (l,b)
# LINFRQ9 = 2459.4 / [N II] line frequency, in GHz
# figured the data axis out by looking at the headers and playing around...
data = hdu[1].data
l = data['GAL_LON']
b = data['GAL_LAT']
f = data['LINE_FLU'][:,3]
f2 = data['LINE_FL2'][:,3]
f3 = data['LINE_FL3'][:,3]
# only use the regions with lower noise
good = f2 < 2
l = l[good]
b = b[good]
f = f[good]
#print(f.min(), f.max())
l_flip = l > 180
l[l_flip] = l[l_flip] - 360
# make into a regular grid
x, y = np.mgrid[-180:180:0.3, -60:60:0.3]
im = griddata((l,b), f, (x,y), method='cubic', fill_value=-100)
im = np.rot90(im, k=3)