{ "cells": [ { "cell_type": "markdown", "id": "d1ec450c-ba6e-47c8-90f2-d22b9a1710f4", "metadata": {}, "source": [ "# Example Gallery" ] }, { "cell_type": "code", "execution_count": null, "id": "466710c8-a924-4159-9cbe-2a05c7be4d70", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from urllib.request import urlretrieve\n", "\n", "from plotastrodata.analysis_utils import AstroData, AstroFrame\n", "from plotastrodata.plot_utils import PlotAstroData as pad\n", "\n", "\n", "data_url = ('https://raw.githubusercontent.com/'\n", " 'yusukeaso-astron/plotastrodata/main/example_data')\n", "\n", "files = ['test2D.fits',\n", " 'test2D_2.fits',\n", " 'test2Damp.fits',\n", " 'test2Dang.fits',\n", " 'test3D.fits',\n", " 'testPV.fits']\n", "\n", "for s in files:\n", " urlretrieve(f'{data_url}/{s}', s)" ] }, { "cell_type": "markdown", "id": "107f2297-852c-4d30-aab3-7bc6413f08f9", "metadata": {}, "source": [ "## 2D image" ] }, { "cell_type": "code", "execution_count": null, "id": "96237208-46a8-41d7-b633-a6330a6701ae", "metadata": {}, "outputs": [], "source": [ "d = AstroData(fitsimage='test2D.fits', Tb=True, sigma=5e-3)\n", "f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')\n", "f.read(d)\n", "p = pad(rmax=0.8, center='ICRS 04h04m43.07s 26d18m56.20s')\n", "p.add_color(**d.todict(), cblabel='Tb (K)')\n", "p.add_contour(fitsimage='test2D_2.fits', colors='r', sigma=5e-3)\n", "p.add_contour(fitsimage='test2D.fits', xskip=2, yskip=2, sigma=5e-3)\n", "p.add_segment(ampfits='test2Damp.fits',\n", " angfits='test2Dang.fits',\n", " xskip=3, yskip=3, sigma=None)\n", "p.add_scalebar(length=50 / 140, label='50 au')\n", "p.add_text([0.3, 0.3], slist='text')\n", "p.add_marker('04h04m43.07s 26d18m56.20s')\n", "p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5])\n", "p.add_arrow([0.4, 0.4], anglelist=150, rlist=0.5)\n", "p.add_region('ellipse', [0.2, 0.8], majlist=0.4, minlist=0.2, palist=45)\n", "p.set_axis_radec(nticksminor=5, title={'label': '2D image', 'loc': 'right'})\n", "p.savefig('test2D.png', show=True)" ] }, { "cell_type": "markdown", "id": "754470b5-1eef-4648-a383-dc5ee473ec1c", "metadata": {}, "source": [ "## 3D channel maps" ] }, { "cell_type": "code", "execution_count": null, "id": "79f64a2d-be44-48a2-92c2-416a0db3c9f3", "metadata": {}, "outputs": [], "source": [ "p = pad(rmax=0.8, fitsimage='test3D.fits', vmin=-5, vmax=5, vskip=2)\n", "p.add_color(fitsimage='test3D.fits',\n", " stretch='log', sigma='hist,edge')\n", "p.add_contour(fitsimage='test3D.fits',\n", " colors='r', sigma='hist,edge')\n", "p.add_contour(fitsimage='test2D.fits', colors='b', sigma=5e-3)\n", "p.add_segment(ampfits='test2Damp.fits',\n", " angfits='test2Dang.fits',\n", " xskip=3, yskip=3, sigma=None)\n", "p.add_scalebar(length=50 / 140, label='50 au')\n", "p.add_text([[0.3, 0.3]], slist=['text'], include_chan=[0, 1, 2])\n", "p.add_marker([0.7, 0.7])\n", "p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5],\n", " include_chan=[6, 7, 8])\n", "p.add_arrow([[0.4, 0.4]], anglelist=[150], rlist=[0.5],\n", " include_chan=[9, 10, 11])\n", "p.add_region('rectangle', [[0.2, 0.8]],\n", " majlist=[0.4], minlist=[0.2], palist=[45],\n", " include_chan=[12, 13, 14])\n", "p.set_axis(grid={}, title='3D channel maps')\n", "p.savefig('test3D.png', show=True)" ] }, { "cell_type": "markdown", "id": "ac70e517-6750-4d05-8f5f-2157f609f41c", "metadata": {}, "source": [ "## PV diagram" ] }, { "cell_type": "code", "execution_count": null, "id": "3651c10f-b0f0-4c1b-9bd7-d91f67d33d36", "metadata": {}, "outputs": [], "source": [ "d = AstroData(fitsimage='testPV.fits', Tb=True, sigma=None, pvpa=60)\n", "f = AstroFrame(pv=True)\n", "f.read(d)\n", "p = pad(rmax=0.8, pv=True, swapxy=True, vmin=-5, vmax=5, figsize=(6, 7))\n", "p.add_color(**d.todict(), cblabel='Tb (K)', cblocation='top')\n", "p.add_contour(fitsimage='testPV.fits', sigma=1e-3, pvpa=60,\n", " colors='r')\n", "p.add_text([0.3, 0.3], slist='text')\n", "p.add_marker([[0.5, 0.5]])\n", "p.set_axis(title='PV diagram')\n", "p.savefig('testPV.png', show=True)" ] }, { "cell_type": "markdown", "id": "547bee11-fc91-402e-8818-8fbb2515e5cd", "metadata": {}, "source": [ "## log log PV" ] }, { "cell_type": "code", "execution_count": null, "id": "7c840b8b-d336-4d33-b508-18c53984a156", "metadata": {}, "outputs": [], "source": [ "p = pad(rmax=0.8 * 140, pv=True, quadrants='13', vmin=-5, vmax=5, dist=140)\n", "p.add_color(fitsimage='testPV.fits', Tb=True, sigma=None,\n", " cblabel='Tb (K)', show_beam=False, pvpa=60)\n", "p.add_contour(fitsimage='testPV.fits', colors='r',\n", " sigma=1e-3, show_beam=False, pvpa=60)\n", "p.set_axis(title='loglog PV diagram', loglog=20)\n", "p.savefig('testloglogPV.png', show=True)" ] }, { "cell_type": "markdown", "id": "3a2f7ac5-1907-4174-923e-c94d359eb8c4", "metadata": {}, "source": [ "## RGB" ] }, { "cell_type": "code", "execution_count": null, "id": "95401b24-0635-45ca-93d8-8207c21d2c57", "metadata": {}, "outputs": [], "source": [ "d = AstroData(fitsimage='test3D.fits', Tb=True, sigma=5e-3)\n", "f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')\n", "f.read(d)\n", "dblue = np.sum(d.data[0:20], axis=0) * d.dv\n", "dgreen = np.sum(d.data[20:41], axis=0) * d.dv\n", "dred = np.sum(d.data[41:61], axis=0) * d.dv\n", "d.data = [dred, dgreen, dblue]\n", "d.sigma = [d.sigma * d.dv * np.sqrt(20)] * 3\n", "p = pad(rmax=0.8, center='ICRS 04h04m43.07s 26d18m56.20s')\n", "p.add_rgb(**d.todict())\n", "p.add_contour(fitsimage='test2D_2.fits', colors='r', sigma=5e-3)\n", "p.add_contour(fitsimage='test2D.fits', xskip=2, yskip=2, sigma=5e-3)\n", "p.add_segment(ampfits='test2Damp.fits',\n", " angfits='test2Dang.fits',\n", " xskip=3, yskip=3, sigma=None)\n", "p.add_scalebar(length=50 / 140, label='50 au')\n", "p.add_text([0.3, 0.3], slist='text')\n", "p.add_marker('04h04m43.07s 26d18m56.20s')\n", "p.add_line([[0.5, 0.5], [0.6, 0.6]], anglelist=[60, 60], rlist=[0.5, 0.5])\n", "p.add_arrow([0.4, 0.4], anglelist=150, rlist=0.5)\n", "p.add_region('ellipse', [0.2, 0.8], majlist=0.4, minlist=0.2, palist=45)\n", "p.set_axis_radec(nticksminor=5, title={'label': '2D RGB', 'loc': 'right'})\n", "p.savefig('test2Drgb.png', show=True)" ] }, { "cell_type": "markdown", "id": "cb02606b-0861-4f26-98d8-acf4fd059104", "metadata": {}, "source": [ "## Line profile" ] }, { "cell_type": "code", "execution_count": null, "id": "05746870-af9e-4c59-bf36-3c72a941b5ac", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.plot_utils import plotprofile\n", "\n", "plotprofile(fitsimage='test3D.fits', ellipse=[[0.2, 0.2, 0]] * 2,\n", " flux=True, sigma='hist,edge',\n", " coords=['04h04m43.045s 26d18m55.766s',\n", " '04h04m43.109s 26d18m56.704s'],\n", " gaussfit=True, savefig='testprofile.png', show=True, width=2)" ] }, { "cell_type": "markdown", "id": "a49a6f9e-b047-4a22-bb0f-c9463b37a2a8", "metadata": {}, "source": [ "## Spatial slice" ] }, { "cell_type": "code", "execution_count": null, "id": "a2c60d38-041e-4840-a52e-8dae24d70dfe", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.plot_utils import plotslice\n", "\n", "plotslice(length=1.6, pa=270, fitsimage='test2D.fits',\n", " center='04h04m43.07s 26d18m56.20s', sigma=5e-3,\n", " savefig='testslice.png', show=True)" ] }, { "cell_type": "markdown", "id": "72bcee23", "metadata": {}, "source": [ "## 3D html figure (Plotly)" ] }, { "cell_type": "code", "execution_count": null, "id": "0b3a2a99", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.plot_utils import plot3d\n", "\n", "d = AstroData(fitsimage='test3D.fits', Tb=True, sigma=5e-3)\n", "f = AstroFrame(rmax=0.8, center='B1950 04h01m40.5705s +26d10m47.285s')\n", "f.read(d)\n", "d.data = np.sum(d.data, axis=0) * d.dv\n", "d.sigma = d.sigma * d.dv * np.sqrt(len(d.v))\n", "vplus = {'data': d.data, 'sigma': d.sigma, 'cmap': 'gray'}\n", "plot3d(rmax=0.8, vmin=-5, vmax=5, fitsimage='test3D.fits',\n", " outname='./_static/plotly/test3D.html', levels=[3, 6, 9],\n", " sigma='hist,edge', vplus=vplus, show=False)" ] }, { "cell_type": "markdown", "id": "6968f63a", "metadata": {}, "source": [ "\n", "" ] }, { "cell_type": "markdown", "id": "f62edbab", "metadata": {}, "source": [ "## Animation" ] }, { "cell_type": "code", "execution_count": null, "id": "04bfcb8f", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib.animation as animation\n", "\n", "nchans = 16\n", "\n", "def update_plot(i):\n", " p = pad(rmax=0.8, fitsimage='test3D.fits',\n", " vmin=-5, vmax=5, vskip=4,\n", " channelnumber=i, fig=fig)\n", " p.add_color(fitsimage='test3D.fits',\n", " stretch='log', sigma='hist,edge')\n", " p.add_scalebar(length=50 / 140, label='50 au')\n", " p.set_axis_radec(grid={}, title='3D channel maps')\n", " p.fig.tight_layout()\n", "\n", "fig = plt.figure(figsize=(7, 5)) # Giving figsize may help ffmpeg.\n", "ani = animation.FuncAnimation(fig, update_plot, frames=nchans)\n", "#Writer = animation.writers['ffmpeg'] # for mp4\n", "Writer = animation.writers['pillow'] # for gif\n", "writer = Writer(fps=1) # frame per second\n", "#ani.save('./_static/video/test_animation.mp4', writer=writer, dpi=64)\n", "ani.save('./_static/video/test_animation.gif', writer=writer, dpi=64)\n", "plt.close()" ] }, { "cell_type": "markdown", "id": "9eab5a88-af48-45c3-9fe6-69c1133986f5", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "2a61d241", "metadata": {}, "source": [ "## Noise estimate" ] }, { "cell_type": "code", "execution_count": null, "id": "b243e907", "metadata": {}, "outputs": [], "source": [ "\n", "from plotastrodata.noise_utils import Noise\n", "\n", "x = np.linspace(-1.5, 1.5, 301)\n", "y = np.linspace(-1.5, 1.5, 301)\n", "x, y = np.meshgrid(x, y)\n", "r = np.hypot(x, y)\n", "data = np.random.randn(*np.shape(r))\n", "data[r > 1.5] = np.nan\n", "data = data * np.exp2(r**2)\n", "n = Noise(data=data, sigma='hist-pbcor')\n", "n.gen_histogram()\n", "n.fit_histogram()\n", "print('Noise (mean, std, Rout):',\n", " [n.mean, n.std, float(n.popt[2])])\n", "n.plot_histogram(savefig='noise.png', show=True)" ] }, { "cell_type": "markdown", "id": "08579523-3973-4396-833e-e6d49c4ce501", "metadata": {}, "source": [ "## Line-of-sight velocity with 3D rotation" ] }, { "cell_type": "code", "execution_count": null, "id": "7c5b52c3-edc3-46fd-88c3-1a95c3520f51", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.los_utils import sys2obs, polarvel2losvel\n", "\n", "incl = 60\n", "theta0 = 60\n", "\n", "xscale = np.sin(np.radians(theta0))**2\n", "vscale = 1 / np.sin(np.radians(theta0))\n", "\n", "xsys = np.linspace(0, 3 / xscale, 200)\n", "ysys = np.sqrt(2 * xsys + 1)\n", "zsys = np.zeros_like(xsys)\n", "\n", "r = np.hypot(xsys, ysys)\n", "theta = np.zeros_like(r) + np.pi / 2\n", "phi = np.arctan2(ysys, xsys)\n", "\n", "v_r = -np.sqrt(1 / r * (2 - 1 / r))\n", "v_theta = np.zeros_like(v_r)\n", "v_phi = 1 / r\n", "\n", "xsys = xsys * xscale\n", "ysys = ysys * xscale\n", "zsys = zsys * xscale\n", "\n", "v_r = v_r * vscale\n", "v_theta = v_theta * vscale\n", "v_phi = v_phi * vscale\n", "\n", "xlist, ylist, vlist = [], [], []\n", "for phi0 in np.linspace(0, 360, 9):\n", " xobs, yobs, zobs = sys2obs(xsys=xsys, ysys=ysys, zsys=zsys,\n", " incl=incl, phi0=phi0, theta0=theta0)\n", " vlos = polarvel2losvel(v_r=v_r, v_theta=v_theta, v_phi=v_phi,\n", " theta=theta, phi=phi,\n", " incl=incl, phi0=phi0, theta0=theta0)\n", " xlist.append(xobs)\n", " ylist.append(yobs)\n", " vlist.append(vlos)\n", "\n", "fig, ax = plt.subplots()\n", "m = ax.scatter(xlist, ylist, c=vlist, cmap='coolwarm', vmin=-1.5, vmax=1.5)\n", "fig.colorbar(m, ax=ax, label=r'$V_\\mathrm{los} / \\sqrt{GM_{*}/r_{c}}$')\n", "ax.set_xlim(3, -3)\n", "ax.set_ylim(-3, 3)\n", "ax.set_aspect(1)\n", "ax.set_xlabel(r'$x / r_{c}$')\n", "ax.set_ylabel(r'$y / r_{c}$')\n", "ax.grid()\n", "fig.savefig('streamer.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "be1afdf5-1c64-4ed2-be69-7791cbdee1b4", "metadata": {}, "source": [ "## FFT with a given center" ] }, { "cell_type": "code", "execution_count": null, "id": "8fd0eee0-1d1a-4727-bf61-189f7c8fbcee", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.fft_utils import fftcentering\n", "\n", "x = np.linspace(-99.5, 99.5, 200)\n", "f = np.where(np.abs(x)<10, 1, 0)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(x, f)\n", "ax.set_xlabel('x')\n", "ax.set_ylabel('f')\n", "fig.tight_layout()\n", "fig.savefig('boxcar.png')\n", "plt.show()\n", "\n", "u = np.fft.fftshift(np.fft.fftfreq(len(x), d=x[1] - x[0]))\n", "F = np.fft.fftshift(np.fft.fft(f))\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(u, np.real(F), label='real')\n", "ax.plot(u, np.imag(F), label='imag')\n", "ax.set_xlabel('u')\n", "ax.set_ylabel('numpy.fft')\n", "ax.legend()\n", "fig.tight_layout()\n", "fig.savefig('numpyfft.png')\n", "plt.show()\n", "\n", "F, u = fftcentering(f=f, x=x, xcenter=0)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(u, np.real(F), label='real')\n", "ax.plot(u, np.imag(F), label='imag')\n", "ax.set_xlabel('u')\n", "ax.set_ylabel('fftcentering')\n", "ax.legend()\n", "fig.tight_layout()\n", "fig.savefig('fftcentering.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "870ef0b5-7407-40f1-b8f0-5942077a18eb", "metadata": {}, "source": [ "## MCMC (emcee, corner, dynesty)" ] }, { "cell_type": "code", "execution_count": null, "id": "20c524af-24ed-4abd-b78f-a9ade8cedf34", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.fitting_utils import EmceeCorner\n", "from plotastrodata.plot_utils import set_rcparams\n", "\n", "set_rcparams(fontsize=12)\n", "\n", "def logl(p):\n", " x1, x2, x3 = p\n", " chi2 = (x1 / 1)**2 + (x2 / 2)**2 + (x3 / 4)**2\n", " return -0.5 * chi2\n", "\n", "fitter = EmceeCorner(bounds=[[-5, 5], [-10, 10], [-20, 20]],\n", " logl=logl, progressbar=False, percent=[16, 84])\n", "\n", "fitter.fit(nwalkersperdim=30, nsteps=11000, nburnin=1000,\n", " #savechain='chain.npy'\n", " )\n", "print('best:', fitter.popt)\n", "print('lower percentile:', fitter.plow)\n", "print('50 percentile:', fitter.pmid)\n", "print('higher percentile:', fitter.phigh)\n", "fitter.getDNSevidence()\n", "print('evidence:', fitter.evidence)\n", "fitter.plotcorner(show=True, savefig='corner.png',\n", " labels=['par1', 'par2', 'par3'],\n", " cornerrange=[[-4, 4], [-8, 8], [-16, 16]])\n", "fitter.plotchain(show=True, savefig='chain.png',\n", " labels=['par1', 'par2', 'par3'],\n", " ylim=[[-2, 2], [-4, 4], [-8, 8]])" ] }, { "cell_type": "markdown", "id": "ad7e33ab-b1ec-44c8-a1c9-708ef67be7b0", "metadata": {}, "source": [ "## Likelihood on a parameter grid" ] }, { "cell_type": "code", "execution_count": null, "id": "974fefc9-1673-48ef-9513-fa203a289846", "metadata": {}, "outputs": [], "source": [ "fitter.posteriorongrid(ngrid=[101, 201, 401])\n", "print('best:', fitter.popt)\n", "print('lower percentile:', fitter.plow)\n", "print('50th percentile:', fitter.pmid)\n", "print('higher percentile:', fitter.phigh)\n", "print('evidence:', fitter.evidence)\n", "fitter.plotongrid(show=True, savefig='grid.png',\n", " labels=['par1', 'par2', 'par3'],\n", " cornerrange=[[-4, 4], [-8, 8], [-16, 16]])" ] }, { "cell_type": "markdown", "id": "11713e93-c805-4117-b7ef-b3c936ea7bac", "metadata": {}, "source": [ "## Other functions" ] }, { "cell_type": "code", "execution_count": null, "id": "f6ba9125-2535-42ce-98d7-2f81f01237e8", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.coord_utils import xy2coord, coord2xy\n", "\n", "coord = xy2coord(xy=[[30, 90], [0, 0]], coordorg='00h00m00s 60d00m00s')\n", "print(coord)\n", "xy = coord2xy(coords=coord, coordorg='00h00m00s 60d00m00s')\n", "print(np.round(xy, 2))" ] }, { "cell_type": "code", "execution_count": null, "id": "ea947c40-3486-4e9e-b567-8ddb5cd9c22e", "metadata": {}, "outputs": [], "source": [ "import plotastrodata.const_utils as cu\n", "\n", "print(cu.pc)\n", "print(cu.M_sun)\n", "print(cu.centi)" ] }, { "cell_type": "code", "execution_count": null, "id": "c6f2a4d7-2dfe-46b3-9670-f07ce1688bc8", "metadata": {}, "outputs": [], "source": [ "from plotastrodata.ext_utils import BnuT, JnuT\n", "\n", "print(BnuT(T=30, nu=230e9))\n", "print(JnuT(T=30, nu=230e9))" ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.3" } }, "nbformat": 4, "nbformat_minor": 5 }