diff --git a/tutorials/visualization-with-gaia/requirements.txt b/tutorials/visualization-with-gaia/requirements.txt new file mode 100644 index 00000000..29476e6a --- /dev/null +++ b/tutorials/visualization-with-gaia/requirements.txt @@ -0,0 +1,4 @@ +numpy +matplotlib +astropy +astroquery diff --git a/tutorials/visualization-with-gaia/visualization-with-gaia.ipynb b/tutorials/visualization-with-gaia/visualization-with-gaia.ipynb new file mode 100644 index 00000000..84f977f7 --- /dev/null +++ b/tutorials/visualization-with-gaia/visualization-with-gaia.ipynb @@ -0,0 +1,463 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6f85a9d1-f7b6-4b6a-b9d3-c0cfa04e8b17", + "metadata": {}, + "source": [ + "# Vizualising high proper motion sources using Gaia\n", + "\n", + "## Author\n", + "C. E. Brasseur\n", + "\n", + "## Learning Goals\n", + "\n", + "1. Use Astroquery to access GAIA data\n", + "2. Properly format an ADQL query\n", + "3. Use the Astropy SkyCoord class to represent stars with known proper motion and radial velocity values\n", + "4. Visualize in 2- and 3D star positions on the sky\n", + "5. Using the matplotlib animation library, visualize the changing position of stars through time \n", + "\n", + "## Keywords\n", + "\n", + "astropy.coordinates, astroquery.gaia, matplotlib.animation\n", + "\n", + "\n", + "## Companion Content\n", + "\n", + "- [Gaia Archive](https://gea.esac.esa.int/archive/)\n", + "- [Gaia Astroquery Module](https://astroquery.readthedocs.io/en/latest/gaia/gaia.html)\n", + "- [Gaia source catalog column definitions](https://gea.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_source_catalogue/ssec_dm_gaia_source.html)\n", + "
\n", + "\n", + "- [SkyCoord Documentation](https://docs.astropy.org/en/stable/coordinates/skycoord.html)\n", + "- [Matplotlib 3D Plots](https://matplotlib.org/stable/users/explain/toolkits/mplot3d.html)\n", + "- [Matplotlib Animation](https://matplotlib.org/stable/users/explain/animations/animations.html)\n", + "\n", + "## Summary\n", + "\n", + "In this tutorial we will use the astroquery package to get a list of stars in the Gaia catalogue with high proper motions and measured radial velocities. We will then visualize the resulting star positions on a 2D Aitoff-projected grid, and in 3 dimensions. Next we will use the functionality built into the Astropy coordinates framework to project the path of our collection of stars forward in time. Finally, we will visualize the movement of the stars by creating an animation within matplotlib.\n", + "\n", + "1. [Imports](#Imports)\n", + "2. [Query Gaia](#Query-Gaia)\n", + "3. [Visualize star positions](#Visualize-star-positions)\n", + "4. [Animate stellar trajectories](#Animate-stellar-trajectories)" + ] + }, + { + "cell_type": "markdown", + "id": "8d1cc023", + "metadata": {}, + "source": [ + "## Imports\n", + "\n", + "We will use ``astroquery.gaia`` for data access, ``astropy.coordinates`` for sky position representation and manipulation, and ``matplotlib`` for visualization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bc7cc9a-2a99-467a-8eff-8e8074d4beb6", + "metadata": {}, + "outputs": [], + "source": [ + "import warnings # So we can suppress expected warnings\n", + "import numpy as np\n", + "\n", + "from astroquery.gaia import Gaia # For data access\n", + "\n", + "from astropy.coordinates import SkyCoord # For storing a manipulation object sky positions\n", + "import astropy.units as u\n", + "\n", + "# For plotting\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.animation as animation\n", + "\n", + "plt.rcParams[\"animation.html\"] = \"jshtml\" # To make the animations render correctly in the notebook\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "e193bd69", + "metadata": {}, + "source": [ + "## Query Gaia\n", + "\n", + "The [``astroquery.gaia``](https://astroquery.readthedocs.io/en/latest/gaia/gaia.html) module uses the Astronomical Data Query Language ([ADQL](https://www.ivoa.net/documents/ADQL/2.0)) to query its various databases. Here we query the ``gaia_source`` table in the Gaia Data Release 3 database, selecting a subset of [columns](https://gea.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_source_catalogue/ssec_dm_gaia_source.html) to be returned and adding some conditions for the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b35cfad-088c-4e7a-b888-ede3d3f4b558", + "metadata": {}, + "outputs": [], + "source": [ + "adql_query = (\"SELECT source_id,ra,dec,pmra,pmdec,radial_velocity,distance_gspphot \" # The columns we want\n", + " \"FROM gaiadr3.gaia_source \" # The table we are querying\n", + " \"WHERE pm>=1000 \" # Only return rows where the proper motion is >= 1000 mas/yr\n", + " \"AND distance_gspphot<=40 \" # Only return rows where photomentric distance is <= 40 pc\n", + " \"AND radial_velocity IS NOT NULL)\" # Only return rows where there is a radial velocity measurement\n", + " )\n", + "\n", + "job = Gaia.launch_job(adql_query)\n", + "\n", + "gaia_table = job.get_results()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44eacde8", + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Number of objects found: {len(gaia_table)}\\n\")\n", + "\n", + "print(\"First 5 rows:\")\n", + "gaia_table[:5]" + ] + }, + { + "cell_type": "markdown", + "id": "5875404e", + "metadata": {}, + "source": [ + "For the rest of this notebook we are going to take advantage of the functionality built into Astropy's [``SkyCoord``](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html) class. In order to do this we are going to combine not only the star positions (including distance), but their proper motions and radial velocities to create a ``SkyCoord`` object that not only knows where the stars are at present but how they are moving." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97f69e16", + "metadata": {}, + "outputs": [], + "source": [ + "gaia_table[\"position\"] = SkyCoord(gaia_table['ra'], gaia_table['dec'], # star sky positions\n", + " distance=gaia_table[\"distance_gspphot\"], # star distances\n", + "\n", + " # proper motion, conventional format multiplies the ra proper motion \n", + " # by the cosine of the declination\n", + " pm_ra_cosdec=u.Quantity(gaia_table['pmra'])*np.cos(gaia_table['dec']), \n", + " pm_dec=gaia_table['pmdec'],\n", + " \n", + " # Radial velocity\n", + " radial_velocity=gaia_table['radial_velocity'])" + ] + }, + { + "cell_type": "markdown", + "id": "2c746cd8", + "metadata": {}, + "source": [ + "## Visualize star positions\n", + "\n", + "In this section we use the stellar positions at the present time to visualize the spread of stars on this sky.\n", + "To do this, we write a function that takes in a list of object coordinatates and builds a figure that visualized the stars in two ways: projected onto the plane of the sky, and in three dimensional space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16240bae", + "metadata": {}, + "outputs": [], + "source": [ + "def star_plot(star_coords, vmin=None, vmax=None):\n", + " \"\"\"\n", + " Take a list of stars as SkyCoord objects and build a figure showing the stars projected onto\n", + " the plane of the sky using the Aitoff projection, and in three dimensional space centered \n", + " on the Earth. In both plots the points are colored by distance (from Earth).\n", + " \n", + " Parameters\n", + " ----------\n", + " star_coords : Astropy SkyCoord object\n", + " The coordinates of the objects to be plotted\n", + " vmin, vmax : float\n", + " Optional min and max values in pc for the colormap used to color the plot by distance. \n", + " If not set they will be the minimum and maximum distances of star_coords.\n", + " \n", + " Returns\n", + " -------\n", + " response : matplotlib figure\n", + " \n", + " \"\"\"\n", + " \n", + " # Set the min and max values for the colormap if not given\n", + " if not vmin:\n", + " vmin = star_coords.distance.min().value\n", + " \n", + " if not vmax:\n", + " vmax = star_coords.distance.max().value\n", + " \n", + " # Initialize the figure\n", + " fig = plt.figure(figsize=(14, 7))\n", + " \n", + " # Sky Projection plot\n", + " ax1 = fig.add_subplot(1,2,1, projection='aitoff') # Add left subplot\n", + " \n", + " # Turn off the axis ticks and labels, and turn on the plot grid\n", + " ax1.set_xticklabels([])\n", + " ax1.set_yticklabels([])\n", + " ax1.grid(True)\n", + " \n", + " # Plot the stars, colored by distance\n", + " ax1.scatter(star_coords.ra.wrap_at(180 * u.deg).radian, star_coords.dec.radian,\n", + " c=star_coords.distance, vmin=vmin, vmax=vmax, marker=\"*\", s=200, cmap=\"plasma\")\n", + " \n", + " # 3D plot\n", + " ax2 = fig.add_subplot(1,2,2, projection='3d') # Add right subplot\n", + " \n", + " # Turn off axis ticks and labels\n", + " ax2.set_xticklabels([])\n", + " ax2.set_yticklabels([])\n", + " ax2.set_zticklabels([])\n", + " \n", + " # Plotting the Earth as a black circle\n", + " ax2.scatter([0],[0],[0],s=100,c='k',marker=\"o\")\n", + "\n", + " # Plot the stars, colored by distance\n", + " pc = ax2.scatter(star_coords.cartesian.x, star_coords.cartesian.y, star_coords.cartesian.z, \n", + " c=star_coords.distance, s=200, vmin=vmin, vmax=vmax, marker=\"*\", cmap=\"plasma\")\n", + " \n", + " # Adding the colorbar\n", + " cbar = fig.colorbar(pc, shrink=0.6, location=\"right\")\n", + " cbar.ax.tick_params(labelsize=14)\n", + " cbar.ax.set_ylabel('Distance (pc)', fontsize=18)\n", + " \n", + " # Remove extra space between the subplots\n", + " fig.subplots_adjust(wspace=0)\n", + " \n", + " # This prevents getting an extra copy of the plot when you call the function\n", + " plt.close()\n", + " \n", + " return fig" + ] + }, + { + "cell_type": "markdown", + "id": "7f039f92", + "metadata": {}, + "source": [ + "First we plot all the stars in our query result.\n", + "\n", + "In the below plot we can see that we have good sky coverage, and the 3D plot clearly shows that we are indeed coloring the points by distance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "955d7d6f", + "metadata": {}, + "outputs": [], + "source": [ + "star_plot(gaia_table[\"position\"])" + ] + }, + { + "cell_type": "markdown", + "id": "80c1046e", + "metadata": {}, + "source": [ + "We can also plot a subset of the stars, for example, below we plot only a single quadrant of the sky by imposing the condition that stars must have RA$~< 180^\\circ$ and Dec $> 0^\\circ$.\n", + "\n", + "Here we can clearly see in the plot on the left the single quandrant we are plotting, and on the right how those stars spread out from the Earth." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2cca46ee", + "metadata": {}, + "outputs": [], + "source": [ + "star_plot(gaia_table[\"position\"][(gaia_table[\"ra\"] < 180) & (gaia_table[\"dec\"] > 0)])" + ] + }, + { + "cell_type": "markdown", + "id": "37e3e5bb", + "metadata": {}, + "source": [ + "## Animate stellar trajectories\n", + "\n", + "In the previous section we plotted the stars in their current positions. However, we also included in our ``SkyCoord`` object information about how the stars are moving, and in this section we will use that infromation to evolve their positions through time. \n", + "\n", + "The function we use is [``apply_space_motion``](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord.apply_space_motion), which allows us to input a length of time (``dt``) and get back the coordinates for the object(s) after that span of time. \n", + "\n", + "We will write a function to take the plot we created in the previous section and animate it, showing the trajectory of our stars as time passes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "198d9709", + "metadata": {}, + "outputs": [], + "source": [ + "def star_animation(star_coords, evolution_time, steps, vmin=None, vmax=None):\n", + " \"\"\"\n", + " Take a list of stars as a SkyCoord object and build an animation showing how the star's\n", + " positions evolve over time. \n", + " \n", + " Two subplots are produced, the left showing the stars projected onto the plane of the sky \n", + " using the Aitoff projection, and the right in three dimensional space centered on the Earth. \n", + " In both plots the points are colored by distance (from Earth).\n", + " \n", + " Parameters\n", + " ----------\n", + " star_coords : Astropy SkyCoord object\n", + " The coordinates of the objects you wish to plot\n", + " evolution_time : float or Astropy Quantity\n", + " The amount of time to evolve the stellar positions each step. \n", + " If not specified as a quantity, years will be assumed.\n", + " steps : int\n", + " The number of time steps to take. \n", + " The total amount of time the system will be eveolved for is evolution_time * steps\n", + " vmin, vmax : float\n", + " Optional min and max values for the colormap used to color the plot by distance. \n", + " If not set they will be the minimum and maximum distances of the input stars at present time.\n", + " \n", + " Returns\n", + " -------\n", + " response : matplotlib figure\n", + " \n", + " \"\"\"\n", + " \n", + " # Make sure the evolution time is in years\n", + " if not isinstance(evolution_time, u.Quantity):\n", + " evolution_time *= u.yr\n", + " \n", + " # Create array of time deltas (all in relation to present time)\n", + " dt_array = np.linspace(0, evolution_time*steps, steps, endpoint=False)\n", + " \n", + " # If vmin/vmax were not set we need to set them (if they are not set at all, \n", + " # each step of the animation will have different colormap boundaries).\n", + " if not vmin:\n", + " vmin = star_coords.distance.min().value\n", + " \n", + " if not vmax:\n", + " vmax = star_coords.distance.max().value\n", + " \n", + "\n", + " def animate(iteration, dt_array, aitoff_scatter, cube_scatter):\n", + " \"\"\"\n", + " This function handles updating the plot for each frame of the animation.\n", + " \n", + " Parameters\n", + " ----------\n", + " iteration: int\n", + " Current iteration (frame number)\n", + " dt_array : Quantity\n", + " Array of time deltas\n", + " aitoff_scatter : matplotlib PathCollection\n", + " A matplotlib collection object that holds the 2D (left) plot data points.\n", + " cube_scatter : matplotlib Path3DCollection\n", + " A matplotlib collection object that holds the 3D (right) plot data points.\n", + " \"\"\"\n", + "\n", + " with warnings.catch_warnings():\n", + " # Projecting more than 5 years into the future gives a \"dubious year\" warning \n", + " # due to the unpredictability of leap seconds, we suppress this warning becuase \n", + " # we are not concerned with temporal accuracy to the second\n", + " warnings.filterwarnings(\"ignore\", message=\"ERFA function \")\n", + " new_pos = star_coords.apply_space_motion(dt=dt_array[iteration])\n", + " \n", + " # set_offsets sets the point locations to the newly calculated coordinates\n", + " aitoff_scatter.set_offsets(list(zip(new_pos.ra.wrap_at(180 * u.deg).radian, new_pos.dec.radian)))\n", + " \n", + " # set_array sets the colors of the points based on their new distances\n", + " aitoff_scatter.set_array(new_pos.distance)\n", + "\n", + " # For the 3D plot we have to use the private property _offsets3d to update the point locations\n", + " cube_scatter._offsets3d = (new_pos.cartesian.x, new_pos.cartesian.y, new_pos.cartesian.z)\n", + " \n", + " # set_array sets the colors of the points based on their new distances\n", + " cube_scatter.set_array(new_pos.distance)\n", + "\n", + " return aitoff_scatter, cube_scatter,\n", + " \n", + " # Call our star_plot function from the previous section to initialize the figure\n", + " fig = star_plot(star_coords, vmin=vmin, vmax=vmax)\n", + " \n", + " # Get the matplotlib Collection objects that correspond to our sets of points\n", + " aitoff_scatter = fig.axes[0].collections[0]\n", + " cube_scatter = fig.axes[1].collections[1]\n", + " \n", + " # Build the animation\n", + " ani = animation.FuncAnimation(fig, # The figure we are animating\n", + " animate, # The function that updates the plot for each frame\n", + " fargs=(dt_array, aitoff_scatter, cube_scatter), # Args for the animate function\n", + " frames=steps, # The number of frames in our animation\n", + " interval=100) # Delay between frames in milliseconds\n", + " \n", + " # This prevents getting an extra copy of the first frame of the plot when you call the function\n", + " plt.close()\n", + " \n", + " return ani" + ] + }, + { + "cell_type": "markdown", + "id": "3e7ebece", + "metadata": {}, + "source": [ + "Here we show the evolution of the stellar positions over 20,000 years by animating 20 steps of 1000 years each." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abbf3d25", + "metadata": {}, + "outputs": [], + "source": [ + "star_animation(gaia_table[\"position\"], 1000*u.yr, 20)" + ] + }, + { + "cell_type": "markdown", + "id": "7ce1db86", + "metadata": {}, + "source": [ + "As we did for the still plots, we also plot a single quadrant of the sky by imposing the condition that stars must have RA$~< 180^\\circ$ and Dec $> 0^\\circ$. We also increase the time step from 1,000 years to 10,000 years, meaning the full evolution time we animate is 200,000 years.\n", + "\n", + "Because some of the stars are moving away from us and we are pushing further forward in time we set the colormap range to $4-60 \\mathrm{~pc}$ to cover the larger range of distances we will need to represent as the stars move." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe2b8dbb", + "metadata": {}, + "outputs": [], + "source": [ + "star_animation(gaia_table[\"position\"][(gaia_table[\"ra\"] < 180) & (gaia_table[\"dec\"] > 0)], 10_000*u.yr, \n", + " 20, vmin=5, vmax=60)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05b1d75c", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython" + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}