+
+
+
+
+
+
+
+
+
+
+
+
+filename = 'gd1_results.fits'
+from astropy.table import Table
+
+results = Table.read(filename)
+
+
+
+
+
+
+from pyia import GaiaData
+
+gaia_data = GaiaData(results)
+type(gaia_data)
+
+
+
+
+pyia.data.GaiaData
+
+
+
+
+
+
+import astropy.units as u
+
+skycoord = gaia_data.get_skycoord(
+ distance=8*u.kpc,
+ radial_velocity=0*u.km/u.s)
+
+
+
+
+
+
+import gala.coordinates as gc
+
+transformed = skycoord.transform_to(gc.GD1Koposov10)
+type(transformed)
+
+
+
+
+astropy.coordinates.sky_coordinate.SkyCoord
+
+
+
+
+
+
+gd1_coord = gc.reflex_correct(transformed)
+
+type(gd1_coord)
+
+
+
+
+astropy.coordinates.sky_coordinate.SkyCoord
+
+
+
+
+
+
+phi1 = gd1_coord.phi1
+phi2 = gd1_coord.phi2
+
+
+
+
+
+
+import matplotlib.pyplot as plt
+
+plt.plot(phi1, phi2, 'ko', markersize=0.1, alpha=0.2)
+
+plt.xlabel('ra (degree GD1)')
+plt.ylabel('dec (degree GD1)');
+
+
+
+
+
+
+
+
+
+
+import pandas as pd
+
+df = results.to_pandas()
+df.shape
+
+
+
+
+(140340, 8)
+
+
+
+
+
+
+df['phi1'] = gd1_coord.phi1
+df['phi2'] = gd1_coord.phi2
+df.shape
+
+
+
+
+(140340, 10)
+
+
+
+
+
+
+df['pm_phi1'] = gd1_coord.pm_phi1_cosphi2
+df['pm_phi2'] = gd1_coord.pm_phi2
+df.shape
+
+
+
+
+(140340, 12)
+
+
+
+
+
+
+phi2_min = -1.0 * u.deg
+phi2_max = 1.0 * u.deg
+
+phi_mask = ((df['phi2'] > phi2_min) &
+ (df['phi2'] < phi2_max))
+
+phi_mask.sum()
+
+
+
+
+25084
+
+
+
+
+
+
+centerline = df[phi_mask]
+
+
+
+
+
+
+pm1 = centerline['pm_phi1']
+pm2 = centerline['pm_phi2']
+
+plt.plot(pm1, pm2, 'ko', markersize=0.3, alpha=0.3)
+
+plt.xlabel('Proper motion phi1 (GD1 frame)')
+plt.ylabel('Proper motion phi2 (GD1 frame)')
+
+plt.xlim(-12, 8)
+plt.ylim(-10, 10);
+
+
+
+
+
+
+
+
+
+
+pm1_min = -8.9
+pm1_max = -6.9
+pm2_min = -2.2
+pm2_max = 1.0
+
+
+
+
+
+
+pm1_rect = [pm1_min, pm1_min, pm1_max, pm1_max, pm1_min] * u.mas/u.yr
+pm2_rect = [pm2_min, pm2_max, pm2_max, pm2_min, pm2_min] * u.mas/u.yr
+
+
+
+
+
+
+plt.plot(pm1, pm2, 'ko', markersize=0.3, alpha=0.3)
+plt.plot(pm1_rect, pm2_rect, '-')
+
+plt.xlabel('Proper motion phi1 (GD1 frame)')
+plt.ylabel('Proper motion phi2 (GD1 frame)')
+
+plt.xlim(-12, 8)
+plt.ylim(-10, 10);
+
+
+
+
+
+
+
+
+
+
+def between(series, low, high):
+ """Make a Boolean Series.
+
+ series: Pandas Series
+ low: lower bound
+ high: upper bound
+
+ returns: Boolean Series
+ """
+ return (series > low) & (series < high)
+
+
+
+
+
+
+pm_mask = (between(df['pm_phi1'], pm1_min, pm1_max) &
+ between(df['pm_phi2'], pm2_min, pm2_max))
+
+
+
+
+
+
+pm_mask.sum()
+
+
+
+
+1049
+
+
+
+
+
+
+selected = df[pm_mask]
+
+
+
+
+
+
+phi1 = selected['phi1']
+phi2 = selected['phi2']
+
+plt.plot(phi1, phi2, 'ko', markersize=0.5, alpha=0.5)
+
+plt.xlabel('ra (degree GD1)')
+plt.ylabel('dec (degree GD1)');
+
+
+
+
+
+
+
+
+
+
+filename = 'gd1_dataframe.hdf5'
+
+df.to_hdf(filename, 'df')
+centerline.to_hdf(filename, 'centerline')
+selected.to_hdf(filename, 'selected')
+
+
+
+
+
+
+!ls -lh gd1_dataframe.hdf5
+
+
+
+
+-rw-rw-r-- 1 downey downey 17M Nov 18 19:06 gd1_dataframe.hdf5
+
+
+
+
+
+
+!rm gd1_dataframe.hdf5
+
+
+
+
+
+
+!wget https://github.com/AllenDowney/AstronomicalData/raw/main/data/gd1_dataframe.hdf5
+
+
+
+
+--2020-11-18 19:31:21-- https://github.com/AllenDowney/AstronomicalData/raw/main/data/gd1_dataframe.hdf5
+Resolving github.com... 140.82.114.3
+Connecting to github.com|140.82.114.3|:443... connected.
+HTTP request sent, awaiting response... 302 Found
+Location: https://raw.githubusercontent.com/AllenDowney/AstronomicalData/main/data/gd1_dataframe.hdf5 [following]
+--2020-11-18 19:31:21-- https://raw.githubusercontent.com/AllenDowney/AstronomicalData/main/data/gd1_dataframe.hdf5
+Resolving raw.githubusercontent.com... 151.101.116.133
+Connecting to raw.githubusercontent.com|151.101.116.133|:443... connected.
+HTTP request sent, awaiting response... 200 OK
+Length: 17335888 (17M) [application/octet-stream]
+Saving to: ‘gd1_dataframe.hdf5’
+
+gd1_dataframe.hdf5 100%[===================>] 16.53M 8.77MB/s in 1.9s
+
+2020-11-18 19:31:24 (8.77 MB/s) - ‘gd1_dataframe.hdf5’ saved [17335888/17335888]
+
+
+
+
+
+
+import pandas as pd
+
+filename = 'gd1_dataframe.hdf5'
+df = pd.read_hdf(filename, 'df')
+centerline = pd.read_hdf(filename, 'centerline')
+selected = pd.read_hdf(filename, 'selected')
+
+
+
+
+
+
+df.shape
+
+
+
+
+(140340, 12)
+
+
+
+
+
+
+centerline.shape
+
+
+
+
+(25084, 12)
+
+
+
+
+
+
+selected.shape
+
+
+
+
+(1049, 12)
+
+
+
+
+
+
+import numpy as np
+
+
+
+
+
+
+points = selected[['pmra', 'pmdec']].to_numpy()
+
+
+
+
+
+
+points
+
+
+
+
+array([[ -3.65859609, -12.48641977],
+ [ -3.57235105, -13.19584836],
+ [ -4.52678061, -12.74505387],
+ ...,
+ [ -3.7313674 , -13.06684162],
+ [ -3.63378135, -12.93358913],
+ [ -3.44742277, -13.22808211]])
+
+
+
+
+
+
+points.shape
+
+
+
+
+(1049, 2)
+
+
+
+
+
+
+from scipy.spatial import ConvexHull
+
+
+
+
+
+
+hull = ConvexHull(points)
+
+
+
+
+
+
+hull
+
+
+
+
+<scipy.spatial.qhull.ConvexHull at 0x7fac07b1b8e0>
+
+
+
+
+
+
+hull.vertices
+
+
+
+
+array([ 692, 873, 141, 303, 42, 622, 45, 83, 127, 182, 1006,
+ 971, 967, 1001, 969, 940], dtype=int32)
+
+
+
+
+
+
+pm_vertices = points[hull.vertices]
+
+
+
+
+
+
+type(pm_vertices)
+
+
+
+
+numpy.ndarray
+
+
+
+
+
+
+pm_vertices.shape
+
+
+
+
+(16, 2)
+
+
+
+
+
+
+pm_vertices
+
+
+
+
+array([[ -4.05037121, -14.75623261],
+ [ -3.41981085, -14.72365546],
+ [ -3.03521988, -14.44357135],
+ [ -2.26847919, -13.7140236 ],
+ [ -2.61172203, -13.24797471],
+ [ -2.73471401, -13.09054471],
+ [ -3.19923146, -12.5942653 ],
+ [ -3.34082546, -12.47611926],
+ [ -5.67489413, -11.16083338],
+ [ -5.95159272, -11.10547884],
+ [ -6.42394023, -11.05981295],
+ [ -7.09631023, -11.95187806],
+ [ -7.30641519, -12.24559977],
+ [ -7.04016696, -12.88580702],
+ [ -6.00347705, -13.75912098],
+ [ -4.42442296, -14.74641176]])
+
+
+
+
+
+
+xs, ys = np.transpose(pm_vertices)
+
+
+
+
+
+
+xs
+
+
+
+
+array([-4.05037121, -3.41981085, -3.03521988, -2.26847919, -2.61172203,
+ -2.73471401, -3.19923146, -3.34082546, -5.67489413, -5.95159272,
+ -6.42394023, -7.09631023, -7.30641519, -7.04016696, -6.00347705,
+ -4.42442296])
+
+
+
+
+
+
+ys
+
+
+
+
+array([-14.75623261, -14.72365546, -14.44357135, -13.7140236 ,
+ -13.24797471, -13.09054471, -12.5942653 , -12.47611926,
+ -11.16083338, -11.10547884, -11.05981295, -11.95187806,
+ -12.24559977, -12.88580702, -13.75912098, -14.74641176])
+
+
+
+
+
+
+pm1 = centerline['pmra']
+pm2 = centerline['pmdec']
+plt.plot(pm1, pm2, 'ko', markersize=0.3, alpha=0.3)
+
+pm1 = selected['pmra']
+pm2 = selected['pmdec']
+plt.plot(pm1, pm2, 'gx', markersize=1, alpha=0.3)
+
+plt.plot(xs, ys)
+
+plt.xlabel('Proper motion ra (ICRS frame)')
+plt.ylabel('Proper motion dec (ICRS frame)')
+
+plt.xlim([-10, 5])
+plt.ylim([-20, 5]);
+
+
+
+
+
+
+
+
+
+
+pm_vertices
+
+
+
+
+array([[ -4.05037121, -14.75623261],
+ [ -3.41981085, -14.72365546],
+ [ -3.03521988, -14.44357135],
+ [ -2.26847919, -13.7140236 ],
+ [ -2.61172203, -13.24797471],
+ [ -2.73471401, -13.09054471],
+ [ -3.19923146, -12.5942653 ],
+ [ -3.34082546, -12.47611926],
+ [ -5.67489413, -11.16083338],
+ [ -5.95159272, -11.10547884],
+ [ -6.42394023, -11.05981295],
+ [ -7.09631023, -11.95187806],
+ [ -7.30641519, -12.24559977],
+ [ -7.04016696, -12.88580702],
+ [ -6.00347705, -13.75912098],
+ [ -4.42442296, -14.74641176]])
+
+
+
+
+
+
+t = [str(x) for x in pm_vertices.flatten()]
+len(t)
+
+
+
+
+32
+
+
+
+
+
+
+pm_point_list = ', '.join(t)
+pm_point_list
+
+
+
+
+'-4.050371212154984, -14.75623260987968, -3.4198108491382455, -14.723655456335619, -3.035219883740934, -14.443571352854612, -2.268479190206636, -13.714023598831554, -2.611722027231764, -13.247974712069263, -2.7347140078529106, -13.090544709622938, -3.199231461993783, -12.594265302440828, -3.34082545787549, -12.476119260818695, -5.674894125178565, -11.160833381392624, -5.95159272432137, -11.105478836426514, -6.423940229776128, -11.05981294804957, -7.096310230579248, -11.951878058650085, -7.306415190921692, -12.245599765990594, -7.040166963232815, -12.885807024935527, -6.0034770546523735, -13.759120984106968, -4.42442296194263, -14.7464117578883'
+
+
+
+
+
+
+phi1_min = -70
+phi1_max = -20
+phi2_min = -5
+phi2_max = 5
+
+
+
+
+
+
+import astropy.units as u
+
+phi1_rect = [phi1_min, phi1_min, phi1_max, phi1_max] * u.deg
+phi2_rect = [phi2_min, phi2_max, phi2_max, phi2_min] * u.deg
+
+
+
+
+
+
+import gala.coordinates as gc
+import astropy.coordinates as coord
+
+corners = gc.GD1Koposov10(phi1=phi1_rect, phi2=phi2_rect)
+corners_icrs = corners.transform_to(coord.ICRS)
+
+
+
+
+
+
+point_base = "{point.ra.value}, {point.dec.value}"
+
+t = [point_base.format(point=point)
+ for point in corners_icrs]
+
+point_list = ', '.join(t)
+point_list
+
+
+
+
+'135.30559858565638, 8.398623940157561, 126.50951508623503, 13.44494195652069, 163.0173655836748, 54.24242734020255, 172.9328536286811, 46.47260492416258'
+
+
+
+
+
+
+query_base = """SELECT
+{columns}
+FROM gaiadr2.gaia_source
+WHERE parallax < 1
+ AND bp_rp BETWEEN -0.75 AND 2
+ AND 1 = CONTAINS(POINT(ra, dec),
+ POLYGON({point_list}))
+"""
+
+
+
+
+
+
+columns = 'source_id, ra, dec, pmra, pmdec, parallax, parallax_error, radial_velocity'
+Exercise 1: Use format to format query_base and define query, filling in the values of columns and point_list.
+
+
+
+
+
+query = query_base.format(columns=columns, point_list=point_list)
+
+
+
+
+
+
+print(query)
+
+
+
+
+SELECT
+source_id, ra, dec, pmra, pmdec, parallax, parallax_error, radial_velocity
+FROM gaiadr2.gaia_source
+WHERE parallax < 1
+ AND bp_rp BETWEEN -0.75 AND 2
+ AND 1 = CONTAINS(POINT(ra, dec),
+ POLYGON(135.30559858565638, 8.398623940157561, 126.50951508623503, 13.44494195652069, 163.0173655836748, 54.24242734020255, 172.9328536286811, 46.47260492416258))
+Exercise 2: Modify query_base by adding a new clause to select stars whose coordinates of proper motion, pmra and pmdec, fall within the polygon defined by pm_point_list.
+
+
+
+
+
+query_base = """SELECT
+{columns}
+FROM gaiadr2.gaia_source
+WHERE parallax < 1
+ AND bp_rp BETWEEN -0.75 AND 2
+ AND 1 = CONTAINS(POINT(ra, dec),
+ POLYGON({point_list}))
+ AND 1 = CONTAINS(POINT(pmra, pmdec),
+ POLYGON({pm_point_list}))
+"""
+
+
+
+
+
+
+query = query_base.format(columns=columns,
+ point_list=point_list,
+ pm_point_list=pm_point_list)
+
+
+
+
+
+
+print(query)
+
+
+
+
+SELECT
+source_id, ra, dec, pmra, pmdec, parallax, parallax_error, radial_velocity
+FROM gaiadr2.gaia_source
+WHERE parallax < 1
+ AND bp_rp BETWEEN -0.75 AND 2
+ AND 1 = CONTAINS(POINT(ra, dec),
+ POLYGON(135.30559858565638, 8.398623940157561, 126.50951508623503, 13.44494195652069, 163.0173655836748, 54.24242734020255, 172.9328536286811, 46.47260492416258))
+ AND 1 = CONTAINS(POINT(pmra, pmdec),
+ POLYGON(-4.050371212154984, -14.75623260987968, -3.4198108491382455, -14.723655456335619, -3.035219883740934, -14.443571352854612, -2.268479190206636, -13.714023598831554, -2.611722027231764, -13.247974712069263, -2.7347140078529106, -13.090544709622938, -3.199231461993783, -12.594265302440828, -3.34082545787549, -12.476119260818695, -5.674894125178565, -11.160833381392624, -5.95159272432137, -11.105478836426514, -6.423940229776128, -11.05981294804957, -7.096310230579248, -11.951878058650085, -7.306415190921692, -12.245599765990594, -7.040166963232815, -12.885807024935527, -6.0034770546523735, -13.759120984106968, -4.42442296194263, -14.7464117578883))
+
+
+
+
+
+
+from astroquery.gaia import Gaia
+
+job = Gaia.launch_job_async(query)
+print(job)
+
+
+
+
+Created TAP+ (v1.2.1) - Connection:
+ Host: gea.esac.esa.int
+ Use HTTPS: True
+ Port: 443
+ SSL Port: 443
+Created TAP+ (v1.2.1) - Connection:
+ Host: geadata.esac.esa.int
+ Use HTTPS: True
+ Port: 443
+ SSL Port: 443
+INFO: Query finished. [astroquery.utils.tap.core]
+<Table length=7346>
+ name dtype unit description n_bad
+--------------- ------- -------- ------------------------------------------------------------------ -----
+ source_id int64 Unique source identifier (unique within a particular Data Release) 0
+ ra float64 deg Right ascension 0
+ dec float64 deg Declination 0
+ pmra float64 mas / yr Proper motion in right ascension direction 0
+ pmdec float64 mas / yr Proper motion in declination direction 0
+ parallax float64 mas Parallax 0
+ parallax_error float64 mas Standard error of parallax 0
+radial_velocity float64 km / s Radial velocity 7295
+Jobid: 1605731697505O
+Phase: COMPLETED
+Owner: None
+Output file: async_20201118153457.vot
+Results: None
+
+
+
+
+
+
+candidate_table = job.get_results()
+type(candidate_table)
+
+
+
+
+astropy.table.table.Table
+
+
+
+
+
+
+x = candidate_table['ra']
+y = candidate_table['dec']
+
+plt.plot(x, y, 'ko', markersize=0.3, alpha=0.3)
+
+plt.xlabel('ra (deg ICRS)')
+plt.ylabel('dec (deg ICRS)');
+
+
+
+
+
+
+
+
+
+
+from pyia import GaiaData
+
+def make_dataframe(table):
+ """Transform coordinates from ICRS to GD-1 frame.
+
+ table: Astropy Table
+
+ returns: Pandas DataFrame
+ """
+ gaia_data = GaiaData(table)
+
+ c_sky = gaia_data.get_skycoord(distance=8*u.kpc,
+ radial_velocity=0*u.km/u.s)
+ c_gd1 = gc.reflex_correct(
+ c_sky.transform_to(gc.GD1Koposov10))
+
+ df = table.to_pandas()
+ df['phi1'] = c_gd1.phi1
+ df['phi2'] = c_gd1.phi2
+ df['pm_phi1'] = c_gd1.pm_phi1_cosphi2
+ df['pm_phi2'] = c_gd1.pm_phi2
+ return df
+
+
+
+
+
+
+candidate_df = make_dataframe(candidate_table)
+candidate_df.shape
+
+
+
+
+(7346, 12)
+
+
+
+
+
+
+x = candidate_df['phi1']
+y = candidate_df['phi2']
+
+plt.plot(x, y, 'ko', markersize=0.3, alpha=0.3)
+
+plt.xlabel('phi1 (deg GD1)')
+plt.ylabel('phi2 (deg GD1)');
+
+
+
+
+Exercise Write candidate_df to an HDF file named gd1_candidates.
Save it as a dataset named candidate_df.
+
+
+
+
+
+candidate_df.to_hdf('gd1_candidates.hdf', 'candidate_df')
+
+
+
+
+
+
+!ls -lh gd1_candidates.hdf
+
+
+
+
+-rw-rw-r-- 1 downey downey 756K Nov 18 15:54 gd1_candidates.hdf
+
+
+
+
+
+
+candidate_df.to_csv('gd1_candidates.csv')
+
+
+
+
+
+
+!ls -lh gd1_candidates.csv
+
+
+
+
+-rw-rw-r-- 1 downey downey 1.6M Nov 18 15:51 gd1_candidates.csv
+
+
+
+
+
+
+
+
+
+
+!head gd1_candidates.csv
+
+
+
+
+,source_id,ra,dec,pmra,pmdec,parallax,parallax_error,radial_velocity,phi1,phi2,pm_phi1,pm_phi2
+0,635559124339440000,137.58671691646745,19.1965441084838,-3.770521900009566,-12.490481778113859,0.7913934419894347,0.2717538145759051,,-59.63048941944396,-1.21648525150429,-7.3613627125974945,-0.5926328820644919
+1,635860218726658176,138.5187065217173,19.09233926905897,-5.941679495793577,-11.346409129876392,0.30745551377348623,0.19946557779138105,,-59.247329893833296,-2.0160784008206476,-7.52712608464053,1.7487794924176672
+2,635674126383965568,138.8428741026386,19.031798198627634,-3.8970011609340207,-12.702779525389634,0.7794631695848825,0.22369151932421705,,-59.13339098769217,-2.3069007451798287,-7.560607655557421,-0.7417999555980245
+3,635535454774983040,137.8377518255436,18.864006786112604,-4.335040664412791,-14.492308604905652,0.31451429349160354,0.10277456217796248,,-59.78529963350303,-1.5945692849622026,-9.35753561702402,-1.2184922339592716
+4,635497276810313600,138.0445160213759,19.00947118796605,-7.1729306406216615,-12.291499169815987,0.42540354480291037,0.3376892347379532,,-59.55774434692722,-1.6821467937242753,-9.000831056712947,2.3344066165829673
+5,635614168640132864,139.59219748145836,18.807955539071433,-3.309602916796381,-13.708904908478631,0.5833821573394844,0.5746175895352217,,-58.938113464510934,-3.0241923218278397,-8.062761531706451,-1.8690819375362338
+6,635821843194387840,139.88094034815086,19.62185456718988,-6.544201177153814,-12.55978220563274,0.3788707393484704,0.2134264091820104,,-58.10569437830526,-2.8105428528649257,-8.814581159026405,1.451032881430119
+7,635551706931167104,138.04665586038192,19.248909662830798,-6.224595114220405,-12.224246333795001,0.7832805254544403,0.12699359319614956,,-59.354285349557074,-1.5557649698784075,-8.427716582116249,1.5738026079151446
+8,635518889086133376,137.2374229207837,18.7428630711791,-3.3186800714801046,-12.710314902969365,0.17932476739336733,0.7503378631682915,,-60.19048364007102,-1.1779749872253444,-7.3401827273605385,-1.064062501817137
+