Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 9d88fc8

Browse files
luthienliuluthienliu
committedAug 6, 2022
Change to 4_Coordinates_Crossmatch.ipynb
1 parent d1dfb7c commit 9d88fc8

File tree

2 files changed

+557
-4
lines changed

2 files changed

+557
-4
lines changed
 

‎tutorials/astropy-coordinates/4_Coordinates_Crossmatch.ipynb

Lines changed: 557 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,557 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {
6+
"id": "pvmxjgnDopS6"
7+
},
8+
"source": [
9+
"# Astronomical Coordinates 4: Cross-matching Catalogs Using astropy.coordinates and astroquery\n",
10+
"\n",
11+
"## Authors\n",
12+
"Adrian Price-Whelan, Lúthien Liu, Zihao Chen, Saima Siddiqui\n",
13+
"\n",
14+
"## Learning Goals\n",
15+
"* Demonstrate how to retrieve a catalog from Vizier using astroquery\n",
16+
"* Show how to perform positional cross-matches between catalogs of sky coordinates\n",
17+
"\n",
18+
"## Keywords\n",
19+
"coordinates, OOP, astroquery, gaia\n",
20+
"\n",
21+
"\n",
22+
"## Summary\n",
23+
"\n",
24+
"In the previous tutorials in this series, we introduced many of the key concepts underlying how to represent and transform astronomical coordinates using `astropy.coordinates`, including how to work with both position and velocity data within the coordinate objects.\n",
25+
"\n",
26+
"In this tutorial, we will explore how the `astropy.coordinates` package can be used to cross-match two catalogs that contain overlapping sources that may have been observed at different times. You may find it helpful to keep [the Astropy documentation for the coordinates package](http://docs.astropy.org/en/stable/coordinates/index.html) open alongside this tutorial for reference or additional reading. In the text below, you may also see some links that look like ([docs](http://docs.astropy.org/en/stable/coordinates/index.html)). These links will take you to parts of the documentation that are directly relevant to the cells from which they link. \n",
27+
"\n",
28+
"*Note: This is the 4th tutorial in a series of tutorials about astropy.coordinates. If you are new to astropy.coordinates, you may want to start from the beginning or an earlier tutorial.*\n",
29+
"- [Previous tutorial: Astronomical Coordinates 3: Working with Velocity Data](3-Coordinates-Velocities)"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {
35+
"id": "sKhydue_opS-"
36+
},
37+
"source": [
38+
"## Imports\n",
39+
"\n",
40+
"We start by importing some general packages that we will need below:"
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": null,
46+
"metadata": {
47+
"colab": {
48+
"base_uri": "https://localhost:8080/"
49+
},
50+
"id": "741ulTSBopS_",
51+
"outputId": "1fea5dc8-ffd2-4aae-b9b6-d5c3ea534a2f"
52+
},
53+
"outputs": [],
54+
"source": [
55+
"import warnings\n",
56+
"\n",
57+
"import matplotlib.pyplot as plt\n",
58+
"%matplotlib inline\n",
59+
"import numpy as np\n",
60+
"\n",
61+
"from astropy import units as u\n",
62+
"from astropy.coordinates import SkyCoord, Distance\n",
63+
"from astropy.io import fits\n",
64+
"from astropy.table import QTable\n",
65+
"from astropy.time import Time\n",
66+
"from astropy.utils.data import download_file\n",
67+
"\n",
68+
"from astroquery.vizier import Vizier"
69+
]
70+
},
71+
{
72+
"cell_type": "markdown",
73+
"metadata": {
74+
"id": "ZKJjNwxiopTB"
75+
},
76+
"source": [
77+
"## Cross-matching and comparing catalogs\n",
78+
"\n",
79+
"In this tutorial, we are going to return to a set of data that we downloaded from the *Gaia* archive back in [Tutorial 1](1-Coordinates-Intro) of this series.\n",
80+
"\n",
81+
"Let's recap what we did in that tutorial: We defined a `SkyCoord` object to represent the center of an open cluster (NGC 188), we queried the *Gaia* DR2 catalog to select stars that are close (on the sky) to the center of the cluster, and we used the parallax values from *Gaia* to select stars that are near NGC 188 in 3D position. Here, we will briefly reproduce those selections so that we can start here with a catalog of sources that are likely members of NGC 188 (see [Tutorial 1](1-Coordinates-Intro) for more information):"
82+
]
83+
},
84+
{
85+
"cell_type": "code",
86+
"execution_count": null,
87+
"metadata": {
88+
"colab": {
89+
"base_uri": "https://localhost:8080/"
90+
},
91+
"id": "pEWYbW8UopTC",
92+
"outputId": "0fed4911-59bc-4980-d4f2-9319ffe76dec"
93+
},
94+
"outputs": [],
95+
"source": [
96+
"ngc188_table = QTable.read('gaia_results.fits')\n",
97+
"ngc188_table = ngc188_table[ngc188_table['parallax'] > 0.25*u.mas]\n",
98+
"\n",
99+
"ngc188_center_3d = SkyCoord(12.11*u.deg, 85.26*u.deg,\n",
100+
" distance=1.96*u.kpc,\n",
101+
" pm_ra_cosdec=-2.3087*u.mas/u.yr,\n",
102+
" pm_dec=-0.9565*u.mas/u.yr)\n",
103+
"\n",
104+
"# Deal with masked quantity data in a backwards-compatible way:\n",
105+
"parallax = ngc188_table['parallax']\n",
106+
"if hasattr(parallax, 'mask'):\n",
107+
" parallax = parallax.filled(np.nan)\n",
108+
" \n",
109+
"velocity_data = {\n",
110+
" 'pm_ra_cosdec': ngc188_table['pmra'],\n",
111+
" 'pm_dec': ngc188_table['pmdec'],\n",
112+
" 'radial_velocity': ngc188_table['radial_velocity']\n",
113+
"}\n",
114+
"for k, v in velocity_data.items():\n",
115+
" if hasattr(v, 'mask'):\n",
116+
" velocity_data[k] = v.filled(0.)\n",
117+
" velocity_data[k][np.isnan(velocity_data[k])] = 0.\n",
118+
"\n",
119+
"ngc188_coords_3d = SkyCoord(\n",
120+
" ra=ngc188_table['ra'], \n",
121+
" dec=ngc188_table['dec'],\n",
122+
" distance=Distance(parallax=parallax),\n",
123+
" obstime=Time('J2015.5'),\n",
124+
" **velocity_data\n",
125+
")\n",
126+
"\n",
127+
"sep3d = ngc188_coords_3d.separation_3d(ngc188_center_3d)\n",
128+
"pm_diff = np.sqrt(\n",
129+
" (ngc188_table['pmra'] - ngc188_center_3d.pm_ra_cosdec)**2 + \n",
130+
" (ngc188_table['pmdec'] - ngc188_center_3d.pm_dec)**2)\n",
131+
"\n",
132+
"ngc188_members_mask = (sep3d < 50*u.pc) & (pm_diff < 1.5*u.mas/u.yr)\n",
133+
"ngc188_members = ngc188_table[ngc188_members_mask]\n",
134+
"ngc188_members_coords = ngc188_coords_3d[ngc188_members_mask]\n",
135+
"len(ngc188_members)"
136+
]
137+
},
138+
{
139+
"cell_type": "markdown",
140+
"metadata": {
141+
"id": "HH-OAdGiopTD"
142+
},
143+
"source": [
144+
"From the selections above, the table `ngc188_members` and the `SkyCoord` instance `ngc188_members_coords` contain 216 sources that, based on their 3D positions and proper motions, are consistent with being members of the open cluster NGC 188."
145+
]
146+
},
147+
{
148+
"cell_type": "markdown",
149+
"metadata": {
150+
"id": "j3vx-xpzopTE"
151+
},
152+
"source": [
153+
"Let's assume that we now want to cross-match our catalog of candidate members of NGC 188 — here, based on *Gaia* data — to some other catalog. In this tutorial, we will demonstrate how to manually cross-match these *Gaia* sources with the 2MASS photometric catalog to retrieve infrared magnitudes for these stars, and then we will plot a color–magnitude diagram. To do this, we first need to query the 2MASS catalog to retrieve all sources in a region around the center of NGC 188, as we did for *Gaia*. Here, we will also take into account the fact that the *Gaia* data release 2 reference epoch is J2015.5, whereas the 2MASS coordinates are likely reported at their time of observation (in the late 1990's). \n",
154+
"\n",
155+
"*Note that some data archives, like the Gaia science archive, support running cross-matches at the database level and even support epoch propagation. If you need to perform a large cross-match, it will be much more efficient to use these services!*\n",
156+
"\n",
157+
"We will again use `astroquery` to execute this query. This will again require an internet connection, but we have included the results of this query in a file along with this notebook in case you are not connected to the internet. To query 2MASS, we will use the `astroquery.vizier` module ([docs](https://astroquery.readthedocs.io/en/latest/vizier/vizier.html)) to run a cone search centered on the sky position of NGC 188 with a search radius of 0.5º:"
158+
]
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": null,
163+
"metadata": {
164+
"id": "nKlO1UBmopTF"
165+
},
166+
"outputs": [],
167+
"source": [
168+
"# NOTE: skip this cell if you do not have an internet connection\n",
169+
"\n",
170+
"# II/246 is the catalog name for the main 2MASS photometric catalog\n",
171+
"v = Vizier(catalog=\"II/246\", columns=['*', 'Date']) \n",
172+
"v.ROW_LIMIT = -1\n",
173+
"\n",
174+
"result = v.query_region(ngc188_center_3d, radius=0.5*u.deg)\n",
175+
"tmass_table = result[0]"
176+
]
177+
},
178+
{
179+
"cell_type": "markdown",
180+
"metadata": {
181+
"id": "et4ax7AqopTF"
182+
},
183+
"source": [
184+
"Alternatively, we can read the 2MASS table provided along with this tutorial:"
185+
]
186+
},
187+
{
188+
"cell_type": "code",
189+
"execution_count": null,
190+
"metadata": {
191+
"id": "_KNyVfsjopTG"
192+
},
193+
"outputs": [],
194+
"source": [
195+
"# the .read() below produces some warnings that we can safely ignore\n",
196+
"with warnings.catch_warnings(): \n",
197+
" warnings.simplefilter('ignore', UserWarning)\n",
198+
" \n",
199+
" tmass_table = QTable.read('2MASS_results.ecsv')"
200+
]
201+
},
202+
{
203+
"cell_type": "markdown",
204+
"metadata": {
205+
"id": "gH3GmuxEopTH"
206+
},
207+
"source": [
208+
"As with the *Gaia* results table, we can now create a single `SkyCoord` object to represent all of the sources returned from our query to the 2MASS catalog. Let's look at the column names in this table by displaying the first few rows:"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"metadata": {
215+
"colab": {
216+
"base_uri": "https://localhost:8080/",
217+
"height": 165
218+
},
219+
"id": "rfe_mIhGopTI",
220+
"outputId": "50270998-1eb3-459d-89c4-22da8078b643"
221+
},
222+
"outputs": [],
223+
"source": [
224+
"tmass_table[:3]"
225+
]
226+
},
227+
{
228+
"cell_type": "markdown",
229+
"metadata": {
230+
"id": "m3bAC2kCopTI"
231+
},
232+
"source": [
233+
"From looking at the column names, the two relevant sky coordinate columns are `RAJ2000` for `ra` and `DEJ2000` for `dec`:"
234+
]
235+
},
236+
{
237+
"cell_type": "code",
238+
"execution_count": null,
239+
"metadata": {
240+
"colab": {
241+
"base_uri": "https://localhost:8080/"
242+
},
243+
"id": "LWFl1JM9opTJ",
244+
"outputId": "d201a13c-7659-4262-9b75-b6bc26f0a19c"
245+
},
246+
"outputs": [],
247+
"source": [
248+
"tmass_coords = SkyCoord(tmass_table['RAJ2000'], \n",
249+
" tmass_table['DEJ2000'])\n",
250+
"len(tmass_coords)"
251+
]
252+
},
253+
{
254+
"cell_type": "markdown",
255+
"metadata": {
256+
"id": "cWGVnrWRopTJ"
257+
},
258+
"source": [
259+
"Note also that the table contains a \"Date\" column that specifies the epoch of the coordinates. Are all of these epochs the same?"
260+
]
261+
},
262+
{
263+
"cell_type": "code",
264+
"execution_count": null,
265+
"metadata": {
266+
"colab": {
267+
"base_uri": "https://localhost:8080/",
268+
"height": 58
269+
},
270+
"id": "djfwkEGbopTK",
271+
"outputId": "9ab50944-fc07-49e4-9fb2-d7e3f359de7c"
272+
},
273+
"outputs": [],
274+
"source": [
275+
"np.unique(tmass_table['Date'])"
276+
]
277+
},
278+
{
279+
"cell_type": "markdown",
280+
"metadata": {
281+
"id": "2HVHIOzPopTK"
282+
},
283+
"source": [
284+
"It looks like all of the sources in our 2MASS table have the same epoch, so let's create an `astropy.time.Time` object to represent this date:"
285+
]
286+
},
287+
{
288+
"cell_type": "code",
289+
"execution_count": null,
290+
"metadata": {
291+
"id": "1KeO3N9AopTK"
292+
},
293+
"outputs": [],
294+
"source": [
295+
"tmass_epoch = Time(np.unique(tmass_table['Date']))"
296+
]
297+
},
298+
{
299+
"cell_type": "markdown",
300+
"metadata": {
301+
"id": "BTV_a1GlopTL"
302+
},
303+
"source": [
304+
"We now want to cross-match our *Gaia*-selected candidate members of NGC 188, `ngc_members_coords`, with this table of photometry from 2MASS. However, as noted previously, the *Gaia* coordinates are given at a different epoch J2015.5, which is nearly ~16 years after the 2MASS epoch of the data we downloaded (1999-10-19 or roughly J1999.88). We will therefore first use the `SkyCoord.apply_space_motion()` method ([docs](http://docs.astropy.org/en/latest/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord.apply_space_motion)) to transform the *Gaia* positions back to the 2MASS epoch before we do the cross-match:"
305+
]
306+
},
307+
{
308+
"cell_type": "code",
309+
"execution_count": null,
310+
"metadata": {
311+
"id": "Krr18s9gopTM"
312+
},
313+
"outputs": [],
314+
"source": [
315+
"# you can ignore the warning raised here\n",
316+
"ngc188_members_coords_1999 = ngc188_members_coords.apply_space_motion(\n",
317+
" tmass_epoch)"
318+
]
319+
},
320+
{
321+
"cell_type": "markdown",
322+
"metadata": {
323+
"id": "6RVxtRLdopTM"
324+
},
325+
"source": [
326+
"The object `ngc188_members_coords_1999` now contains the coordinate information for our *Gaia*-selected members of NGC 188, as we think they would appear if observed on 1999-10-19.\n",
327+
"\n",
328+
"We can now use the ``SkyCoord.match_to_catalog_sky`` method to match these two catalogs ([docs](http://docs.astropy.org/en/latest/coordinates/matchsep.html#astropy-coordinates-matching)), using the `ngc188_members_coords_1999` as our NGC 188 members coordinates. \n",
329+
"\n",
330+
"Note that order matters with this method: Here we will match *Gaia* to 2MASS. `SkyCoord.match_to_catalog_sky` returns three objects: (1) the indices into `tmass_coords` that get the closest matches in `ngc188_members_coords_1999`, (2) the angular separation between each `ngc188_members_coords_1999` coordinate and the closest source in `tmass_coords`, and (3) the 3D distance between each `ngc188_members_coords_1999` coordinate and the closest source in `tmass_coords`. Here, the 3D distances will not be useful because the 2MASS coordinates do not have associated distance information, so we will ignore these quantities:"
331+
]
332+
},
333+
{
334+
"cell_type": "code",
335+
"execution_count": null,
336+
"metadata": {
337+
"id": "UZHtgJ00opTM"
338+
},
339+
"outputs": [],
340+
"source": [
341+
"idx_gaia, sep2d_gaia, _ = ngc188_members_coords_1999.match_to_catalog_sky(\n",
342+
" tmass_coords)"
343+
]
344+
},
345+
{
346+
"cell_type": "markdown",
347+
"metadata": {
348+
"id": "t1qUzhHWopTN"
349+
},
350+
"source": [
351+
"Let's now look at the distribution of separations (in arcseconds) for all of the cross-matched sources:"
352+
]
353+
},
354+
{
355+
"cell_type": "code",
356+
"execution_count": null,
357+
"metadata": {
358+
"colab": {
359+
"base_uri": "https://localhost:8080/",
360+
"height": 297
361+
},
362+
"id": "Acm50SzTopTN",
363+
"outputId": "490adeed-9d80-4720-9576-c9293aafd15b"
364+
},
365+
"outputs": [],
366+
"source": [
367+
"plt.hist(sep2d_gaia.arcsec, histtype='step', \n",
368+
" bins=np.logspace(-2, 2., 64))\n",
369+
"plt.xlabel('separation [arcsec]')\n",
370+
"plt.xscale('log')\n",
371+
"plt.yscale('log')\n",
372+
"plt.tight_layout()"
373+
]
374+
},
375+
{
376+
"cell_type": "markdown",
377+
"metadata": {
378+
"id": "SG0nUr7HopTN"
379+
},
380+
"source": [
381+
"From this, it looks like all of sources in our *Gaia* NGC 188 member list cross-match to another sources within a few arcseconds, so these all seem like they are correctly matches to a 2MASS source!"
382+
]
383+
},
384+
{
385+
"cell_type": "code",
386+
"execution_count": null,
387+
"metadata": {
388+
"colab": {
389+
"base_uri": "https://localhost:8080/"
390+
},
391+
"id": "5bk543X9opTN",
392+
"outputId": "5af58652-f62b-41fc-9dbf-172dc0a76149"
393+
},
394+
"outputs": [],
395+
"source": [
396+
"(sep2d_gaia < 2*u.arcsec).sum(), len(ngc188_members)"
397+
]
398+
},
399+
{
400+
"cell_type": "markdown",
401+
"metadata": {
402+
"id": "kXdqAdfnopTO"
403+
},
404+
"source": [
405+
"With our cross-match done, we can now make `Gaia`+2MASS color–magnitude diagrams of our candidate NGC 188 members using the information returned by the cross-match:"
406+
]
407+
},
408+
{
409+
"cell_type": "code",
410+
"execution_count": null,
411+
"metadata": {
412+
"id": "OHI0B-TYopTO"
413+
},
414+
"outputs": [],
415+
"source": [
416+
"Jmag = tmass_table['Jmag'][idx_gaia] # note that we use the index array returned above\n",
417+
"Gmag = ngc188_members['phot_g_mean_mag']\n",
418+
"Bmag = ngc188_members['phot_bp_mean_mag']"
419+
]
420+
},
421+
{
422+
"cell_type": "code",
423+
"execution_count": null,
424+
"metadata": {
425+
"colab": {
426+
"base_uri": "https://localhost:8080/",
427+
"height": 369
428+
},
429+
"id": "fuePUx8uopTO",
430+
"outputId": "46515b44-4c91-4615-b830-3211c9252843"
431+
},
432+
"outputs": [],
433+
"source": [
434+
"fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
435+
"\n",
436+
"ax = axes[0]\n",
437+
"ax.scatter(Gmag - Jmag, Gmag, \n",
438+
" marker='o', color='k', \n",
439+
" linewidth=0, alpha=0.5)\n",
440+
"ax.set_xlabel('$G - J$')\n",
441+
"ax.set_ylabel('$G$')\n",
442+
"ax.set_xlim(0, 3)\n",
443+
"ax.set_ylim(19, 10) # backwards because magnitudes!\n",
444+
"\n",
445+
"ax = axes[1]\n",
446+
"ax.scatter(Bmag - Gmag, Jmag, \n",
447+
" marker='o', color='k', \n",
448+
" linewidth=0, alpha=0.5)\n",
449+
"ax.set_xlabel('$G_{BP} - G$')\n",
450+
"ax.set_ylabel('$J$')\n",
451+
"ax.set_xlim(0.2, 1)\n",
452+
"ax.set_ylim(17, 8) # backwards because magnitudes!\n",
453+
"\n",
454+
"fig.tight_layout()"
455+
]
456+
},
457+
{
458+
"cell_type": "markdown",
459+
"metadata": {
460+
"id": "01UZ3GFVopTP"
461+
},
462+
"source": [
463+
"Those both look like color-magnitude diagrams of a main sequence + red giant branch of an intermediate-age stellar cluster, so it looks like our selection and cross-matching has worked!\n",
464+
"\n",
465+
"For more on what matching options are available, check out the [separation and matching section of the Astropy documentation](https://astropy.readthedocs.io/en/stable/coordinates/matchsep.html). Or for more on what you can do with `SkyCoord`, see [its API documentation](http://astropy.readthedocs.org/en/stable/api/astropy.coordinates.SkyCoord.html)."
466+
]
467+
},
468+
{
469+
"cell_type": "markdown",
470+
"metadata": {
471+
"id": "95uKx2_E-Lfd"
472+
},
473+
"source": [
474+
"## Exercises"
475+
]
476+
},
477+
{
478+
"cell_type": "markdown",
479+
"metadata": {
480+
"id": "J4TWmYPf-O-z"
481+
},
482+
"source": [
483+
"Instead of transforming the Gaia positions back to the 2MASS epoch for cross-matching, try to move the 2MASS positions to the Gaia's epoch (J2015.5), and name it as `tmass_coords_2015`."
484+
]
485+
},
486+
{
487+
"cell_type": "code",
488+
"execution_count": null,
489+
"metadata": {
490+
"id": "EfNzCmP1opTP"
491+
},
492+
"outputs": [],
493+
"source": []
494+
},
495+
{
496+
"cell_type": "markdown",
497+
"metadata": {
498+
"id": "pZyezlr--Uyo"
499+
},
500+
"source": [
501+
"Now apply the `SkyCoord.match_to_catalog_sky` method to do a new cross-matching between these 2 catalogs, using the `tmass_coords_2015` as our NGC 188 members coordinates this time. Plot a histogram to show the distribution of separation. Is there any difference between this trial and the result above?"
502+
]
503+
},
504+
{
505+
"cell_type": "code",
506+
"execution_count": null,
507+
"metadata": {
508+
"id": "Th2k_vW_-XiP"
509+
},
510+
"outputs": [],
511+
"source": []
512+
},
513+
{
514+
"cell_type": "markdown",
515+
"metadata": {
516+
"id": "a0pCmK01-349"
517+
},
518+
"source": [
519+
"Instead of picking all stars within 50 pc of the cluster center, try to enlarge to a larger range of 80 pc for cross-match and make a comparison between the result above."
520+
]
521+
},
522+
{
523+
"cell_type": "code",
524+
"execution_count": null,
525+
"metadata": {
526+
"id": "gCKRdYix-704"
527+
},
528+
"outputs": [],
529+
"source": []
530+
}
531+
],
532+
"metadata": {
533+
"colab": {
534+
"name": "4-Coordinates-Crossmatch.ipynb",
535+
"provenance": []
536+
},
537+
"kernelspec": {
538+
"display_name": "Python 3 (ipykernel)",
539+
"language": "python",
540+
"name": "python3"
541+
},
542+
"language_info": {
543+
"codemirror_mode": {
544+
"name": "ipython",
545+
"version": 3
546+
},
547+
"file_extension": ".py",
548+
"mimetype": "text/x-python",
549+
"name": "python",
550+
"nbconvert_exporter": "python",
551+
"pygments_lexer": "ipython3",
552+
"version": "3.9.12"
553+
}
554+
},
555+
"nbformat": 4,
556+
"nbformat_minor": 1
557+
}

‎tutorials/astropy-coordinates/requirements.txt

Lines changed: 0 additions & 4 deletions
This file was deleted.

0 commit comments

Comments
 (0)
Please sign in to comment.