pylal.galaxyutils

1 # 2 # Copyright (C) 2007-2010 Nickolas Fotopoulos, Larry Price 3 # 4 # This program is free software; you can redistribute it and/or modify it 5 # under the terms of the GNU General Public License as published by the 6 # Free Software Foundation; either version 2 of the License, or (at your 7 # option) any later version. 8 # 9 # This program is distributed in the hope that it will be useful, but 10 # WITHOUT ANY WARRANTY; without even the implied warranty of 11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 12 # Public License for more details. 13 # 14 # You should have received a copy of the GNU General Public License along 15 # with this program; if not, write to the Free Software Foundation, Inc., 16 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 17 18 from __future__ import division 19 20 __author__ = "Nickolas Fotopoulos <nvf@gravity.phys.uwm.edu>" 21 22 import itertools 23 import math 24 25 import numpy 26 import warnings 27 28 ############################################################################## 29 # math 30 ############################################################################## 31

32 -def _make_within_pi(angles):

33 """ 34 Make an array of angles completely between -pi and pi. 35 """ 36 return angles - numpy.sign(angles) * \ 37 numpy.round(abs(angles) / (2 * numpy.pi)) * 2 * numpy.pi

38

39 -def is_inside_polygon(point, vertices):

40 """ 41 Return True if the (2-D) point is inside the (2-D) polygon defined by the 42 vertices. 43 44 Warning: Result is undefined for points lying along the edge of the polygon. 45 46 Adapted from: 47 http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/ (solution 2) 48 """ 49 point = numpy.array(point) 50 centered_vertices = numpy.empty(shape=(len(vertices) + 1, 2), dtype=float) 51 centered_vertices[:-1] = vertices - point 52 centered_vertices[-1] = centered_vertices[0] 53 angles = numpy.arctan2(centered_vertices[:, 1], centered_vertices[:, 0]) 54 diff_angles = _make_within_pi(angles[1:] - angles[:-1]) 55 angle_sum = diff_angles.sum() 56 return abs(angle_sum) >= numpy.pi

57

58 -def is_within_distances(gal_dist, dmin, dmax):

59 """ 60 Return True if dmax > gal_dist > dmin 61 """ 62 if (dmin > dmax): 63 raise ValueError, "minimum distance is greater than maximum distance "\ 64 + str(dmin) + " > "+ str(dmax) 65 if (0 > dmin) or (0 > dmax): 66 raise ValueError, "negative distance " + str(dmin) + " " + str(dmax) 67 return (dmax > gal_dist) and (dmin < gal_dist)

68 69 70 ############################################################################## 71 # unit conversion 72 ############################################################################## 73

74 -def hms2rad(ra_sex):

75 """ 76 Convert right ascension and from a sexagesimal string to floating point 77 radians. Handles h:m:s or h:m. 78 """ 79 tup = ra_sex.split(":") 80 h = int(tup[0]) 81 m = int(tup[1]) 82 s = float(tup[2]) 83 84 if (h < 0 or h > 23) or (m < 0 or m > 60) or (s < 0 or s >= 60): 85 raise ValueError, "hour, minute, or second out of bounds " + ra_sex 86 return 2 * numpy.pi * (h + (m + s / 60) / 60) / 24

87

88 -def dms2rad(dec_sex):

89 """ 90 Convert declination from a colon-delimited sexagesimal string to floating 91 point radians. Handles d:m:s. 92 """ 93 tup = dec_sex.split(":") 94 if tup[0].startswith("-"): 95 d = int(tup[0][1:]) 96 sign = -1 97 else: 98 d = int(tup[0]) 99 sign = +1 100 m = int(tup[1]) 101 s = float(tup[2]) 102 103 if (d > 89) or (m < 0 or m > 60) or (s < 0 or s >= 60): 104 raise ValueError, "degree, minute, or second out of bounds: " + dec_sex 105 return numpy.pi * sign * (d + (m + s / 60) / 60) / 180

106

107 -def hm2rad(ra_sex):

108 """ 109 Convert right ascension and from a sexagesimal string to floating point 110 radians. Handles h:m. 111 """ 112 tup = ra_sex.split(":") 113 h = int(tup[0]) 114 m = float(tup[1]) 115 116 if (h < 0 or h > 23) or (m < 0 or m > 60): 117 raise ValueError, "hour or minute out of bounds " + ra_sex 118 return 2 * numpy.pi * (h + m / 60) / 24

119

120 -def dm2rad(dec_sex):

121 """ 122 Convert declination from a colon-delimited sexagesimal string to floating 123 point radians. Handles d:m. 124 """ 125 tup = dec_sex.split(":") 126 if tup[0].startswith("-"): 127 d = int(tup[0][1:]) 128 sign = -1 129 else: 130 d = int(tup[0]) 131 sign = 1 132 m = float(tup[1]) 133 134 135 if (d > 89) or (m < 0 or m > 60): 136 raise ValueError, "degree or minute out of bounds: " + dec_sex 137 return numpy.pi * sign * (d + m / 60) / 180

138

139 -def amin2rad_or_tilde(amins):

140 """ 141 convert arcminutes to radians 142 """ 143 if amins == '~': 144 return numpy.nan 145 else: 146 return 10800*float(amins)/numpy.pi

147

148 -def deg2rad_or_tilde(degs):

149 """ 150 convert degrees to radians 151 """ 152 if degs == '~': 153 return numpy.nan 154 else: 155 return float(degs)*numpy.pi/180

156

157 -def h2rad(hours):

158 """ 159 convert hours to radians 160 """ 161 return float(hours)*numpy.pi/12

162

163 -def float_or_tilde(num):

164 """ 165 deal with the use of tildes (and other strings) for unknown float quantities 166 in the GWGC catalog 167 """ 168 if num == '~': 169 return numpy.nan 170 else: 171 return float(num)

172

173 -def int_or_tilde(num):

174 """ 175 deal with the use of tildes for unknown int quantities 176 in the GWGC catalog 177 """ 178 if num == '~': 179 return None 180 else: 181 return int(num)

182 183 184 ############################################################################## 185 # galaxy and galaxy catalog representations 186 ############################################################################## 187

188 -class CBCGC(list):

189 """ 190 class for working with the galaxy catalog created and maintained by the 191 CBC group 192 193 Current catalog: 194 http://www.lsc-group.phys.uwm.edu/cgit/lalsuite/plain/lalapps/src/inspiral/inspsrcs100Mpc.errors 195 196 Literature reference: 197 http://arxiv.org/pdf/0706.1283 198 """ 199 valid_columns = { 200 "name": (0, str), 201 "ra": (1, hm2rad), 202 "dec": (2, dm2rad), 203 # distance measured in kpc 204 "distance_kpc": (3, float), 205 # luminosity measured in milky way equivalent galaxies 206 "luminosity_mwe": (4, float), 207 "metal_correction": (5, float), 208 "magnitude_error": (6, float), 209 "distance_error": (7, float), 210 } 211

212 - def entry_from_line(cls, line, load_columns):

213 # create blank entry 214 row = cls.entry_class() 215 216 # parse line 217 tup = line.split() 218 219 # fill the entry 220 for col_name in load_columns: 221 col_index, col_type = cls.valid_columns[col_name] 222 setattr(row, col_name, col_type(tup[col_index])) 223 224 return row

225 entry_from_line = classmethod(entry_from_line) 226

227 - def from_file(cls, fileobj, load_columns=None):

228 # set/validate columns to load 229 if load_columns is None: 230 load_columns = cls.valid_columns.keys() 231 else: 232 for col in load_columns: 233 if col not in cls.valid_columns: 234 raise ValueError, "no such column exists" 235 cls.entry_class.__slots__ = load_columns 236 237 # load them 238 return cls([cls.entry_from_line(line, load_columns) for line \ 239 in fileobj if not line.startswith("#")])

240 from_file = classmethod(from_file) 241

242 - def within_polygon(self, vertices):

243 return self.__class__([gal for gal in self if \ 244 is_inside_polygon(gal.coords, vertices)])

245

246 - def within_distances(self, dmin, dmax):

247 return self.__class__([gal for gal in self if \ 248 is_within_distances(gal.distance_kpc, dmin, dmax)])

249

250 - def __repr__(self):

251 return "\n".join(itertools.imap(str, self))

252

253 -class GalaxyCatalog(CBCGC):

254 """ 255 left here to maintain compatibility with existing codes 256 """

257 - def __init__(self,*args,**kwargs):

258 warnings.warn( \ 259 'The GalaxyCatalog class has been replaced by the CBCGC class.', 260 DeprecationWarning, stacklevel=3) 261 CBCGC.__init__(self,*args,**kwargs)

262

263 -class GWGC(CBCGC):

264 """ 265 useful class for dealing with the gravitational wave galaxy catalog 266 267 Current catalog: 268 https://www.lsc-group.phys.uwm.edu/cgi-bin/pcvs/viewcvs.cgi/bursts/collabs/DO_proposal/gwgc/GWGCCatalog.txt?rev=1.4&content-type=text/vnd.viewcvs-markup 269 """ 270 valid_columns = { 271 #hyperleda identifier 272 "pgc": (0,int_or_tilde), 273 "name": (1, str), 274 "ra": (2, h2rad), 275 "dec": (3, deg2rad_or_tilde), 276 #morphological type 277 "mtype": (4, str), 278 #apparent blue magnitude 279 "app_mag": (5, float_or_tilde), 280 #major diameter 281 "maj_diam": (6, amin2rad_or_tilde), 282 #error in major diameter 283 "maj_diam_error": (7, amin2rad_or_tilde), 284 #minor diameter 285 "min_diam": (8, amin2rad_or_tilde), 286 #error in minor diameter 287 "min_diam_error": (9, amin2rad_or_tilde), 288 #ratio of minor to major diameters 289 "ratio_diams": (10, float_or_tilde), 290 #error in ratio of diameters 291 "ratio_diams_error": (11, float_or_tilde), 292 #position angle of galaxy 293 "pos_ang": (12, deg2rad_or_tilde), 294 #absolute blue magnitude 295 "abs_mag": (13, float_or_tilde), 296 # distance measured in mpc 297 "distance_mpc": (14, float_or_tilde), 298 #distance error in mpc 299 "distance_error": (15, float_or_tilde), 300 #apparent magnitude error 301 "app_mag_error": (16, float_or_tilde), 302 #absolute magnitude error 303 "abs_mag_error": (17, float_or_tilde) 304 } 305

306 - def entry_from_line(cls, line, load_columns):

307 # create blank entry 308 row = cls.entry_class() 309 310 # parse line 311 tup = line.split('|') 312 313 # fill the entry 314 for col_name in load_columns: 315 col_index, col_type = cls.valid_columns[col_name] 316 setattr(row, col_name, col_type(tup[col_index])) 317 return row

318 entry_from_line = classmethod(entry_from_line) 319

320 - def from_file(cls, fileobj, load_columns=None):

321 # set/validate columns to load 322 if load_columns is None: 323 load_columns = cls.valid_columns.keys() 324 else: 325 for col in load_columns: 326 if col not in cls.valid_columns: 327 raise ValueError, "no such column exists" 328 cls.entry_class.__slots__ = load_columns 329 330 # load them 331 return cls([cls.entry_from_line(line, load_columns) for line \ 332 in fileobj if not line.startswith("#")])

333 from_file = classmethod(from_file) 334

335 - def within_distances(self, dmin, dmax):

336 return self.__class__([gal for gal in self if \ 337 is_within_distances(gal.distance_mpc, dmin, dmax)])

338

339 -class CBCGGalaxy(object):

340 """ 341 A galaxy object that knows how to initialize itself from a line in a text 342 file and consumes a minimum of memory. 343 """ 344 __slots__ = GalaxyCatalog.valid_columns.keys() 345

346 - def __str__(self):

347 return "\t".join([str(getattr(self, slot)) for slot in self.__slots__])

348

349 - def __repr__(self):

350 return "Galaxy(\"" + str(self) + "\")"

351

352 - def _coords_getter(self):

353 return (self.ra, self.dec)

354 coords = property(fget=_coords_getter)

355

356 -class GWGCgalaxy(CBCGGalaxy):

357 __slots__ = GWGC.valid_columns.keys()

358 359 GalaxyCatalog.entry_class = CBCGGalaxy 360 CBCGC.entry_class = CBCGGalaxy 361 GWGC.entry_class = GWGCgalaxy 362

Source Code for Module pylal.galaxyutils