Example Gallery
import numpy as np
import matplotlib.pyplot as plt
from urllib.request import urlretrieve
from plotastrodata.analysis_utils import AstroData, AstroFrame
from plotastrodata.plot_utils import PlotAstroData as pad
data_url = ('https://raw.githubusercontent.com/'
'yusukeaso-astron/plotastrodata/main/example_data')
files = ['test2D.fits',
'test2D_2.fits',
'test2Damp.fits',
'test2Dang.fits',
'test3D.fits',
'testPV.fits']
for s in files:
urlretrieve(f'{data_url}/{s}', s)
Matplotlib is building the font cache; this may take a moment.
2D image
d = AstroData(fitsimage='test2D.fits', Tb=True, sigma=5e-3)
f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')
f.read(d)
p = pad(rmax=0.8, center='ICRS 04h04m43.07s 26d18m56.20s')
p.add_color(**d.todict(), cblabel='Tb (K)')
p.add_contour(fitsimage='test2D_2.fits', colors='r', sigma=5e-3)
p.add_contour(fitsimage='test2D.fits', xskip=2, yskip=2, sigma=5e-3)
p.add_segment(ampfits='test2Damp.fits',
angfits='test2Dang.fits',
xskip=3, yskip=3, sigma=None)
p.add_scalebar(length=50 / 140, label='50 au')
p.add_text([0.3, 0.3], slist='text')
p.add_marker('04h04m43.07s 26d18m56.20s')
p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5])
p.add_arrow([0.4, 0.4], anglelist=150, rlist=0.5)
p.add_region('ellipse', [0.2, 0.8], majlist=0.4, minlist=0.2, palist=45)
p.set_axis_radec(nticksminor=5, title={'label': '2D image', 'loc': 'right'})
p.savefig('test2D.png', show=True)
3D channel maps
p = pad(rmax=0.8, fitsimage='test3D.fits', vmin=-5, vmax=5, vskip=2)
p.add_color(fitsimage='test3D.fits',
stretch='log', sigma='hist,edge')
p.add_contour(fitsimage='test3D.fits',
colors='r', sigma='hist,edge')
p.add_contour(fitsimage='test2D.fits', colors='b', sigma=5e-3)
p.add_segment(ampfits='test2Damp.fits',
angfits='test2Dang.fits',
xskip=3, yskip=3, sigma=None)
p.add_scalebar(length=50 / 140, label='50 au')
p.add_text([[0.3, 0.3]], slist=['text'], include_chan=[0, 1, 2])
p.add_marker([0.7, 0.7])
p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5],
include_chan=[6, 7, 8])
p.add_arrow([[0.4, 0.4]], anglelist=[150], rlist=[0.5],
include_chan=[9, 10, 11])
p.add_region('rectangle', [[0.2, 0.8]],
majlist=[0.4], minlist=[0.2], palist=[45],
include_chan=[12, 13, 14])
p.set_axis(grid={}, title='3D channel maps')
p.savefig('test3D.png', show=True)
PV diagram
d = AstroData(fitsimage='testPV.fits', Tb=True, sigma=None, pvpa=60)
f = AstroFrame(pv=True)
f.read(d)
p = pad(rmax=0.8, pv=True, swapxy=True, vmin=-5, vmax=5, figsize=(6, 7))
p.add_color(**d.todict(), cblabel='Tb (K)', cblocation='top')
p.add_contour(fitsimage='testPV.fits', sigma=1e-3, pvpa=60,
colors='r')
p.add_text([0.3, 0.3], slist='text')
p.add_marker([[0.5, 0.5]])
p.set_axis(title='PV diagram')
p.savefig('testPV.png', show=True)
log log PV
p = pad(rmax=0.8 * 140, pv=True, quadrants='13', vmin=-5, vmax=5, dist=140)
p.add_color(fitsimage='testPV.fits', Tb=True, sigma=None,
cblabel='Tb (K)', show_beam=False, pvpa=60)
p.add_contour(fitsimage='testPV.fits', colors='r',
sigma=1e-3, show_beam=False, pvpa=60)
p.set_axis(title='loglog PV diagram', loglog=20)
p.savefig('testloglogPV.png', show=True)
RGB
d = AstroData(fitsimage='test3D.fits', Tb=True, sigma=5e-3)
f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')
f.read(d)
dblue = np.sum(d.data[0:20], axis=0) * d.dv
dgreen = np.sum(d.data[20:41], axis=0) * d.dv
dred = np.sum(d.data[41:61], axis=0) * d.dv
d.data = [dred, dgreen, dblue]
d.sigma = [d.sigma * d.dv * np.sqrt(20)] * 3
p = pad(rmax=0.8, center='ICRS 04h04m43.07s 26d18m56.20s')
p.add_rgb(**d.todict())
p.add_contour(fitsimage='test2D_2.fits', colors='r', sigma=5e-3)
p.add_contour(fitsimage='test2D.fits', xskip=2, yskip=2, sigma=5e-3)
p.add_segment(ampfits='test2Damp.fits',
angfits='test2Dang.fits',
xskip=3, yskip=3, sigma=None)
p.add_scalebar(length=50 / 140, label='50 au')
p.add_text([0.3, 0.3], slist='text')
p.add_marker('04h04m43.07s 26d18m56.20s')
p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5])
p.add_arrow([0.4, 0.4], anglelist=150, rlist=0.5)
p.add_region('ellipse', [0.2, 0.8], majlist=0.4, minlist=0.2, palist=45)
p.set_axis_radec(nticksminor=5, title={'label': '2D RGB', 'loc': 'right'})
p.savefig('test2Drgb.png', show=True)
Line profile
from plotastrodata.plot_utils import plotprofile
plotprofile(fitsimage='test3D.fits', ellipse=[[0.2, 0.2, 0]] * 2,
flux=True, sigma='hist,edge',
coords=['04h04m43.045s 26d18m55.766s',
'04h04m43.109s 26d18m56.704s'],
gaussfit=True, savefig='testprofile.png', show=True, width=2)
sigma has been divided by sqrt(2) because of binning in the v-axis.
Gauss (peak, center, FWHM): [0.02376936 1.71860519 2.02614765]
Gauss uncertainties: [0.00194454 0.08494428 0.15898626]
Estimated sigma: 0.003254768057067125
Gauss (peak, center, FWHM): [ 0.02583402 -1.7780023 1.70315156]
Gauss uncertainties: [0.00210303 0.0745902 0.14631028]
Estimated sigma: 0.0034813829219451893
Spatial slice
from plotastrodata.plot_utils import plotslice
plotslice(length=1.6, pa=270, fitsimage='test2D.fits',
center='04h04m43.07s 26d18m56.20s', sigma=5e-3,
savefig='testslice.png', show=True)
3D html figure (Plotly)
from plotastrodata.plot_utils import plot3d
d = AstroData(fitsimage='test3D.fits', Tb=True, sigma=5e-3)
f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')
f.read(d)
d.data = np.sum(d.data, axis=0) * d.dv
d.sigma = d.sigma * d.dv * np.sqrt(len(d.v))
vplus = {'data': d.data, 'sigma': d.sigma, 'cmap': 'gray'}
plot3d(rmax=0.8, vmin=-5, vmax=5, fitsimage='test3D.fits',
outname='./_static/plotly/test3D.html', levels=[3, 6, 9],
sigma='hist,edge', vplus=vplus, show=False)
Animation
import matplotlib.pyplot as plt
import matplotlib.animation as animation
nchans = 16
def update_plot(i):
p = pad(rmax=0.8, fitsimage='test3D.fits',
vmin=-5, vmax=5, vskip=4,
channelnumber=i, fig=fig)
p.add_color(fitsimage='test3D.fits',
stretch='log', sigma='hist,edge')
p.add_scalebar(length=50 / 140, label='50 au')
p.set_axis_radec(grid={}, title='3D channel maps')
p.fig.tight_layout()
fig = plt.figure(figsize=(7, 5)) # Giving figsize may help ffmpeg.
ani = animation.FuncAnimation(fig, update_plot, frames=nchans)
#Writer = animation.writers['ffmpeg'] # for mp4
Writer = animation.writers['pillow'] # for gif
writer = Writer(fps=1) # frame per second
#ani.save('./_static/video/test_animation.mp4', writer=writer, dpi=64)
ani.save('./_static/video/test_animation.gif', writer=writer, dpi=64)
plt.close()
Noise estimate
from plotastrodata.noise_utils import Noise
x = np.linspace(-1.5, 1.5, 301)
y = np.linspace(-1.5, 1.5, 301)
x, y = np.meshgrid(x, y)
r = np.hypot(x, y)
data = np.random.randn(*np.shape(r))
data[r > 1.5] = np.nan
data = data * np.exp2(r**2)
n = Noise(data=data, sigma='hist-pbcor')
n.gen_histogram()
n.fit_histogram()
print('Noise (mean, std, Rout):',
[n.mean, n.std, float(n.popt[2])])
n.plot_histogram(savefig='noise.png', show=True)
Noise (mean, std, Rout): [0.006263354870539588, 0.9992137678658358, 1.5017394455616297]
Line-of-sight velocity with 3D rotation
from plotastrodata.los_utils import sys2obs, polarvel2losvel
incl = 60
theta0 = 60
xscale = np.sin(np.radians(theta0))**2
vscale = 1 / np.sin(np.radians(theta0))
xsys = np.linspace(0, 3 / xscale, 200)
ysys = np.sqrt(2 * xsys + 1)
zsys = np.zeros_like(xsys)
r = np.hypot(xsys, ysys)
theta = np.zeros_like(r) + np.pi / 2
phi = np.arctan2(ysys, xsys)
v_r = -np.sqrt(1 / r * (2 - 1 / r))
v_theta = np.zeros_like(v_r)
v_phi = 1 / r
xsys = xsys * xscale
ysys = ysys * xscale
zsys = zsys * xscale
v_r = v_r * vscale
v_theta = v_theta * vscale
v_phi = v_phi * vscale
xlist, ylist, vlist = [], [], []
for phi0 in np.linspace(0, 360, 9):
xobs, yobs, zobs = sys2obs(xsys=xsys, ysys=ysys, zsys=zsys,
incl=incl, phi0=phi0, theta0=theta0)
vlos = polarvel2losvel(v_r=v_r, v_theta=v_theta, v_phi=v_phi,
theta=theta, phi=phi,
incl=incl, phi0=phi0, theta0=theta0)
xlist.append(xobs)
ylist.append(yobs)
vlist.append(vlos)
fig, ax = plt.subplots()
m = ax.scatter(xlist, ylist, c=vlist, cmap='coolwarm', vmin=-1.5, vmax=1.5)
fig.colorbar(m, ax=ax, label=r'$V_\mathrm{los} / \sqrt{GM_{*}/r_{c}}$')
ax.set_xlim(3, -3)
ax.set_ylim(-3, 3)
ax.set_aspect(1)
ax.set_xlabel(r'$x / r_{c}$')
ax.set_ylabel(r'$y / r_{c}$')
ax.grid()
fig.savefig('streamer.png')
plt.show()
FFT with a given center
from plotastrodata.fft_utils import fftcentering
x = np.linspace(-99.5, 99.5, 200)
f = np.where(np.abs(x)<10, 1, 0)
fig, ax = plt.subplots()
ax.plot(x, f)
ax.set_xlabel('x')
ax.set_ylabel('f')
fig.tight_layout()
fig.savefig('boxcar.png')
plt.show()
u = np.fft.fftshift(np.fft.fftfreq(len(x), d=x[1] - x[0]))
F = np.fft.fftshift(np.fft.fft(f))
fig, ax = plt.subplots()
ax.plot(u, np.real(F), label='real')
ax.plot(u, np.imag(F), label='imag')
ax.set_xlabel('u')
ax.set_ylabel('numpy.fft')
ax.legend()
fig.tight_layout()
fig.savefig('numpyfft.png')
plt.show()
F, u = fftcentering(f=f, x=x, xcenter=0)
fig, ax = plt.subplots()
ax.plot(u, np.real(F), label='real')
ax.plot(u, np.imag(F), label='imag')
ax.set_xlabel('u')
ax.set_ylabel('fftcentering')
ax.legend()
fig.tight_layout()
fig.savefig('fftcentering.png')
plt.show()
MCMC (emcee, corner, dynesty)
from plotastrodata.fitting_utils import EmceeCorner
from plotastrodata.plot_utils import set_rcparams
set_rcparams(fontsize=12)
def logl(p):
x1, x2, x3 = p
chi2 = (x1 / 1)**2 + (x2 / 2)**2 + (x3 / 4)**2
return -0.5 * chi2
fitter = EmceeCorner(bounds=[[-5, 5], [-10, 10], [-20, 20]],
logl=logl, progressbar=False, percent=[16, 84])
fitter.fit(nwalkersperdim=30, nsteps=11000, nburnin=1000,
#savechain='chain.npy'
)
print('best:', fitter.popt)
print('lower percentile:', fitter.plow)
print('50 percentile:', fitter.pmid)
print('higher percentile:', fitter.phigh)
fitter.getDNSevidence()
print('evidence:', fitter.evidence)
fitter.plotcorner(show=True, savefig='corner.png',
labels=['par1', 'par2', 'par3'],
cornerrange=[[-4, 4], [-8, 8], [-16, 16]])
fitter.plotchain(show=True, savefig='chain.png',
labels=['par1', 'par2', 'par3'],
ylim=[[-2, 2], [-4, 4], [-8, 8]])
best: [-0.0085724 -0.00404063 -0.02137008]
lower percentile: [-0.98518159 -1.97382879 -3.95667345]
50 percentile: [0.00519443 0.01654306 0.00702374]
higher percentile: [1.00380292 1.99999229 3.97098905]
evidence: 0.016730135126799137
Likelihood on a parameter grid
fitter.posteriorongrid(ngrid=[101, 201, 401])
print('best:', fitter.popt)
print('lower percentile:', fitter.plow)
print('50th percentile:', fitter.pmid)
print('higher percentile:', fitter.phigh)
print('evidence:', fitter.evidence)
fitter.plotongrid(show=True, savefig='grid.png',
labels=['par1', 'par2', 'par3'],
cornerrange=[[-4, 4], [-8, 8], [-16, 16]])
best: [0. 0. 0.]
lower percentile: [-1. -2. -4.]
50th percentile: [0. 0. 0.]
higher percentile: [1. 2. 4.]
evidence: 0.015477376407909737
Other functions
from plotastrodata.coord_utils import xy2coord, coord2xy
coord = xy2coord(xy=[[30, 90], [0, 0]], coordorg='00h00m00s 60d00m00s')
print(coord)
xy = coord2xy(coords=coord, coordorg='00h00m00s 60d00m00s')
print(np.round(xy, 2))
['03h16m25.58528421s +48d35m25.36040662s', '06h00m00s +00d00m00s']
[[30. 90.]
[ 0. -0.]]
import plotastrodata.const_utils as cu
print(cu.pc)
print(cu.M_sun)
print(cu.centi)
3.085677581491367e+16
1.988409870698051e+30
0.01
from plotastrodata.ext_utils import BnuT, JnuT
print(BnuT(T=30, nu=230e9))
print(JnuT(T=30, nu=230e9))
4.033705316844142e-16
24.818562479617647