-
Notifications
You must be signed in to change notification settings - Fork 7
/
get_galaxies.py
291 lines (248 loc) · 9.39 KB
/
get_galaxies.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
# Author: Igor Andreoni
# email: [email protected]
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from astroquery.vizier import Vizier
from astropy.table import Table, vstack
from astropy.cosmology import Planck15 as cosmo
def query_glade(ra, dec, rad=1, dist_min=0, dist_max=np.inf,
sep_max_kpc=np.inf, catalog='VII/281/glade2'):
"""
Fetch galaxies from the GLADE catalog
---
Parameters
ra float
Right ascension (deg)
dec float
Declination (deg)
rad float
Search radius (deg)
dist_min float
Min distance of the galaxies (Mpc)
dist_max float
Max distance of the galaxies (Mpc)
sep_max_kpc float
Min projected distance of the galaxies (kpc)
catalog str
Exact version of the catalog
---
Return
glade_select list
Selected galaxies, all info
sep_select list
Separations for the selected galaxies (arcsec)
dist_kpc_select list
Separations for the selected galaxies (kpc)
"""
# Create a SkyCoord object
coord = SkyCoord(ra=ra, dec=dec,unit=(u.deg, u.deg))
# Query GLADE via Vizier
glade = Vizier.query_region(coord, radius=rad*u.deg,
catalog="glade")
# Empty result
if len(glade) == 0:
return None, None, None
# Check if any of the galaxies found are within an
# acceptable distance
glade_select = []
dist_kpc_select = []
sep_select = []
for g in glade[0]:
# Ignore galaxies too nearby
if g['Dist'] < dist_min or g['Dist'] > dist_max:
continue
sep = SkyCoord(ra=g['RAJ2000']*u.deg,
dec=g['DEJ2000']*u.deg).separation(coord)
# Projected distance (kpc)
dist_kpc = g['Dist']*(10**3)*np.sin(sep)/np.cos(sep)
if dist_kpc >= 0 and dist_kpc < sep_max_kpc:
glade_select.append(g)
dist_kpc_select.append(dist_kpc)
sep_select.append(sep)
# No selected galaxies
if len(glade_select) == 0:
return None, None, None
else:
return glade_select, sep_select, dist_kpc_select
def query_6df(ra, dec, rad=1, dist_min=0, dist_max=np.inf,
sep_max_kpc=np.inf, catalog='VII/259/spectra'):
"""
Fetch galaxies from the 6df catalog
---
Parameters
ra float
Right ascension (deg)
dec float
Declination (deg)
rad float
Search radius (deg)
dist_min float
Min distance of the galaxies (Mpc)
dist_max float
Max distance of the galaxies (Mpc)
sep_max_kpc float
Min projected distance of the galaxies (kpc)
catalog str
Exact version of the catalog
---
Return
gal_select list
Selected galaxies, all info
sep_select list
Separations for the selected galaxies (arcsec)
dist_kpc_select list
Separations for the selected galaxies (kpc)
"""
# Create a SkyCoord object
coord = SkyCoord(ra=ra, dec=dec,unit=(u.deg, u.deg))
# Query GLADE via Vizier
gal = Vizier.query_region(coord, radius=rad*u.deg,
catalog=catalog)
# Empty result
if len(gal) == 0:
return None, None, None
# Check if any of the galaxies found are within an
# acceptable distance
gal_select = []
dist_kpc_select = []
sep_select = []
# Get distances from redshifts
print("Note: distances are luminosity distances (Mpc) \
using Planck+15 cosmology")
for g in gal[0]:
# Ignore galaxies too nearby
dist_g = cosmo.luminosity_distance(g['z']).value
if dist_g < dist_min or dist_g > dist_max:
continue
ra_deg = Angle(g['RAJ2000'] + " hours").deg
dec_deg = Angle(g['DEJ2000'] + " degrees").deg
sep = SkyCoord(ra=ra_deg*u.deg,
dec=dec_deg*u.deg).separation(coord)
# Projected distance (kpc)
dist_kpc = dist_g*(10**3)*np.sin(sep)/np.cos(sep)
if dist_kpc >= 0 and dist_kpc < sep_max_kpc:
gal_select.append(g)
dist_kpc_select.append(dist_kpc)
sep_select.append(sep)
# No selected galaxies
if len(gal_select) == 0:
return None, None, None
else:
return gal_select, sep_select, dist_kpc_select
def query_generic(ra, dec, rad=1, catalog=None):
"""
Fetch galaxies from a generic input catalog
---
Parameters
ra float
Right ascension (deg)
dec float
Declination (deg)
rad float
Search radius (deg)
catalog str
Exact version of the catalog
---
Return
gal_select list
Selected galaxies, all info
"""
# Create a SkyCoord object
coord = SkyCoord(ra=ra, dec=dec,unit=(u.deg, u.deg))
# Query GLADE via Vizier
gal = Vizier.query_region(coord, radius=rad*u.deg,
catalog=catalog)
# Empty result
if len(gal) == 0:
return None
else:
return gal
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description='Query galaxy catalogs')
parser.add_argument('--ra', dest='ra', type=float,
required=True,
help='Right Ascension of the center of \
the query (deg)')
parser.add_argument('--dec', dest='dec', type=float,
required=True,
help='Declination of the center of the \
query (deg)')
parser.add_argument('--r', dest='rad', type=float,
required=False,
help='Radius of the query (deg), \
default=1deg', default=1.)
parser.add_argument('--c', dest='catalog', type=str,
required=False, default="glade",
help='Catalog name. \
The default is GLADE2.3 (VII/281/glade2). \
Available catalogs with distances & separations: \n \
GLADE v2.3 (VII/281/glade2); \
6dF DR3 spec (VII/259/spectra). \
Other catalogs will not have calculated separations. \n \
Some examples are: \
2MASS extended sources (VII/233/xsc); \
HYPERLEDA (VII/237/pgc); \
SDSS DR12 (V/147/sdss12); \
GAIA S/G class (VII/285/gdr2ext)')
parser.add_argument('--dist-min', dest='dist_min', type=float,
required=False, default=0,
help='Minimum distance (Mpc)')
parser.add_argument('--dist-max', dest='dist_max', type=float,
required=False, default=np.inf,
help='Maximum distance (Mpc)')
parser.add_argument('--sep-max', dest='sep_max_kpc', type=float,
required=False, default=np.inf,
help='Maximum projected separation (kpc)')
parser.add_argument('--out', dest='out', type=str,
required=False, default=None,
help='Output file name (CSV)')
args = parser.parse_args()
# Query the GLADE catalog
Vizier.ROW_LIMIT = -1
print("Querying the GLADE catalog")
print(f"Central coordinates: RA, Dec = {args.ra}, {args.dec} (deg)")
print(f"Search radius: {args.rad} deg")
print(f"Minimum distance: {args.dist_min} Mpc")
print(f"Maximum distance: {args.dist_max} Mpc")
print(f"Maximum projected separation: {args.sep_max_kpc} kpc")
print("")
if args.catalog == 'VII/281/glade2':
galaxies, sep, dist_kpc = query_glade(args.ra, args.dec,
rad=args.rad,
dist_min=args.dist_min,
dist_max=args.dist_max,
sep_max_kpc=args.sep_max_kpc,
catalog=args.catalog)
elif args.catalog == 'VII/259/spectra':
galaxies, sep, dist_kpc = query_6df(args.ra, args.dec,
rad=args.rad,
dist_min=args.dist_min,
dist_max=args.dist_max,
sep_max_kpc=args.sep_max_kpc,
catalog=args.catalog)
else:
galaxies = query_generic(args.ra, args.dec,
rad=args.rad, catalog=args.catalog)
if galaxies is not None:
galaxies = galaxies[0]
if galaxies is None:
print("No galaxies found with the given search parameters.")
print("(double-check: does the input catalog exist in Vizier?)")
exit()
else:
print(f"Found {len(galaxies)} galaxies")
print("")
# Create a table with the results
tbl = vstack(galaxies)
if args.catalog == 'VII/281/glade2' or args.catalog == 'VII/259/spectra':
tbl["sep_arcsec"] = sep
tbl["dist_kpc"] = dist_kpc
print(tbl)
if args.out is None:
out = f"galaxies_{str(args.ra)}_{str(args.dec)}_{args.catalog.split('/')[-1]}.csv"
else:
out = args.out
print(f"The results were written in CSV format in {out}")
tbl.write(out, format='csv', overwrite=True)