In the previous lesson, we wrote a query to select stars from the region of the sky where we expect GD-1 to be, and saved the results in a FITS file.
Now we'll read that data back and implement the next step in the analysis, identifying stars with the proper motion we expect for GD-1.
Here are the steps in this lesson:
We'll read back the results from the previous lesson, which we saved in a FITS file.
Then we'll transform the coordinates and proper motion data from ICRS back to the coordinate frame of GD-1.
We'll put those results into a Pandas
DataFrame, which we'll use to select stars near the centerline of GD-1.
Plotting the proper motion of those stars, we'll identify a region of proper motion for stars that are likely to be in GD-1.
Finally, we'll select and plot the stars whose proper motion is in that region.
After completing this lesson, you should be able to
Select rows and columns from an Astropy
Use Matplotlib to make a scatter plot.
Use Gala to transform coordinates.
Make a Pandas
DataFrame and use a Boolean
Series to select rows.
DataFrame in an HDF5 file.
If you are running this notebook on Colab, you can run the following cell to install the libraries we'll use.
If you are running this notebook on your own computer, you might have to install these libraries yourself. See the instructions in the preface.
# If we're running on Colab, install libraries import sys IN_COLAB = 'google.colab' in sys.modules if IN_COLAB: !pip install astroquery astro-gala
In the previous lesson, we ran a query on the Gaia server and downloaded data for roughly 140,000 stars. We saved the data in a FITS file so that now, picking up where we left off, we can read the data from a local file rather than running the query again.
If you ran the previous lesson successfully, you should already have a file called
gd1_results.fits that contains the data we downloaded.
If not, you can download the file or run the following cell.
from os.path import basename, exists def download(url): filename = basename(url) if not exists(filename): from urllib.request import urlretrieve local, _ = urlretrieve(url, filename) print('Downloaded ' + local) download('https://github.com/AllenDowney/AstronomicalData/raw/main/' + 'data/gd1_results.fits')
Now here's how we can read the data from the file back into an Astropy
from astropy.table import Table filename = 'gd1_results.fits' results = Table.read(filename)
The result is an Astropy
We can use
info to refresh our memory of the contents.
And select an individual column like this:
The result is a
Column object that contains the data, and also the data type, units, and name of the column.
The rows in the
Table are numbered from 0 to
n is the number of rows. We can select the first row like this:
As you might have guessed, the result is a
Notice that the bracket operator selects both columns and rows. You might wonder how it knows which to select. If the expression in brackets is a string, it selects a column; if the expression is an integer, it selects a row.
If you apply the bracket operator twice, you can select a column and then an element from the column.
Or you can select a row and then an element from the row.
You get the same result either way.
To see what the results look like, we'll use a scatter plot. The library we'll use is Matplotlib, which is the most widely-used plotting library for Python. The Matplotlib interface is based on MATLAB (hence the name), so if you know MATLAB, some of it will be familiar.
We'll import like this.
import matplotlib.pyplot as plt
Pyplot is part of the Matplotlib library. It is conventional to import it using the shortened name
In recent versions of Jupyter, plots appear "inline"; that is, they are part of the notebook. In some older versions, plots appear in a new window. In that case, you might want to run the following Jupyter magic command in a notebook cell:
scatter is more versatile; for example, you can make every point in a scatter plot a different color.
plot is more limited, but for simple cases, it can be substantially faster.
Jake Vanderplas explains these differences in The Python Data Science Handbook.
Since we are plotting more than 100,000 points and they are all the same size and color, we'll use
Here's a scatter plot with right ascension on the x-axis and declination on the y-axis, both ICRS coordinates in degrees.
x = results['ra'] y = results['dec'] plt.plot(x, y, 'ko') plt.xlabel('ra (degree ICRS)') plt.ylabel('dec (degree ICRS)');
The arguments to
y, and a string that specifies the style. In this case, the letters
ko indicate that we want a black, round marker (
k is for black because
b is for blue).
ylabel put labels on the axes.
Looking at this plot, we can see that the region we selected, which is a rectangle in GD-1 coordinates, is a non-rectanglar region in ICRS coordinates.
However, this scatter plot has a problem. It is "overplotted", which means that there are so many overlapping points, we can't distinguish between high and low density areas.
To fix this, we can provide optional arguments to control the size and transparency of the points.
In the call to
plt.plot, use the keyword argument
markersize to make the markers smaller.
Then add the keyword argument
alpha to make the markers partly transparent.
Adjust these arguments until you think the figure shows the data most clearly.
Note: Once you have made these changes, you might notice that the figure shows stripes with lower density of stars. These stripes are caused by the way Gaia scans the sky, which you can read about here. The dataset we are using, Gaia Data Release 2, covers 22 months of observations; during this time, some parts of the sky were scanned more than others.
# Solution goes here
Remember that we selected data from a rectangle of coordinates in the
GD1Koposov10 frame, then transformed them to ICRS when we constructed the query.
The coordinates in
results are in ICRS.
To plot them, we will transform them back to the
GD1Koposov10 frame; that way, the axes of the figure are aligned with the orbit of GD-1, which is useful for two reasons:
By transforming the coordinates, we can identify stars that are likely to be in GD-1 by selecting stars near the centerline of the stream, where $\phi_2$ is close to 0.
By transforming the proper motions, we can identify stars with non-zero proper motion along the $\phi_1$ axis.
To do the transformation, we'll put the results into a
SkyCoord object. In a previous lesson we created a
SkyCoord object like this:
from astropy.coordinates import SkyCoord skycoord = SkyCoord(ra=results['ra'], dec=results['dec'])
Now we're going to do something similar, but in addition to
dec, we'll also include:
pmdec, which are proper motion in the ICRS frame, and
radial_velocity, which we explain below.
import astropy.units as u distance = 8 * u.kpc radial_velocity= 0 * u.km/u.s skycoord = SkyCoord(ra=results['ra'], dec=results['dec'], pm_ra_cosdec=results['pmra'], pm_dec=results['pmdec'], distance=distance, radial_velocity=radial_velocity)
For the first four arguments, we use columns from
radial_velocity we use constants, which we explain below.
The result is an Astropy
SkyCoord object, which we can transform to the GD-1 frame.
from gala.coordinates import GD1Koposov10 gd1_frame = GD1Koposov10() transformed = skycoord.transform_to(gd1_frame)
The result is another
SkyCoord object, now in the
The next step is to correct the proper motion measurements for the effect of the motion of our solar system around the Galactic center.
When we created
skycoord, we provided constant values for
radial_velocity rather than measurements from Gaia.
That might seem like a strange thing to do, but here's the motivation:
Because the stars in GD-1 are so far away, the distance estimates we get from Gaia, which are based on parallax, are not very precise. So we replace them with our current best estimate of the mean distance to GD-1, about 8 kpc. See Koposov, Rix, and Hogg, 2010.
For the other stars in the table, this distance estimate will be inaccurate, so reflex correction will not be correct. But that should have only a small effect on our ability to identify stars with the proper motion we expect for GD-1.
The measurement of radial velocity has no effect on the correction for proper motion, but we have to provide a value to avoid errors in the reflex correction calculation. So we provide
0 as an arbitrary place-keeper.
With this preparation, we can use
reflex_correct from Gala (documentation here) to correct for the motion of the solar system.
from gala.coordinates import reflex_correct skycoord_gd1 = reflex_correct(transformed)
The result is a
SkyCoord object that contains
phi2, which represent the transformed coordinates in the
pm_phi2, which represent the transformed and corrected proper motions.
We can select the coordinates and plot them like this:
x = skycoord_gd1.phi1 y = skycoord_gd1.phi2 plt.plot(x, y, 'ko', markersize=0.1, alpha=0.1) plt.xlabel('phi1 (degree GD1)') plt.ylabel('phi2 (degree GD1)');
Remember that we started with a rectangle in the GD-1 frame. When transformed to the ICRS frame, it's a non-rectangular region. Now, transformed back to the GD-1 frame, it's a rectangle again.
At this point we have two objects containing different subsets of the data.
results is the Astropy
Table we downloaded from Gaia.
skycoord_gd1 is a
SkyCoord object that contains the transformed coordinates and proper motions.
On one hand, this division of labor makes sense because each object provides different capabilities. But working with multiple object types can be awkward.
It will be more convenient to choose one object and get all of the data into it. We'll choose a Pandas
DataFrame, for two reasons:
It provides capabilities that (almost) a superset of the other data structures, so it's the all-in-one solution.
Pandas is a general-purpose tool that is useful in many domains, especially data science. If you are going to develop expertise in one tool, Pandas is a good choice.
However, compared to an Astropy
Table, Pandas has one big drawback: it does not keep the metadata associated with the table, including the units for the columns.
It's easy to convert a
Table to a Pandas
import pandas as pd results_df = results.to_pandas()
shape, which shows the number of rows and columns.
It also provides
head, which displays the first few rows.
head is useful for spot-checking large results as you go along.
shape is an attribute, so we display its value without calling it as a function;
head is a function, so we need the parentheses.
Now we can extract the columns we want from
skycoord_gd1 and add them as columns in the
phi2 contain the transformed coordinates.
results_df['phi1'] = skycoord_gd1.phi1 results_df['phi2'] = skycoord_gd1.phi2 results_df.shape
pm_phi2 contain the components of proper motion in the transformed frame.
results_df['pm_phi1'] = skycoord_gd1.pm_phi1_cosphi2 results_df['pm_phi2'] = skycoord_gd1.pm_phi2 results_df.shape
Detail: If you notice that
SkyCoord has an attribute called
proper_motion, you might wonder why we are not using it.
We could have:
proper_motion contains the same data as
pm_phi2, but in a different format.
One benefit of using Pandas is that it provides functions for exploring the data and checking for problems.
One of the most useful of these functions is
describe, which computes summary statistics for each column.
Review the summary statistics in this table.
Do the values make sense based on what you know about the context?
Do you see any values that seem problematic, or evidence of other data issues?
# Solution goes here
Now we are ready to replicate one of the panels in Figure 1 of the Price-Whelan and Bonaca paper, the one that shows components of proper motion as a scatter plot:
In this figure, the shaded area identifies stars that are likely to be in GD-1 because:
Due to the nature of tidal streams, we expect the proper motion for stars in GD-1 to be along the axis of the stream; that is, we expect motion in the direction of
phi2 to be near 0.
In the direction of
phi1, we don't have a prior expectation for proper motion, except that it should form a cluster at a non-zero value.
By plotting proper motion in the GD-1 frame, we hope to find this cluster. Then we will use the bounds of the cluster to select stars that are more likely to be in GD-1.
The following figure is a scatter plot of proper motion, in the GD-1 frame, for the stars in
x = results_df['pm_phi1'] y = results_df['pm_phi2'] plt.plot(x, y, 'ko', markersize=0.1, alpha=0.1) plt.xlabel('Proper motion phi1 (mas/yr GD1 frame)') plt.ylabel('Proper motion phi2 (mas/yr GD1 frame)');
Most of the proper motions are near the origin, but there are a few extreme values.
Following the example in the paper, we'll use
ylim to zoom in on the region near the origin.
x = results_df['pm_phi1'] y = results_df['pm_phi2'] plt.plot(x, y, 'ko', markersize=0.1, alpha=0.1) plt.xlabel('Proper motion phi1 (mas/yr GD1 frame)') plt.ylabel('Proper motion phi2 (mas/yr GD1 frame)') plt.xlim(-12, 8) plt.ylim(-10, 10);
There is a hint of an overdense region near (-7.5, 0), but if you didn't know where to look, you would miss it.
To see the cluster more clearly, we need a sample that contains a higher proportion of stars in GD-1. We'll do that by selecting stars close to the centerline.
As we can see in the following figure, many stars in GD-1 are less than 1 degree from the line
Stars near this line have the highest probability of being in GD-1.
To select them, we will use a "Boolean mask". We'll start by selecting the
phi2 column from the
phi2 = results_df['phi2'] type(phi2)
The result is a
Series, which is the structure Pandas uses to represent columns.
We can use a comparison operator,
>, to compare the values in a
Series to a constant.
phi2_min = -1.0 * u.degree phi2_max = 1.0 * u.degree mask = (phi2 > phi2_min) type(mask)
The result is a
Series of Boolean values, that is,
To select values that fall between
phi2_max, we'll use the
& operator, which computes "logical AND".
The result is true where elements from both Boolean
Series are true.
mask = (phi2 > phi2_min) & (phi2 < phi2_max)
Python detail: Python's logical operators (
not) don't work with NumPy or Pandas. Both libraries use the bitwise operators (
~) to do elementwise logical operations (explanation here).
Also, we need the parentheses around the conditions; otherwise the order of operations is incorrect.
The sum of a Boolean
Series is the number of
True values, so we can use
sum to see how many stars are in the selected region.
Series is sometimes called a "mask" because we can use it to mask out some of the rows in a
DataFrame and select the rest, like this:
centerline_df = results_df[mask] type(centerline_df)
centerline_df is a
DataFrame that contains only the rows from
results_df that correspond to
True values in
So it contains the stars near the centerline of GD-1.
We can use
len to see how many rows are in
And what fraction of the rows we've selected.
len(centerline_df) / len(results_df)
There are about 25,000 stars in this region, about 18% of the total.
Since we've plotted proper motion several times, let's put that code in a function.
def plot_proper_motion(df): """Plot proper motion. df: DataFrame with `pm_phi1` and `pm_phi2` """ x = df['pm_phi1'] y = df['pm_phi2'] plt.plot(x, y, 'ko', markersize=0.3, alpha=0.3) plt.xlabel('Proper motion phi1 (mas/yr GD1 frame)') plt.ylabel('Proper motion phi2 (mas/yr GD1 frame)') plt.xlim(-12, 8) plt.ylim(-10, 10)
And we can call it like this:
Now we can see more clearly that there is a cluster near (-7.5, 0).
You might notice that our figure is less dense than the one in the paper. That's because we started with a set of stars from a relatively small region. The figure in the paper is based on a region about 10 times bigger.
In the next lesson we'll go back and select stars from a larger region. But first we'll use the proper motion data to identify stars likely to be in GD-1.
The next step is to select stars in the "overdense" region of proper motion, which are candidates to be in GD-1.
In the original paper, Price-Whelan and Bonaca used a polygon to cover this region, as shown in this figure.
We'll use a simple rectangle for now, but in a later lesson we'll see how to select a polygonal region as well.
Here are bounds on proper motion we chose by eye:
pm1_min = -8.9 pm1_max = -6.9 pm2_min = -2.2 pm2_max = 1.0
To draw these bounds, we'll use
make_rectangle to make two lists containing the coordinates of the corners of the rectangle.
def make_rectangle(x1, x2, y1, y2): """Return the corners of a rectangle.""" xs = [x1, x1, x2, x2, x1] ys = [y1, y2, y2, y1, y1] return xs, ys
pm1_rect, pm2_rect = make_rectangle( pm1_min, pm1_max, pm2_min, pm2_max)
Here's what the plot looks like with the bounds we chose.
plot_proper_motion(centerline_df) plt.plot(pm1_rect, pm2_rect, '-');
Now that we've identified the bounds of the cluster in proper motion, we'll use it to select rows from
We'll use the following function, which uses Pandas operators to make a mask that selects rows where
series falls between
def between(series, low, high): """Check whether values are between `low` and `high`.""" return (series > low) & (series < high)
The following mask selects stars with proper motion in the region we chose.
pm1 = results_df['pm_phi1'] pm2 = results_df['pm_phi2'] pm_mask = (between(pm1, pm1_min, pm1_max) & between(pm2, pm2_min, pm2_max))
Again, the sum of a Boolean series is the number of
Now we can use this mask to select rows from
selected_df = results_df[pm_mask] len(selected_df)
These are the stars we think are likely to be in GD-1. Let's see what they look like, plotting their coordinates (not their proper motion).
x = selected_df['phi1'] y = selected_df['phi2'] plt.plot(x, y, 'ko', markersize=1, alpha=1) plt.xlabel('phi1 (degree GD1)') plt.ylabel('phi2 (degree GD1)');
Now that's starting to look like a tidal stream!
At this point we have run a successful query and cleaned up the results; this is a good time to save the data.
To save a Pandas
DataFrame, one option is to convert it to an Astropy
Table, like this:
selected_table = Table.from_pandas(selected_df) type(selected_table)
Then we could write the
Table to a FITS file, as we did in the previous lesson.
But Pandas provides functions to write DataFrames in other formats; to see what they are find the functions here that begin with
One of the best options is HDF5, which is Version 5 of Hierarchical Data Format.
HDF5 is a binary format, so files are small and fast to read and write (like FITS, but unlike XML).
An HDF5 file is similar to an SQL database in the sense that it can contain more than one table, although in HDF5 vocabulary, a table is called a Dataset. (Multi-extension FITS files can also contain more than one table.)
And HDF5 stores the metadata associated with the table, including column names, row labels, and data types (like FITS).
Finally, HDF5 is a cross-language standard, so if you write an HDF5 file with Pandas, you can read it back with many other software tools (more than FITS).
We can write a Pandas
DataFrame to an HDF5 file like this:
filename = 'gd1_data.hdf' selected_df.to_hdf(filename, 'selected_df', mode='w')
Because an HDF5 file can contain more than one Dataset, we have to provide a name, or "key", that identifies the Dataset in the file.
We could use any string as the key, but it will be convenient to give the Dataset in the file the same name as the
mode='w' means that if the file already exists, we should overwrite it.
We're going to need
centerline_df later as well. Write a line of code to add it as a second Dataset in the HDF5 file.
Hint: Since the file already exists, you should not use
# Solution goes here
We can use
getsize to confirm that the file exists and check the size:
from os.path import getsize MB = 1024 * 1024 getsize(filename) / MB
If you forget what the names of the Datasets in the file are, you can read them back like this:
with pd.HDFStore(filename) as hdf: print(hdf.keys())
Python note: We use a
with statement here to open the file before the print statement and (automatically) close it after. Read more about context managers.
The keys are the names of the Datasets. Notice that they start with
/, which indicates that they are at the top level of the Dataset hierarchy, and not in a named "group".
In future lessons we will add a few more Datasets to this file, but not so many that we need to organize them into groups.
In this lesson, we re-loaded the Gaia data we saved from a previous query.
We transformed the coordinates and proper motion from ICRS to a frame aligned with the orbit of GD-1, and stored the results in a Pandas
Then we replicated the selection process from the Price-Whelan and Bonaca paper:
We selected stars near the centerline of GD-1 and made a scatter plot of their proper motion.
We identified a region of proper motion that contains stars likely to be in GD-1.
We used a Boolean
Series as a mask to select stars whose proper motion is in that region.
So far, we have used data from a relatively small region of the sky. In the next lesson, we'll write a query that selects stars based on proper motion, which will allow us to explore a larger region.
When you make a scatter plot, adjust the size of the markers and their transparency so the figure is not overplotted; otherwise it can misrepresent the data badly.
For simple scatter plots in Matplotlib,
plot is faster than
Table and a Pandas
DataFrame are similar in many ways and they provide many of the same functions. They have pros and cons, but for many projects, either one would be a reasonable choice.
To store data from a Pandas
DataFrame, a good option is an HDF file, which can contain multiple Datasets.